Since the directgradient (DG) method uses the ShackHartmann wave front sensor (SHWFS), based on the phaseconjugation principle, for atmospheric compensation in freespace optical (FSO) communication, it cannot effectively correct highorder aberrations. While the stochastic parallel gradient descent (SPGD) can compensate the distorted wave front, it requires more calculations, which is sometimes undesirable for an FSO system. A hybrid compensation (HC) method is proposed by properly using the DG method and SPGD algorithm to improve the performance of FSO communication. Simulations show that this method can well compensate wavefront aberrations and upgrade the coupling efficiency with few computations, preferable correction results, and rapid convergence rate.
I. INTRODUCTION
Freespace optical (FSO) communication is an advanced technology to implement lineofsight transmission of light signals. It transmits laser carrying signals through a freespace channel with high capacity, high bandwidth, and a flexible network
[1]
. FSO is one of the most promising alternative schemes for addressing the ‘last mile’ communication bottleneck in emerging broadbandaccess markets
[2]
. However, atmospheric disturbance can easily affect the laser beam. The refractive index of the atmosphere changes randomly because of variations in temperature, humidity, and wind speed in the atmosphere. This will lead to beam wandering, scattering, scintillation, and power fluctuations
[3]
. The phase and intensity of the laser are distorted, and the coupling efficiency decreases
[4
,
5]
. Consequently the bit error rate (BER) of the communication system is degraded. This arouses the interest of researchers worldwide in studying phase correction, and effective methods are in high demand for FSO communication systems.
Adaptive optical (AO) systems have been successfully applied to FSO to compensate the distorted wave front
[3

9]
. Based on the phaseconjugation principle, traditionally there are two main branches. One in
[6]
is based on a wavefront sensor, such as the SHWFS, to detect the local slope of the wave front and reconstruct the wave front in the Zernike or some other model. With a model, a deformable mirror (DM) can be controlled to construct a conjugated wave front, to offset the aberrations and obtain an approximately planar wave. In this method, the computational cost is very high for large numbers of subapertures in the sensor. Therefore the typically used sensors will fall short when high speed of detection and correction are needed, to implement fast adaptive optics in the FSO system. The other branch is the wavefront sensorless optimization method
[3
,
9]
, which aims to optimize the performance metrics of the received laser, such as Strehl ratio (SR), root mean square (RMS) and image sharpness functions, etc. The AO system searches for the suitable voltage to control the DM to optimize the performance metric. Many algorithms, such as simulated annealing (SA), hill climbing, and stochastic parallel gradient descent (SPGD) have been developed
[5]
. Among them, the SPGD algorithm is widely considered for its simple mechanism and rapid convergence, although hundreds of iterations are needed.
Recently, some new methods have been proposed. Reference
[7]
proposes a method to process ShackHartmann data by a focalplane approach. This method is more favorable in noise propagation, compared to classical ShackHartmann, and senses more phase modes with fewer subapertures under a comparable computation burden. A phaseretrieval method is used instead of the wavefrontslope method to reconstruct the wave front in reference
[8]
. The method provides more accurate estimation of aberrations in nearly flat wave fronts. Interestingly, reference
[10]
proposes a combined approach involving SPGD and DMmodelbased algorithms, to achieve similar correction results to those of SPGD with many fewer iterations. But in this method a SHWFS is not effective, and both SPGD and DMmodelbased algorithms need parameter settings that will increase the complexity of the method. In addition, a trustregion method has been proposed that is superior to both SA and SPGD algorithms, with respect to convergence rate for slowly changing wavefront aberrations, in reference
[11]
.
Considering the rapidly changing atmospheric environment in an FSO system, we propose an HC scheme for wavefront corrections. This method combines the DG method
[12]
and SPGD algorithm to improve both the results of wavefront compensation and convergence rate.
This paper is organized as follows: In section II, the principles of the DG method and SPGD algorithm are briefly introduced and their deficiencies discussed, then we propose an HC scheme based on a suitable combination of the methods discussed above. In section III, computer simulations using MATLAB are carried out and the results are analyzed, to investigate the compensation capability of the HC scheme. Finally, we give our conclusion in section IV.
II. PROPOSED HYBRID COMPENSATION METHOD
The fundamental scheme of FSO communication is shown in
Fig. 1
. We adapt optical intensity modulation, and the data signals directly modulate the light source to generate optical signals. When the optical signals are transmitted to an atmospheric channel, they may suffer from atmospheric turbulence, and the receiver may receive optical signals with distorted wave fronts. Therefore, at the receiver we use an AO system to compensate for the aberrations of the distorted wave front. Then the optical signals with corrected wave fronts are coupled into the fiber and detected by the photodetector, to recover the original data signals.
Fundamental scheme for FSO communication.
As we know, a wave front distorted by atmospheric turbulence is commonly described by the Zernike polynomials
[13]
, which are a set of polynomials defined on a unit circle. It is convenient to use polar coordinates, so that the polynomials are a product of angular functions and radial polynomials. The wave front is described as
[14]
where (
r
,
θ
) are the polar coordinates of the pupil,
Z_{i}
(
r
,
θ
) is the
i
th order Zernike polynomial,
a
i is the
i
th coefficient, and
q
is the highest order Zernike polynomial considered. The expression
ϕ
(
r
,
θ
) can be transformed into rectangular coordinates
ϕ
(
x
,
y
) for calculation.
To compensate a distorted wave front, the DG method is frequently used to detect and calculate a wavefront gradient matrix
G
using SHWFS and the optimal control
V
of DM, which is the leastsquares solution of
where
G
= [
G
_{1x}
,
G
_{2x}
,
⋯
,
G_{Mx}
,
G
_{1y}
,
G
_{2y}
,
⋯
G_{My}
,]
^{T}
contains the local slopes in the horizontal direction
x
and vertical direction
y
, detected in every subaperture, and the subscript
M
is the number of subapertures of the SHWFS. The matrix K = (
G
_{1}
,
G
_{2}
,
⋯
,
G_{N}
) is defined as the gradient response matrix of the DM, and
N
is the number of control elements of the DM.
For later use, we give a brief interpretation of Eq. (2). We suppose that the voltages applied to the actuators of the DM are linear, and let
Ψ
(
x
,
y
) be the wave front generated by the
N
actuators, i.e.
where (
x
,
y
) are the coordinates of the wave front generated by the DM and
I_{n}
(
x
,
y
), n = 1, 2,
⋯
,
N
is the influence function
[5]
of the
n
th actuator, and
where (
x_{n}
,
y_{n}
) are the central coordinates of the
n
th actuator,
ω
is the coupling coefficient,
b
is the normalized interval between the adjacent actuators, and
α
is the Gaussian index.
On the other hand, the mean local gradient of the distorted wave front in the
m
th (
m
=1, 2,
⋯
,
M
) subaperture is regarded as the mean of the partial derivative of the wave front
ϕ
(
x
,
y
) in the area
S_{m}
, which is the local area of the
m
th subaperture
[6]
. Based on Eq. (3) and (4), we obtain
where
Since this method directly calculates the optimal control voltage by the gradient of the wave front, without wavefront reconstruction, it can reduce the computational cost of the AO system. From Eq. (2) we can see that the influence function of DM is a key factor in compensation. To investigate the compensation capacity of the DM, we introduce the definition of compensation error as
where RMS
_{o}
is the RMS value of the original distorted wave front, and RMS
_{c}
is the RMS value of the corrected wave front. To test the general correction ability of the DM, the first 35 Zernike polynomials are taken as the targets for the correction, and all of their amplitudes are normalized. Here we consider a DM with 61 elements for the analysis. The compensation error of the 61element DM for each Zernike order from 3 to 35 is shown in
Fig. 2
. From this figure, we can see that the compensation error increases with increasing Zernike order, which means the DG method is relatively incapable of correcting highorder aberrations.
Compensation error of a 61element DM for each Zernike order.
For this reason, we introduce the SPGD algorithm for the sensorless wavefront compensation given in
[3]
to correct highorder aberrations. In this algorithm the object of optimization is a performance metric
J
for the received laser. This algorithm searches for optimal control
u
= (
u
_{1}
,
u
_{2}
,
⋯
,
u_{N}
) and obtains the optimal
J
iteratively. In addition, the SPGD controller generates a set of statistically independent control perturbations {
δu_{j}
^{(k)}
},
j
= 1, 2,
⋯
,
N
for each iteration
k
to update the control
u
. The search procedure is carried out as follows:
where
and
γ
is the gain coefficient. If the objective of optimization is to achieve the maximum value of
J
, then
γ
is positive, else
γ
is negative.
For SPGD, given
γ
and
J
, the gradients of the wave front follow the direction of descent with iterations. When the atmospheric disturbance increases, leading to a more distorted wave front, the correction by the SPGD algorithm will act to reduce this distortion.
Considering the effect of increasing convergence rate on the correction and the performance metric, we make use of the mature technologies of the DG and SPGD algorithms to propose an HC scheme based on a combination of the properties of these methods. For a distorted wave front, the DG method can compensate loworder aberrations without iteration, while the residual wave front consisting of a large proportion of highorder aberrations will be compensated by SPGD. In this way, the initial wave front corrected by SPGD is a less distorted wave front, with a smoother phase plane than the wave front of the inserted laser, which means much fewer iterations are needed to reach convergence. Moreover, our proposed HC scheme uses the SHWFS to compensate the loworder aberrations, so that the accuracy requirement for the sensor is reduced, and we do not have to choose a SHWFS with many subapertures. Large number of subapertures will divide the inserted laser power into many small parts, which causes difficulties in CCD detection. Furthermore, the computation pressure is relieved, since the size of
G
is 2
M
×1, which is proportional to the number of subapertures.
In this paper, we adapt a hybrid scheme by combining the mature schemes of the DG method and SPGD algorithm. The functional blocks are diagrammed in
Fig. 3
, in which the wavefront corrector is the DM shared by the DirectGradient Block and SPGD Block. The main process of this scheme is given as follows.
Functional block diagram of the HC method.

(1) When the light arrives, the Algorithm Selection Block identifies whether the inserted wave is the original, distorted wave or the residual wave, then selects the corresponding correction block.

(2) If it is the original, distorted wave, the DirectGradient Block begins to detect the wavefrontgradient matrixGand calculatesVof DM using Eq. (2).

(3) The wavefront corrector generates a conjugated wave front to compensate the original wave front, and a residual wave is output.

(4) The residual wave propagating back to the Algorithm Selection Block, is switched to the SPGD Block, and then iteratively corrected based on Eq. (8). With the SPGD algorithm, the wave front is iteratively corrected until the value of the performance metricJconverges.

(5) After correction, the optimized wave is coupled into the singlemode fiber and processed by the receiver. A diagram of this entire process is presented inFig. 4.
Diagram of the HC process.
For FSO, the coupling efficiency of the signals inserted into the fiber in the receiver has significant influence on system performance
[15]
. Thus in the SPGD block the performance metric
J
is set equal to the SR, which is the ratio of the distorted laser’s peak intensity to that of the ideal laser. As shown in Eq. (11)
[16]
, SR is an estimation of the coupling efficiency and can directly reflect the aberration effect
where
A_{r}
(
f
) is the Fourier transform of the singlemode fiber's optical field and
f
is the frequency of the optical wave. The coupling efficiency increases as SR improves.
On the other hand, the BER of FSO communication links with onoff keying systems increases as RMS rises
[17]
. The relationship between RMS and SR is expressed as
[16]
This means that the BER decreases with increasing SR. In other words, improvement of SR both increases coupling efficiency and reduces BER.
In the following section, by observing the SR of the system with our proposed compensation scheme, we present simulation results and analysis of wavefront compensation, to investigate the performance of our HC algorithm.
III. SIMULATION RESULTS AND ANALYSIS
In our simulations, for the DirectGradient Block we use the SHWFS with an 8×8 subaperture array, in which the subapertures are not considered in every corner of the array shown in
Fig. 5
, to avoid detection difficulties in low light
[6]
. For the SPGD Block, the perturbations {
δu_{j}
} obey a Bernoulli distribution with probability 0.5 each for
δu
= +
σ
and
δu
= 
σ
. We define
γ
=
γ
_{0}
/(
J
+
C
), where C is a constant to avoid too large
γ
when
J
is sufficiently small. Parameter settings for the SPGD algorithm and HC method are the same in our simulations. The wavefront corrector is a DM containing 61 control elements, to realize a tradeoff between an accurate influence function and a suitable aperture size in the FSO system. The arrangement of the 61 actuators is shown in
Fig. 6
.
8×8 subaperture array.
Arrangement of the 61element DM.
The distorted wave front that we introduce is the superposition of Zernike polynomials with coefficients
a_{i}
listed in
Table 1
, where
a_{i}
is set by the covariance matrix of Zernike terms
[14]
. We define the coefficient vector
A
= [
a
_{1}
,
a
_{2}
,
⋯
,
a
_{20}
]
^{T}
and use the covariance matrix (Eq. (4) in
[14]
) to evaluate
A
according to Eqs. (8), (12), and (13) presented in
[14]
. The SR of the distorted wavefront is 0.1668.
Coefficients of Zernike polynomials
Coefficients of Zernike polynomials
Using this distortedwavefront sample, we first simulate the DG method; after correction, the corresponding SR increases to 0.6958. The SRs of the iteratively corrected residual wave fronts are given in
Fig. 7
. From this figure, we can see that this method cannot further improve SR with increasing iterations, for highorder aberrations. Thus in our HC scheme the correction is carried out only once, to compensate the loworder aberrations with the DG algorithm and the highorder aberrations with the SPGD method.
Improvement of SR with DG, as a function of number of iterations.
Next we conduct the simulation of our HC scheme and compare its performance to that of the SPGD algorithm in correcting the original, distorted wave front. The improvements in SR for the HC and SPGD algorithms are presented in
Fig. 8
. We find that the HC scheme can converge much faster than the SPGD algorithm. In
Fig. 8
, when the SR is 0.8 the SPGD algorithm requires 144 iterations, but the HC scheme only needs 21. When the SR becomes 0.9, the number of iterations is about 370 by SPGD, while 110 by the HC scheme. Since the loworder aberrations have been corrected, the SR increases to 0.6958, so that further corrections implemented by the SPGD method need much fewer iterations.
Improvement of SR with HC and SPGD algorithms.
Figure 9
shows the differences of the phase plane before and after correction.
Figure 9(a)
is the original, distorted wave front (SR=0.1668). It is obvious that it becomes flatter after correction by DG (SR=0.6958), as shown in
Fig. 9(b)
, and it gets smooth after correction by the HC scheme (SR=0.9), as shown in
Fig. 9(c)
. The range of the wave front narrows correspondingly.
Comparison of different wave fronts: (a) original wave front (SR=0.1668), (b) corrected by DG (SR=0.6958), and (c) corrected by HC (SR=0.9).
Furthermore, we discuss the intensity characteristics of the laser beam to investigate the improvement of the coupling efficiency. The intensity
I
(
x
,
y
) is normalized to the peak intensity of the ideal plane wave, where the coordinates (
x
,
y
) identify a pixel in the focus plane of the detector.
Figure 10
shows the normalized intensity of the original, distorted light, the DGcorrected light, and the HCcorrected light. The peak intensity rises from 0.1098 to 0.4319 after correction by the DG method. With correction by the HC scheme, when SR is 0.9 in
Fig. 9(c)
the peak intensity is 0.5549, as shown in
Fig. 10(c)
. This illustrates that the HC scheme can enhance the light's intensity by almost a factor of 5. With the increase in SR, more signal power will be coupled into the singlemode fiber so that the coupling efficiency increases, which helps the FSO system to correctly detect the signals.
Different light intensities: (a) original, distorted light (SR=0.1668), (b) corrected by DG (SR=0.6958), and (c) corrected by HC (SR=0.9).
To verify the compensation effect of the HC scheme, besides the wave front discussed above, we introduce another seven wave fronts, simulated by Zernike polynomials with different maximum Zernike orders
q
(
q
= 20, 25, 30
and
35) and defined as
ϕ
_{1}
,
ϕ
_{3}

ϕ
_{8}
respectively. The wave front already discussed is
ϕ
_{2}
. The initial SRs of the eight distortedwavefront samples
ϕ
_{1}
~
ϕ
_{8}
are 0.1604, 0.1668, 0.0493, 0.3072, 0.1064, 0.1667, 0.0120, and 0.3760 respectively. We compare the compensation effect of the DG method, SPGD algorithm, and HC scheme for every sample, to get a clear picture of the HC scheme. The improvements in SR with the SPGD and HC schemes of the eight samples
ϕ
_{1}

ϕ
_{8}
are shown in
Fig. 11
.
Improvement of SR for ϕ_{1}ϕ_{8} (a)~(h) with the HC or SPGD method.
From
Fig. 11
we can see that the HC scheme is obviously superior to the SPGD algorithm, because it converges more rapidly than SPGD for every sample, especially for samples
ϕ
_{3}
and
ϕ
_{7}
with quite low initial SR. As shown in
Figs. 11(c)
and
(g)
, by correction with SPDG, when SR is lower than 0.1 the rate of SR improvement with iterations is quite low, but the curve becomes much steeper beyond SR=0.1. This means that the correction speed of SPGD is relative to the value of SR. Since loworder aberrations make up a large fraction of the wavefront aberrations, it is difficult for SPGD to correct the loworder aberrations. The HC scheme can compensate the loworder aberrations by the DGBlock; thus it increases the correction speed of SPGD. Additionally, our HC scheme provides higher SR than the SPGD algorithm within 1000 iterations.
Next we will perform a numerical analysis of the results. We define
J_{O}
as the SR of the original, atmospheric distorted wave front, and
J_{DG}
as the SR of the residual wave front after correction by the DG method.
J_{H}
and
J_{S}
are respectively the correction results of the HC and SPGD algorithms after 1000 iterations, as shown in
Fig. 10
. Furthermore,
k_{H}
^{0.8}
and
k_{S}
^{0.8}
are the minimum iteration numbers that the HC scheme and SPGD algorithm need respectively for SR to reach 0.8. The simulation results for the eight samples are presented in
Table 2
, which significantly shows that for every sample the HC scheme can reach much higher SR than the DG method, especially for samples
ϕ
_{5}

ϕ
_{8}
with higher Zernike orders. This verifies that the HC scheme can utilize the SPGD algorithm to offset the deficiency of the DG method in highorder aberration correction.
Correction results for different distorted wave fronts
Correction results for different distorted wave fronts
To obtain convenient comparisons of correction results for different samples, we define the indice of the correction effect of the HC scheme as
η
and the indice of its correction speed Δ
k
^{0.8}
, precisely
For the samples enjoying better corrections by our HC scheme,
η
is larger, and for the samples with faster correction by the HC scheme than by SPGD, Δ
k
^{0.8}
is larger. The indice values for the eight samples are shown in
Table 3
.
Indices of correction effect and speed for the HC method
Indices of correction effect and speed for the HC method
As shown in
Table 3
, for samples of the same Zernike order, those with lower initial SR can realize higher
η
and Δ
k
^{0.8}
. For example, for samples
ϕ
_{5}
and
ϕ
_{6}
simulated by 30 Zernike polynomials, the initial SRs are 0.1064 and 0.1667 respectively,
η
are 7.43 and 4.52, and Δ
k
^{0.8}
are 266 and 172. The results keep consistent for other samples. This means that under the same wavefront simulation conditions, the HC scheme can perform better at compensation of wave fronts with more severe aberrations. This benefits from the DGBlock’s compensation, because the DG method plays an important role in compensating the loworder aberrations. By DGBlock correction, the initial wave front for the SPGD algorithm to correct is smoother and the SR increases, since the correction speed of SPGD depends on the value of the SR.
In summary, from simulations of different distorted wave fronts and the analysis of the results, we know that through a suitable combination of the DG method and the SPGD algorithm our HC scheme can achieve higher performance metrics than that of the DG method and convergent faster than the SPGD algorithm when loworder aberrations are compensated in advance by the DG method. After correction by the HC scheme, SR increases, so that the light intensity increases and more power is transmitted to the receiver, for higher coupling efficiency and improved communication quality.
IV. CONCLUSION
In this paper we have proposed a hybridmethod HC scheme to compensate wavefront distortions resulting from the turbulent atmosphere in FSO communication. This scheme reduces the computational complexity of the compensation process and accelerates the convergence of the performance metric, to adapt to the realtime requirements of FSO communication. Based on the distorted wave fronts generated by Zernike polynomials, we compare the compensation performance of the SPGD algorithm, DG method, and our HC scheme by computer simulations. We find that with the SPGDBlock’s correction our HC scheme can achieve higher SR than the DG method, especially for wave fronts with higher Zernike orders. And the HC scheme can increase the correction speed of SPGD by increasing SR in the DGBlock, especially for severely distorted wave fronts, because high SR requires fewer iterations than low SR does. Our method can help an FSO system to improve coupling efficiency of the laser to the receiver and decrease BER.
Acknowledgements
This study is supported by the National Natural Science Foundation of China (No.61171079).
Ghassemlooy Z.
,
Popoola W.
,
Rajbhandari S.
2012
Optical Wireless Communications System and Channel Modelling with MATLAB
CRC Press
Boca Raton, FL, USA
Chaudhary S.
,
Amphawan A.
(2014)
“The role and challenges of freespace optical systems,”
Journal of Optical Communications
4
327 
334
Weyrauch T.
,
Vorontsov M. A.
(2005)
“Atmospheric compensation with a speckle beacon in strong scintillation conditions: directed energy and laser communication applications,”
Appl. Opt.
44
6388 
6401
Wei L.
,
Shi W. X.
(2013)
“Free space optical communication performance analysis with focal plane based wavefront measurement,”
Opt. Commun.
309
212 
220
Li Z. K.
,
Cao J. T.
,
Zhao X. H.
(2015)
“Atmospheric compensation in free space optical communication with simulated annealing algorithm,”
Opt. Commun.
338
11 
21
Lane R. G.
,
Tallon M.
(1992)
“Wavefront reconstruction using a ShackHartmann sensor,”
Appl. Opt.
32
6902 
6908
Meimon S.
,
Fusco T.
,
Michau V.
,
Plantet C.
(2014)
“Sensing more modes with fewer subapertures: the LIFTed ShackHartmann wavefront sensor,”
Opt. Lett.
10
2835 
2837
Li J.
,
Gong Y.
,
Chena H. F.
,
Hu X. R.
(2015)
“Wavefront reconstruction with HartmannShack sensor using a phaseretrieval method,”
Opt. Commun.
336
127 
133
Vorontsov M. A.
,
Sivokon V. P.
(1998)
“Stochastic parallelgradientdescent technique for highresolution wavefront phasedistortion correction,”
J. Opt. Soc. Am. A
10
2745 
2758
Dong B.
,
Yu J.
(2015)
“Hybrid approach used for extended imagebased wavefront sensorless adaptive optics,”
Chin. Opt. Lett.
13
041101 
Yang Q. Y.
,
Zhao J. Y.
,
Wang M. H.
,
Jia J. L.
(2015)
“Wavefront sensorless adaptive optics based on the trust region method,”
Opt. Lett.
7
1235 
1237
Jiang W. H.
,
Li H. G.
(1990)
“HartmannShack wavefront sensing and wavefront control algorithm,”
Proc. SPIE
1271
82 
93
Noll R. J.
(1976)
“Zernike polynomials and atmospheric turbulence,”
J. Opt. Soc. Am.
3
207 
211
Roddier N.
(1990)
“Atmospheric wavefront simulation using Zernike polynomials,”
Opt. Eng.
10
1174 
1180
Dikmelik Y.
,
Davidson F. M.
(2005)
“Fibercoupling efficiency for freespace optical communication through atmospheric turbulence,”
Appl. Opt.
29
4946 
4952
Li Z. K.
,
Cao J. T.
,
Zhao X. H.
,
Liu W.
(2014)
“Combinationaldeformablemirror adaptive optics system for atmospheric compensation in free space communication,”
Opt. Commun.
320
162 
168
Yang Y. Q.
,
Han Q. Q.
,
Tan L. L.
(2011)
“Influence of wavefront aberrations on bit error rate in intersatellite laser communications,”
Opt. Commun.
284
3065 
3069