Experimental Channel Statistics of Drone-to-Ground Retro-Reflected FSO Links With Fine-Tracking Systems

The utilization of drones as flying base stations (BSs) in the sixth-generation (6G) cellular networks has attracted much attention recently. In this context, free-space optical (FSO) systems could be deployed to provide high-capacity fronthaul/backhaul links between drones and ground BSs. Particularly, a drone-to-ground FSO communication link can be established by equipping the drone with a modulating retro-reflector (MRR) to modulate the incoming optical beam and reflect the modulated beam back along the same path. This helps to alleviate stringent pointing requirements from the drone while satisfying the limited size, weight, and power (SWaP) consumption requirements. Nevertheless, the underlying physical channel effects of this aerial retro-reflected FSO system have never been experimentally explored in the literature. In this paper, we report, for the first time, the experimental channel statistics of a drone-to-ground retro-reflected FSO link, offering practical insights into the angle-of-arrival (AoA) fluctuations at the receiver, the channel coherence time, the probability of fade, the level crossing rate, the average fade duration, and the time-frequency channel characteristics. Our results are expected to serve as practical sources of reference for the theoretical performance analyses and engineering designs of drone-based retro-reflected FSO systems.


I. INTRODUCTION
The integration of unmanned aerial vehicles (UAVs) into the sixth-generation (6G) networks has recently become a hot topic, due to the rapid advancement of UAV technologies [1]. When multiple UAVs are deployed as flying base stations (BSs), they can either operate independently or cooperatively in a swarm to serve users simultaneously or share the data load with each other through inter-UAV links. Of particular interest at the very-low-altitude platform (VLAP), e.g.
≤50 m, is the utilization of small-size UAVs, e.g. multirotor drones, serving as flying BSs to provide on-demand all-location communication services to ground users while maintaining high-capacity aerial fronthaul/backhaul links to The associate editor coordinating the review of this manuscript and approving it for publication was Miguel López-Benítez . ground BSs. With the heavily congested spectrum in radio frequency (RF), one promising candidate for the aerial fronthaul/backhaul links is free-space optics (FSO), due to its license-free spectrum and ability to deliver high-capacity point-to-point communications [2], [3]. In addition, the FSO links can be utilized to simultaneously transmit the energy to charge the batteries for drone-based BSs, thereby increasing their service time [4]. These aerial FSO links are also paving the way for the employment of free-space quantum key distribution (QKD) protocols to offer the information-theoretic security between terrestrial and aerial platforms in the future quantum wireless networks [5]- [7].
Terrestrial FSO systems are susceptible to misalignments due to transceiver sway, platform vibrations, and beam wander caused by the large-scale eddies in the atmosphere. The impairment of misalignments is even more severe in drone-based FSO systems due to the random fluctuations in the position and orientation of the drone, which may result in the degradation or total loss of the received signal. Therefore, acquisition, tracking, and pointing (ATP) mechanisms stand out as important solutions to reduce the impact of misalignments in mobile platforms [8]. However, the ATP subsystem is typically large, heavy, and high power-consumption, which may not be suitable for the limited size, weight, and power (SWaP) consumption aboard the drone. In addition, the ATP unit may require either one [6] or multiple laser beacons [9] at different wavelengths with the communication laser to be installed on the drone, which is sometimes daunting from the practical implementation perspective.
To relax the SWaP constraint and pointing accuracy requirements, an alternative solution is to replace the ATP unit on the drone with a modulating retro-reflector (MRR), which is assembled by a light modulator, i.e. to modulate the incoming optical beam, and a passive retro-reflector, i.e. to reflect the modulated beam back along its incident path [10]. As a result, the sophisticated ATP subsystem can be moved to the optical transceiver on the ground. However, current MRRs are based on an electro-optic modulation which results in a low data rate, e.g. up to Mbps [11], [12], due to the fundamental bandwidth bottleneck between optics and electronics. Fortunately, in the quest to achieve Gbps data rates, an all-optical retro-modulation where the retro-reflected optical signal is modulated all-optically was proposed [13]. In addition, a retro-reflected FSO link operating at 500 Mbps at a range of 560 m between an octocopter drone and a ground station has been recently reported [14], with the potential of achieving 1∼2 Gbps by using electro-absorption modulators (EAMs). These advancements promise the utilization of MRRs for the future drone-based high-capacity fronthaul/backhaul links.
From the theoretical aspect, the statistical channel models for the retro-reflected terrestrial FSO links have only been investigated by a limited number of research works [15]- [23]. In fact, all previous works merely focused on characterizing the effect of atmospheric turbulence [15]- [21], and the misalignments in the forward link (i.e. from the transmitter to the MRR) [22], [23], without considering the misalignments in the retro-reflected link when the MRR is installed on a hovering drone. Furthermore, from the experimental aspect, the recent experiment in [14] only concentrated on the achievable data rate and bit-error rate (BER). To the best of authors' knowledge, the underlying physical channel characteristics of the atmospheric turbulence and misalignments-induced fading in the drone-based retro-reflected FSO systems have never been experimentally investigated in the literature.
Motivated by the above-mentioned facts, in this paper, we report an optical beam propagation experiment between an optical ground station (OGS) and a hexacopter drone equipped with a passive corner-cube retro-reflector (CCR), over a maximum 204-m roundtrip link. The novel contributions of this paper can be summarized as follows.
• Although the CCR can help to alleviate the pointing requirements, the randomly hovering nature of the drone would still result in severe angle-of-arrival (AoA) fluctuations at the OGS. For high data-rate systems where small detectors or fiber-coupled receivers are required, the AoA fluctuations could possibly cause the received optical beam to fall out of the detector/fiber-tip's field of view (FoV). To mitigate these AoA fluctuations, we present the practical design of a compact (30 cm width × 30 cm length) fine-tracking system at the OGS, which is composed of a fast-steering mirror (FSM) and a quadrant detector (QD) working in a closed-loop operation governed by a proportional-integral-derivative (PID) controller. The tracking accuracies 1 of less than 6 µm and 9 µm are then achieved with the drone's hovering altitudes at 15 m and 20 m, over roundtrip distances of 201 m and 204 m, respectively.
• The drone's random hovering trajectory is, for the first time, experimentally visualized from the AoA fluctuations data at the OGS. This also enables the estimation of the retro-reflected beam-centroid displacements, which can be subsequently fitted to well-known statistical probability density functions (PDFs) characterizing the misalignments. From this, separated effects of the misalignments and the atmospheric turbulence can be comprehensively investigated by fitting the received power histogram data with the corresponding theoretical composite channel PDFs.
• Finally, important insights into the drone-based retroreflected physical channel are analyzed for the first time in the literature, thus revealing the channel coherence time, the probability of fade, the level crossing rate (LCR), the average fade duration (AFD), and the time-frequency channel characteristics. These practical results are also compared with theoretical backgrounds, providing crucial references for the theoretical analyses and engineering designs of the drone-based retroreflected FSO systems. The remainder of this paper is organized as follows. Section II presents the theoretical channel PDFs characterizing the atmospheric turbulence and misalignments-induced fading that will be used to fit with the experimental data. Section III describes the related physical channel characteristics. The practical design of the fine-tracking system at the OGS is detailed in Section IV, followed by the actual experimental setup and results reported in Section V. Finally, the conclusion is drawn in Section VI.

A. ATMOSPHERIC ATTENUATION
The molecular absorption and aerosol scattering suspended in the atmosphere can cause a partial loss in the received 1 The tracking accuracy is defined as the standard deviations of the focused beam-centroid displacements in the horizontal and vertical axes at the receiving plane of a small detector or a fiber tip. VOLUME 9, 2021 signal, which is known as the atmospheric attenuation, given by the Beer's law as h l = exp(−β l L), where β l denotes the attenuation coefficient and L is the optical beam propagation distance [24]. For the retro-reflected link, since the optical beam is propagated twice in the same path, the atmospheric attenuation is simply derived as h l = exp(−β l (2L)).

B. ATMOSPHERIC TURBULENCE-INDUCED FADING
Atmospheric turbulence is an effect caused by the refractive-index variations along the propagation path, which results in received power fluctuations at the receiver, i.e. scintillation or fading [24]. Since the optical beam experiences the atmospheric turbulence of the same path in both the forward and retro-reflected links, it has been proved by experimental measurements [15]- [17] and wave-optics simulations [18] that their fading channels are statistically correlated and the retro-reflected channel should be modeled as the product of two correlated random variables (RVs) [19], [20]. If the transmitter and receiver are co-located at the OGS, i.e. monostatic configuration, the forward and retro-reflected links show a high correlation. On the contrary, if the transmitter and receiver are separated by the order of a Fresnel zone radius, i.e. r F = √ Lλ/2π with λ the optical wavelength, the forward and retro-reflected links can be considered independent. This is also known as the bistatic configuration, which is chosen in our OGS design. However, as the optical beam is widely diverged in our experiment, the overlap of the beams can make the off-axis receiver appear to experience the same correlation effects of a monostatic system [16].
As a result, the turbulent channel coefficient can be written as h t c = h t 1 h t 2 , where h t 1 and h t 2 are the channel coefficients of the forward and retro-reflected links, respectively. To quantify the strength of the atmospheric turbulence, the scintillation index (SI) is usually used and defined by [24] where E[·] denotes the statistical expectation. To ensure that fading does not attenuate or amplify the average power, the fading coefficient is normalized, i.e. E h t c = 1. In general, weak and moderate-to-strong turbulence conditions are characterized by SI < 1 and SI ≥ 1, respectively.

1) LOG-NORMAL TURBULENCE MODEL
For weak turbulence conditions, the fading channel coefficient can be modeled as h t c = exp(2(X 1 + X 2 )) = exp(2X c ) [19], where X 1 and X 2 are respectively the log-amplitudes of the optical intensity in the forward and retro-reflected paths following Gaussian distributions, and X c is the sum of two correlated Gaussian RVs. Hence, the PDF of h t c follows a log-normal (LN) distribution, written as [19] f LN (h t c ) where µ X c and σ X c are the mean and variance of X c . The SI under weak turbulence conditions modeled by the LN distribution can be calculated as [24] SI LN = exp 4σ 2 X c − 1. (3)

2) GAMMA-GAMMA TURBULENCE MODEL
For moderate-to-strong turbulence conditions, the turbulence fading coefficient can be modeled by a Gamma-Gamma (GG) distribution [25]. The retro-reflected turbulence fading coefficient thus can be described as the product of two correlated GG RVs. Although an expression for the joint PDF of the product of two correlated GG RVs was readily derived in [21], it involves complicated infinite summation terms. Nevertheless, it has been verified by experiments in [15] that a single GG distribution could provide a well approximation to the histogram data of the retro-reflected link. Thus, the PDF of h t c can be given as [25] f where (·) denotes the Gamma function defined as (w) = ∞ 0 t w−1 e (−t) dt, K α−β (·) represents the modified Bessel function of the second kind of order α − β ; α > 0 and β > 0 are the effective numbers of large-scale and smallscale atmospheric eddies. The SI under moderate-to-strong turbulence conditions modeled by the GG distribution can be given as [25]

3) α-µ TURBULENCE MODEL
For moderate-to-strong turbulence conditions modeled by a GG distribution, an alternative approach to approximate the product of two correlated GG RVs for characterizing the retro-reflected FSO channel was proposed in [20]. Specifically, the product of two correlated GG RVs is approximated to an α-µ RV using the moment matching method. Thus, the PDF of the channel coefficient can be expressed as [20] where α > 0 is an arbitrary parameter,ĥ t c is an α-root mean value given asĥ t c = α E h t c α , and µ > 0 is the inverse of the normalized variance of h t c α expressed operator. From (1) and [20, (10)], the SI under moderate-tostrong turbulence conditions modeled by the α-µ distribution can be derived as C.

MISALIGNMENTS-INDUCED FADING
In our experiment, the drone carrying the CCR hovers with random trajectories within the static beam footprint in the forward link, hence leading to random movements of the retro-reflected beam with respect to the OGS's receiving telescope. Since the misalignments in the forward link result in the misalignments in the retro-reflected link, the misalignments-induced fading coefficient can be considered as the product of two correlated RVs, written as h p c = h p 1 h p 2 , where h p 1 and h p 2 are two correlated RVs characterizing the misalignments in the forward and retroreflected links, respectively. The misalignments-induced fading coefficient in each link can be readily given as [26], where A 1 and A 2 are the fractions of power falling within the CCR in the forward link and the receiving aperture in the retro-reflected link, respectively, when there are no misalignments; r 1 and r 2 are respectively the displacements of the drone's position with respect to the forward beam centroid and the retro-reflected beam centroid with respect to the center of the receiving aperture; w zeq,1 and w zeq,2 are the equivalent beam-widths of the forward and retro-reflected beams, respectively. The random displacements in both the horizontal x and vertical y directions are commonly modeled as independent Gaussian RVs, with means {µ x , µ y } and variances {σ x , σ y }. In the most general case, both displacements are nonzero mean Gaussian RVs with different jitters, i.e. (µ x = µ y , σ x = σ y ), thus r i follows a four-parameter Beckmann distribution written as [27] Depending on the values of {µ x , µ y , σ x , σ y }, (8) reduces to several well-known and tractable distributions such as the , and the Hoyt distribution (µ x = µ y = 0, σ x = σ y ), as detailed in [27]. These distributions correspond to different special cases of the misalignments characterized by the four-parameter Beckmann distribution. Although the closed-form expression for the PDF of correlated misalignments-induced fading channels was derived in [28] for the case of Rayleigh-distributed displacements, deriving the PDF for the coefficient h p c as the product of two correlated RVs for all cases is generally complicated, if not impossible. For the sake of simplicity, we utilize the misalignments-induced fading PDF of a single RV as an approximation approach for fitting with the experimental histogram data, i.e. r i = r. Since the four-parameter Beckmann distribution can be conveniently approximated by a Rayleigh distribution with a modified variance [29], the PDF of radial displacements can now be written as where is the approximated jitter variance. As a result, the PDF of h p c can be derived as [29] where ϕ mod = w zeq 2σ mod is the ratio between the equivalent beam radius and the jitter variance, 2σ y are the ratios between the equivalent beam radius and the jitter variances in the x and y directions, respectively. To quantify the strength of the fading caused by misalignments, the definition of SI in (1) is also applied. With the help of [27, (23)], the SI for misalignments can be newly derived as As the misalignments-induced fading becomes more severe when ϕ mod → 1 and weaker when ϕ mod > 5, we define that weak, moderate, and strong misalignments correspond to SI p c ≤ 0.001, 0.001 < SI p c < 0.3, and SI p c ≥ 0.3, respectively.

D. AoA FLUCTUATIONS-INDUCED FADING
Due to drone hovering trajectories, the incident angle of the optical beam at the receiver experiences severe fluctuations, leading to jitters on the focal plane of the receiving lens that can fluctuate or interrupt the received power at a detector or a fiber tip with limited FoVs [30]. The total radial AoA deviations can be written as θ AoA = θ 2 X + θ 2 Y , where θ X and θ Y are the AoA fluctuations in the horizontal and vertical directions, respectively. θ X and θ Y can be modeled as zero-mean independent Gaussian RVs with variances σ θ X and σ θ Y , respectively. The corresponding loss due to AoA fluctuations, denoted as h AoA , can take two discrete values ''1'' and ''0'' to determine if the incident beam spot falls within the detector's FoV or not [31]. This can be formulated where θ FoV represents the detector's half-angle FoV. This assumption of h AoA only holds if the power contributions from the side lobes in the Airy pattern are ignored, since some fractions of power could still be detected even when θ AoA > θ FoV . In our experiment, the AoA fluctuations at the OGS are within the tracking detector's FoV, thus enabling the fine-tracking system to compensate for the AoA fluctuations, making the focused beam spot aligned to the center of the detector, and the whole beam spot size is kept within the photosensitive area, resulting VOLUME 9, 2021 in h AoA = 1. Consequently, we could exclude the AoA fluctuations-induced fading in our statistical analysis.

E. COMPOSITE CHANNEL MODELS
The composite channel coefficient of the double-pass retroreflected link can be expressed as h = Rh l h t c h p c , where R is the normalized reflection coefficient of the CCR, which can be considered constant within a certain transmission duration [19].

1) LN AND APPROXIMATED BECKMANN COMPOSITE CHANNEL
From (2) and (10), the composite channel PDF of the retro-reflected link can be expressed as [26] where x exp −t 2 dt is the complementary error function, and µ = 2σ 2 X c 1 + 2ϕ 2 mod . From (1) and given the fact that h t c and h p c are statistically independent processes [27], the SI of the composite fading channel can be calculated as With the help of [27, (10) and (23)], (13) can be newly derived as 2) GG AND APPROXIMATED BECKMANN COMPOSITE CHANNEL From (4) and (10), the composite channel PDF of the retro-reflected link can be expressed as [32] where G m,n p,q [·] represents the Meijer's G-function [33, (9.301)]. From (13) and [27, (14) and (23)], the SI of the composite fading channel can be newly derived as
Remark 2: For statistical verifications, both the data and the composite channel PDFs, i.e. in (12), (15), and (17), are normalized to the mean. The normalized PDFs then contain only the parameters of the fluctuations that could be fitted with the normalized experimental data. Details on the derivations of normalized composite channel PDFs are presented in Appendix B.

III. PHYSICAL CHANNEL CHARACTERISTICS
A. CHANNEL COHERENCE TIME The channel coherence time, denoted as T c , measures the time duration wherein two signals show high amplitude correlation, which can be directly calculated from the normalized auto-covariance function of the received power, denoted aŝ , by considering the coherence time at 1/e of the peak at zero lag, written asˆ ( being the auto-covariance function, wherê P Rx,n = P Rx,n − 1 N N i=1 P Rx,i with P Rx,i the ith sample of the total N samples of the received power. The channel coherence time reveals the time duration over which the channel fading coefficient is considered to be not varying. This metric is of critical importance for the error-correction code and inter-leaver designs in FSO communication systems [36].

B. RECEIVED POWER SI
The SI, i.e. normalized variance of the received power, can be experimentally derived from the received power data as where The experimental SI P Rx can be later compared with the theoretical SIs derived in (14), (16), and (18).

C. PROBABILITY OF FADE
The experimental probability of fade is defined by the probability that the received power level falls below a given fade threshold P th , written as P (P Rx ≤ P th ) = F P Rx (P th ) , where F P Rx (P Rx ) is the empirical cumulative distribution function (CDF) of the received power data. It is generally customary to investigate the experimental probability of fade versus the fade margin, defined as F T = 10 log 10 E P Rx /P th , F T ≥ 0 dB, to find out the power margin above the fade threshold level to achieve a specific fade probability. The experimental fade probability can be compared with the corresponding theoretical CDF of the composite channel, denoted as F h (h 0 ) with h 0 a threshold level and F h (h 0 ) the CDF of the composite channel. The theoretical F T can be expressed as F T = 10 log 10 E h /h 0 . Theoretical closed-form expressions of F h (h 0 ) for the composite channels in Sections II-E1 and II-E2 can be found in [26, (18)] and [32, (15)], respectively.

D. LEVEL CROSSING RATE
The LCR is defined as the average number of crossings of the signal above a certain threshold h 0 per unit of time, mathematically given as where f h,ḣ (·) denotes the joint PDF of h and its time derivativė h. It is also customary to investigate the LCR versus the normalized threshold level defined as ρ 0 = h 0 /E(h). A theoretical closed-form expression of the LCR for an FSO channel considering the LN turbulence and approximated Beckmann misalignments is newly derived in Appendix C.

E. AVERAGE FADE DURATION
The AFD is defined as the average duration during which the signal stays below a certain threshold h 0 , written as where F h (h 0 ) and L h (h 0 ) are the above-defined CDF (Section III-C) and LCR (Section III-D), respectively. Thus, a theoretical closed-form expression of the AFD considering the LN turbulence and approximated Beckmann displacements can be straightforwardly obtained. The LCR and AFD are important second-order statistical metrics that describe the time-varying behavior of a fading channel [37].

F. TIME-FREQUENCY CHARACTERISTICS 1) POWER SPECTRAL DENSITY
For an FSO link, it is essential to quantitatively analyze the temporal behavior of the received power fluctuations, which is particularly useful for evaluating the burst error rate and the design of optimal schemes for the detection and coding.
Based on Taylor's frozen turbulence hypothesis, which permits the conversion of spatial statistics into temporal statistics, the temporal spectrum of irradiance fluctuations, i.e. power spectral density (PSD) denoted as S (ω), of the optical signal for point receivers is defined by the Fourier transform of the temporal covariance function as where ω is the angular frequency and B (τ, L) is the temporal covariance function [24]. Based on this theory, there exists a cut-off frequency f c that divides the spectrum into two parts, where the part below f c has a relatively flat PSD and the part above f c decays following a power-law exponent κ = 8/3, i.e. the PSD decays by f −κ with f a frequency vector. f c can be linked with the channel coherence time T c as f c = √ 2π /T c . In practice, since non-point receivers are typically used for communication systems, a slight deviation in the power-law exponent of the experimental spectrum from the theoretical one is expected. From the received power data, the PSD can be estimated by its periodogram, i.e. the Fourier transform of the biased estimate of the autocorrelation sequence, given aŝ where −1/2 t < f ≤ 1/2 t and t is the sampling interval.

2) SPECTROGRAM
One effective tool to visualize the spectrum of frequencies of a signal as it varies with time is the spectrogram, defined as an intensity plot of the short-time Fourier transform (STFT) magnitude. The STFT is a sequence of fast Fourier transforms of windowed data segments, where the windows are overlapped in time. Given that the received power data P Rx is a discrete time signal with N samples and consecutive segments of length m, let S ∈ R m×(N −m+1) be the matrix with the consecutive segments as consecutive columns, i.e.
[P Rx [0] , P Rx [1] , . . . , P Rx [m − 1]] T is the first column, and [P Rx [1] , P Rx [2] , . . . , P Rx [m]] T is the second column, and so forth, with the rows and columns indexed by time. The spectrogram of P Rx with window size m is the matrixŜ whose columns are the discrete Fourier transform of the columns of S, given asŜ =FS, whereF is the complex conjugate of the Fourier matrix F expressed as Each location onŜ is a point in frequency, i.e. the rows, and time, i.e. the columns.

IV. FINE-TRACKING SYSTEM DESIGN A. OPTICS DESIGN
The optics configuration at the OGS is shown in Fig. 1a and a 3-dimensional (3D) design using an optical ray-tracing software is illustrated in Fig. 1b. We employ the bistatic configuration, where the transmitting laser and the receiving telescope are separated by 2.5 cm. The transceiver consists of a high-power 976-nm laser and a 5-cm telescope aperture. The output beam is transmitted through a diverging lens with the full divergence angle of about 0.04 rad, which makes a beam footprint of 4 m at 100-m line-of-sight (LoS) distance. 2 This divergent beam is designed to sufficiently illuminate the drone in hovering conditions with some margin. Table 1 summarizes the main information of the OGS's components.
To aid the initial link alignment process, two alignment cameras, i.e. a wide-FoV external camera and a narrow-FoV internal camera, are included. The fine-tracking system consists of an FSM, which is placed at 45 • relative to the telescope's optical axis, and a QD working in a closed-loop operation with the FSM, governed by a PID controller. The optical beam at the QD is defocused so that the FSM can steer the beam through all the QD's FoV. The FSM also has a built-in high-precision optical sensor to monitor the AoA 2 It is noted that the maximum power of 9 W at 976-nm wavelength is employed in our experiment to compensate for the large geometrical losses of the forward and retro-reflected beams. The 976-nm wavelength is selected only for the sake of the experimental setup due to the availability of the high-power source. Nevertheless, when deploying MRR-based FSO communication systems, the wavelengths in C-band (1530∼1565 nm) are preferred and smaller beam divergence angle with lower transmitted power can be applied [14]. fluctuations of the received beam. A simplified diagram of the fine-tracking system is depicted in Fig. 1c, where a power meter (PM) is installed at the focal plane of lens 3 to emulate a possible receiving optical fiber or small detector. All the optics, i.e. lenses and beam splitters (BSs), and other devices are compactly housed on a 30 × 30-cm optical breadboard 3 and mounted on top of a 3-axis gimbal tripod to assist the alignment for the entire terminal.

B. POINTING AND TRACKING
At the initial link alignment process, the coarse pointing is performed by using the external camera to point the entire terminal towards the target drone. Once the link is roughly aligned, the internal camera is used to further tune the tracking of the optical beam captured by the 5-cm telescope. The external and internal cameras respectively provide FoV frames of 6.5 × 4.8 m and 3.27 × 2.45 m at approximately 100-m LoS distance. The wide-FoV camera allows to quickly search for the target drone, while the narrow-FoV camera helps to locate the CCR attached on the drone and align with the OGS's receiving telescope. At this stage, the system is in the open-loop operation, thus the received beam spot is affected by the AoA fluctuations resulting in the image dancing on the focal plane.
The fine-tracking system is based on a control loop to control the FSM to correct the AoA of the received beam. The tracking control loop uses the received beam spot incident onto the QD to detect changes in the beam's AoA and compensates the AoA error with the FSM. More specifically, the QD outputs four signals proportional to the received optical power incident on each segment and compares these signals to determine the centroid position of the beam spot. The intersection of the transition region is considered as the origin of the X -Y coordinates. A PID controller is then used 3 The compact size of the OGS makes it suitable for possible installations on pole-mounted BSs, which could be located on the ground or building rooftops in the cellular networks. to close the feedback loop from the QD to the FSM, making the beam spot centered to the detector's photosensitive area.

C. OPTICAL RAY-TRACING SIMULATION RESULTS
An optical ray-tracing software was used in the design stage of the OGS to characterize the behavior of the received optical beam. 4 Fig. 2a shows the simulation result of the beam spot radius at different positions in the optics configuration. The beam spot radius is defined as the radial distance at which the relative intensity has decreased to 1/ exp(2) or 0.135 of its peak value for a Gaussian beam. It is seen that the values of the spot radius after the telescope, and at the QD and PM are 4.2 mm, 1.83 mm and 0.147 mm, respectively. The internal collimated beam spot radius was designed to be 4.2 mm after the telescope so that it is well within the FSM's surface area with some margin for misalignments. In addition, the spot radius at the QD was designed to be approximately half the size of the QD's photosensitive area, which is a common design choice [38]. Finally, the spot radius at the PM is very small because the PM's active area was placed approximately at the focal plane of lens 3 to emulate a possible beam coupling into an optical fiber or a small detector. The stability of the beam spot position at the PM is guaranteed by putting the FSM's surface at the internal pupil of the receiving telescope. The external FoVs of the QD and PM can then be determined by changing the AoA of the simulated incoming beam at the telescope entrance until the corresponding internal beam centroid is detected across the whole detector size, which respectively are 29.32 and 34.9 milliradian (mrad) as seen in Fig. 2b. This makes an FoV frame of about 2.93 × 2.93 m at a 100-m LoS distance for the QD. There also exists a linear relation between the changes in the centroid position at the QD and PM due to the AoA variations at the telescope entrance, depicted in Fig. 2c. This relation can be used to determine the beam displacements at the PM from the data measured at the QD.

A. EXPERIMENT SETUP AND DATA ACQUISITION
The field experiment was conducted on August 27 th , 2020 at NICT Kashima Space Technology Center in Kashima city, 4 In optical engineering, ray tracing is a technique to emulate the propagation of optical wavefronts through a system consisting of various optical components. In our simulations, 861 optical rays were used to simulate an optical beam with a sufficient accuracy. Ibaraki prefecture, Japan. Fig. 3a shows the experiment field where the OGS was located at 100-m distance away from the drone's position while the drone was hovering at 15-m and 20-m altitudes. 5 The commercial drone used in the experiment was the DJI Matrice 600 Pro [39]. A 3-axis DJI Ronin MX gimbal [40] was attached beneath the drone to carry a 63.5-mm CCR and a camera for alignments with the OGS, as shown in Fig. 3b. At the initial alignment process, the drone's camera is used to find the OGS and align with the transmitting laser by remotely controlling the pointing direction of the drone gimbal from a remote controller at the OGS. It is noted that the gimbal was operated in the independent mode, i.e. the pan, tilt, and roll axes can be adjusted independently with the changes of the drone's orientation due to hovering. As a result, the retro-reflected beam is not affected by the drone's orientation but the changes of the drone's hovering position and gimbal's vibration. Fig. 3c shows an example of the image captured from the drone's camera together with the flight information such as the drone's hovering height above the ground and the horizontal distance with respect to the drone's remote controller at the OGS. Once the link is roughly aligned, the OGS's gimbal is fixed and the drone is allowed to randomly hover without any further control. The OGS was designed so that even when the drone is randomly hovering, it is still within the FoV frames of the external and internal cameras, as illustrated in Fig. 3d. 5 These drone's altitudes are equivalent to ground BSs deployed on the rooftops of 5-story and 6-story buildings, respectively. The hovering altitude is defined as the altitude of the drone at the initial stage of the link alignment. VOLUME 9, 2021 For operating the fine-tracking system, a LabVIEW program was built to calibrate and control the PID closedloop operation. All the measurements were done with the closed-loop tracking and experimental data were obtained by a National Instruments TM data-acquisition card at the sampling rate of 4 kHz. In each dataset examined in Sections V-B and V-C, 12,000 samples over 3 seconds were collected. To compare the experimental histogram data with the theoretical PDFs, we utilize the root mean square error (RMSE) and goodness-of-fit R 2 statistical metrics [41]. RMSE determines how well the considered distribution predicts the experimental data, i.e. the data are well predicted by the distribution when RMSE → 0, while R 2 is used to test the distribution's fitness, i.e. R 2 → 1 determines the best fit. To the best of authors' knowledge, in our paper, it is the first time that the theoretical composite PDFs of a retro-reflected FSO channel are verified with experimental data in the literature.

B. DRONE HOVERING AT 15-m ALTITUDE
In this section, we examine a dataset when the drone was hovering at 15-m altitude, making a slanted path of 101 m between the OGS and the drone. The optical beam then propagated over a total 202-m roundtrip distance. 6 Figure 4a shows the AoA measured by the built-in sensor of the FSM, which is the optical deflection angle of the beam exiting the FSM relative to its neutral position. In the closedloop operation, the QD detects the residual beam deviations and sends the error feedback signal to the PID controller. The controller then sends a control signal to the FSM to correct the AoA so that the beam spot is directed to the center of the QD. This can be clearly seen in Fig. 4b, where the FSM tracking signals show the correction of AoA over time in both axes. Fig. 4c depicts the beam-centroid normalized position density measured at the QD in the closed-loop operation, and from the position relation between the QD and PM in Fig. 2c, the normalized position density is estimated at the PM in Fig. 4d. It is seen that the high density of the beam centroid position appears around the centers of the QD and PM with small displacements, i.e. with standard deviations σ X = 11, σ Y = 8 µm at the QD and σ X = 5.91, σ Y = 4.36 µm at the PM. This confirms the tracking accuracy of less than 6 µm at the PM. From Fig. 4d and the beam spot size at the PM in Fig. 2a, it is suggested that the beam coupling into a very small detector, e.g. diameter of typically less than 1 mm, or a multi-mode fiber (MMF) with a core diameter of 300 µm, should be feasible. In the insets of Figs. 4c and 4d, optical ray-tracing simulation results based on the AoA data from Fig. 4a respectively show the beam-spot positions at the QD and PM if the system is in the open-loop operation, i.e. no fine-tracking and the FSM serves as a static mirror. Without the closed-loop tracking, the beam spot is largely deviated from the center position at (0,0), i.e. red spots in the insets, which confirms that the beam coupling into a small detector or a fiber tip at the PM position is not possible. To further investigate the statistics of AoA fluctuations after fine-tracking corrections in both axes, respectively denoted as { X , Y } and {θ X , θ Y }, at the QD and PM, Figs. 4e-4h depict the AoA histogram data versus the fitted PDF following a Gaussian distribution. Both metrics RMSE and R 2 confirm that the theoretical Gaussian PDFs fit well with the data. The means and variances in Figs. 4e-4h could be used as practical reference values for analyzing the AoA fluctuations-induced fading when coupling the received beam into a small detector or a fiber tip with a fine-tracking system.
The AoA fluctuations in both axes measured at the receiving telescope entrance are plotted in Fig. 5a, where the means (µ X = 6.96 mrad and µ Y = 1.32 mrad) and standard deviations (σ X = 1.47 mrad and σ Y = 1.17 mrad) in both axes are provided to serve as practical references for theoretical analyses of the AoA fluctuations and beam misalignments at the receiving telescope. With the known LoS distance and the AoA data in Fig. 5a, the normalized position density of the drone hovering trajectory can be inferred in Fig. 5b with the means µ p X = 0.7 m and µ p Y = 0.13 m, and standard deviations σ p X = 0.15 m and σ p Y = 0.12 m in both hovering directions, respectively. These mean and standard deviation values are important references for characterizing the hovering-induced misalignments in drone-based FSO systems. Since the optical beam is reflected back along the same path from the drone to the OGS, it can be considered that the AoA at the telescope entrance is equivalent to the displacement angle of the retro-reflected beam from the drone. From this, the radial displacements of the retro-reflected beam centroid with respect to the center of the receiving telescope can be determined and plotted in Fig. 5c.
The coherence time of the beam centroid displacement is 700 milliseconds (ms) as shown in Fig. 5d, which can be derived by applying the definition in Section III-A to the data in Fig. 5c. This confirms the long coherence time of the beam misalignments compared to that of the atmospheric turbulence, typically 10∼100 ms. Theoretical PDFs of misalignments in Section II-C are then fitted to the data in Fig. 5c and the results are plotted in Fig. 5e. Details of the statistical metrics and fitting parameters are summarized in Table 2. It is seen from the RMSE and R 2 metrics that the Beckmann distribution is best fit with the data, followed by the Hoyt distribution, and finally the Rayleigh and Rician distributions. To evaluate the small-scale mismatch at the tails of all the fitted PDFs, a log-log plot is provided in Fig. 5f. It is shown that all the fitted PDFs moderately overestimate the displacements at the lower tail and underestimate the displacements at the upper tail. This discrepancy could result from the approximation approach in Section II-C.
To analyze the combine effects of atmospheric turbulenceand misalignments-induced fading, the total received optical power at the QD shown in Fig. 6a is analyzed. The normalized data are then fitted to theoretical composite PDFs described   in Section II-E, as plotted in Fig. 6b. Details of the statistical metrics and fitting parameters are summarized in Table 3. From the RMSE and R 2 metrics, it is suggested that the composite PDFs in (12) and (15) are better fit to the data, compared to (17). A closer look at the tails of the PDFs is illustrated in a log-log scale in Fig. 6c. It is confirmed that (15) is slightly better than (12) in matching the upper tail with the data, while (12) shows a better fit at the lower tail. Nevertheless, both (12) and (15) marginally overestimate both tails of the data PDF. As seen in Table 3, the theoretical SIs of the composite fading channel are calculated based on (14), (16), and (18). Due to the approximation approach, a small difference of about 2∼3% in the theoretical SI and the measured SI P Rx calculated from (20) can be seen. It should be emphasized that our analysis helps to separately evaluate the SIs of the turbulence and misalignments, highlighting that the retro-reflected FSO link is affected by moderate misalignments and weak turbulence conditions, according to definitions in Sections II-B and II-C.
For further investigations of the physical channel characteristics, the channel coherence time can be directly deduced  from the definition in Section III-A, which is found as T c = 180 ms in Fig. 6d. This confirms the slow composite fading characteristic and suggests that the misalignments-induced fading is the dominant factor. The probability of fade is then investigated versus the fade margin F T in Fig. 6e and compared with the theoretical composite PDFs in (12) and (15). It is seen that the data are well predicted at F T ≤ 1.4 dB and gradually mismatched with the theoretical results at higher F T . From the data, to attain a probability of fade at 10 −4 , a 2.35-dB power margin should be guaranteed, while the theoretical results from (12) and (15) suggest the margins of 3.8 dB and 4 dB, respectively. Finally, the LCR and AFD of the physical channel are investigated in Fig. 6f. The experimental LCR and AFD are obtained by simply counting the average number of crossings above a threshold of the data in Fig. 6a and measuring the average duration during which the data stay below the threshold. As expected, the LCR reaches its peak around ρ 0 = 0 dB and decreases when the threshold is set at lower or higher values. The AFD is, however, proportionally increasing with the values of ρ 0 , as the fade duration becomes longer at higher thresholds. Theoretical results in Sections III-D and III-E are also plotted using fitted parameters in Table 3 to compare with the experimental results and a reasonable accuracy is achieved within 1 dB near ρ 0 = 0 dB.
In addition to the magnitude fluctuations of the signal, the speed of the fluctuations is also an important figure of merit. Particularly, the PSD profiles of the radial AoA fluctuations at the telescope entrance and the normalized received power are illustrated in Figs. 7a and 7b, respectively. Since the AoA fluctuations are mainly the consequence of the drone hovering, it is evident from the PSD in Fig. 7a that significant magnitudes are contained mostly in frequencies below 50 Hz. In Fig. 7b, the cutoff frequency is determined from the coherence time as f c ≈ 13.93 Hz, which separates the part with a relatively flat PSD, i.e. below 13.93 Hz, and the part that decays with a power law exponent κ, i.e. above 13.93 Hz. By applying the linear regression method, it can be estimated that κ = 2.3378, which is close to the theoretical value of 8/3 (≈ 2.66). This certifies the theoretical background in Section III-F1. To have a closer look at how the PSD magnitude changes versus the frequency and time, the spectrograms of the radial AoA fluctuations and received power are respectively depicted in Figs. 7c and 7d. It is again confirmed in Fig. 7c that the AoA fluctuations have the main frequency content below 50 Hz and some small components up to 200 Hz. In Fig. 7d, the spectrogram of the received power also confirms the high PSD magnitudes below 50 Hz, indicating the dominant effect of the drone hovering-induced AoA fluctuations. Higher frequency components with  considerable magnitudes mostly appear from 100 Hz to 500 Hz, and sometimes reach to 2 kHz indicating stronger turbulence conditions. 7 It is seen from the spectrogram in Fig. 7d that the turbulence effect appears to be stronger during the first second, which is in agreement with the wide range of power fluctuations in the first second shown in Fig. 6a.

C. DRONE HOVERING AT 20-m ALTITUDE
In this section, following the same approach as in Section V-B, we examine another dataset when the drone was hovering at 20-m altitude, making a LoS distance of about 102 m and a total propagation distance of 204 m. Figure 8 again confirms that the fine-tracking system can compensate for the AoA fluctuations (shown in Fig. 8a) by 7 It is noted that the frequency components with high magnitudes can be interpreted from the yellow stripes that intermittently appear over time in Fig. 7d. controlling the FSM tracking angle (shown in Fig. 8b) to direct the beam centroid to the centers of the QD and PM. To verify this fact, the beam-centroid normalized position density illustrations at the QD and PM are shown, respectively in Figs. 8c and 8d. It is seen that the highest position density concentrates at the centers of the QD and PM and the standard deviations are also very small, i.e. (σ X = 17 µm, σ Y = 16 µ m) at the QD and (σ X = 8.83 µm, σ Y = 8.38 µm) at the PM. This suggests that the tracking accuracy at the PM is less than 9 µm, which is also suitable for coupling into a small detector or an MMF. Similar to Figs   theoretical PDF can be confirmed, where the means and variances are given as practical reference values for theoretical analyses and engineering designs.
In Fig. 9, the drone hovering-induced misalignments are comprehensively analyzed. Specifically, the AoA fluctuations at the receiving telescope entrance are shown in Fig. 9a, where the means the standard deviations in both axes are (µ X = 4.19 mrad, µ Y = 10.82 mrad) and (σ X = 1.17 mrad, σ Y = 2.67 mrad), respectively, which show a stronger displacement in the Y -axis. From the data in Fig. 9a, the normalized position density of the drone hovering trajectory with the means (µ p X = 0.42 m and µ p Y = 1.1 m) and standard deviations (σ p X = 0.12 m and σ p Y = 0.27 m) is depicted in Fig. 9b, which clearly displays the dominance of the random movements in the Y -axis. The beam-centroid radial displacement is also illustrated in Fig. 9c. The coherence time of the radial displacements is also found to be 700 ms in Fig. 9d, which again confirms the long coherence time of the misalignments. In Fig. 9e, the histogram data from Fig. 9c are then fitted to different PDFs in Section II-C, with the statistical metrics and fitting parameters summarized in Table 4. It is again verified that the Beckmann distribution is superior to other distributions in characterizing the beam centroid displacements, followed by the Hoyt, Rayleigh, and Rician distributions with small differences among them. To emphasize the small-scale mismatch at the tails of the PDFs, the log-log plot is provided in Fig. 9f, where all the fitted PDFs slightly overestimate the lower tail and underestimate the upper tail of the histogram data.
The received power data at the QD are shown in Fig. 10a and the normalized power data are then fitted to theoretical PDFs discussed in Section II-E, as depicted in Fig. 10b. Statistical metrics and fitting parameters are summarized in Table 5. It is again confirmed that the composite PDFs in (12) and (15) are better fit to the data than that in (17). A log-log plot is also examined in Fig. 10c, and it is observed that (12) and (15) marginally overestimate the lower tail of the data PDF with (12) showing the smallest discrepancy, while (15) performs better in the upper tail. In Table 5, the SIs of the turbulence and misalignments are both smaller than their counterparts for the case of 15-m hovering altitude in Table 3. This is logical as the turbulence is weaker at higher altitudes and as the total LoS distance is a little longer which consequently broadens the beam footprint at the receiver, the misalignments-induced fading also becomes weaker. In addition, the composite channel theoretical SIs match very well with the measured SI P Rx , with the discrepancy of 1∼1.3%.
Further investigations of the physical channel characteristics are done in Figs. 10d-10f. The channel coherence time of the composite channel can be determined as T c = 45 ms in Fig. 10d. This reveals a shorter channel coherence time,  which corresponds to a faster fading channel compared to Fig. 6d, and suggests that the turbulence has more contributions to the fading. The probability of fade is then plotted versus the fade margin F T in Fig. 10e together with the theoretical results. It is confirmed that the data are well predicted at F T ≤ 1 dB and gradually mismatched at higher values, which is in line with the lower-tail mismatch in Fig. 10c. To achieve a probability of fade at 10 −4 , a 2.5-dB power margin should be guaranteed, while the theoretical results suggest the margins of 3.8 dB and 4 dB, respectively. Finally, the LCR and AFD are investigated in Fig. 10f. As expected, the same trends as in Fig. 6f are witnessed, however, the theoretical lines fit very well with the experimental results when −1.5 ≤ ρ 0 ≤ 1.5 dB, which covers most of the dataset from Fig. 10a.
The PSD profiles of the radial AoA fluctuations and normalized received power are again investigated in Figs. 11a and 11b, respectively. The PSD of the radial AoA fluctuations in Fig. 11a shows the same trends as observed in Fig. 7a. In Fig. 11b, it is seen that the received power fluctuations have the main frequency content below f c = 55.7 Hz. It is also estimated that κ = 2.7449, which is very close to the theoretical value of 8/3 (≈ 2.66). Finally, Figs. 11c and 11d provide a closer look at the spectrograms of the radial AoA fluctuations and received power, respectively. It is confirmed in both figures that the highest PSD magnitudes are mostly found at frequencies below 50 Hz, which is a consequence from the domination of the drone hovering-induced misalignments. However, high PSD magnitudes also appear frequently at high frequency components up to 0.5∼2 kHz in the first 2.3 s in Fig. 11d, indicating more contributions from the turbulence. This is in agreement with the power data fluctuations in Fig. 10a, where stronger fluctuations are also seen in the first 2.3 s.

VI. CONCLUSION
In this paper, we experimentally investigated, for the first time, the physical channel characteristics of a drone-based retro-reflected FSO link, by conducting a beam propagation experiment between an OGS equipped with a fine-tracking system and a hexacopter drone equipped with a CCR over a maximum 204-m roundtrip distance. The experimental results were then compared with corresponding theoretical backgrounds, thereby offering valuable references for theoretical analyses as well as engineering designs of the considered system. It was highlighted that the LN and approximated Beckmann composite channel model provided the best fit to experimental data. In addition, the developed fine-tracking system showed a good tracking accuracy which is capable of coupling signals into a small detector or an MMF for high-speed signal processing techniques, while guaranteeing a compact size (30 × 30 cm) for practical deployments. This paper hence serves as an important guideline for the practical development of MRR-based drone-to-ground communication systems for fronthaul/backhaul links in the future 6G cellular networks. .

APPENDIX A A PROOF OF THE COMPOSITE CHANNEL PDF IN (17)
The PDF of the channel coefficient h can be expressed as where f h|h tc h|h t c denotes the conditional probability given a turbulence state h t c , expressed as Substituting (10) and (27) into (26), we have where = A mod Rh l . Then, plugging (6) into (28), we have Making a change of variable and with the help of [33, (3.383.4)], 8 the integration in (29) can be eventually solved as in (17). This completes the proof.

APPENDIX B NORMALIZED COMPOSITE CHANNEL PDFs
The PDFs in (12), (15), and (17) are normalized to the statistical mean, i.e. E [h] = Rh l E h t c E h p c , where E h t c = 1 and E h p c with h p c of a Beckmann RV can be written as [27] E h p c = A 0 ϕ x ϕ y 1 + ϕ 2