各向异性大气湍流中成像仿真研究

2022.06.15

崔林艳

摘要：大气湍流降低成像质量，一直以来是制约远距离光电成像探测的瓶颈问题。通过对大气湍流中成像仿真进行研究，一方面可以评估光电成像探测系统性能，同时还可以为图像处理算法验证、目标识别等提供测试数据。本文主要针对各向异性非Kolmogorov大气湍流中成像仿真开展研究，通过利用图像处理和各向异性非Kolmogorov大气湍流理论成像模型相结合的方法，实现了综合考虑各向异性系数、成像波长、接收器孔径、传播路径等物理参量的各向异性湍流中成像快速仿真。通过与随机相位屏仿真方法和图像处理成像仿真方法分别从仿真时间和仿真精度两个方面进行对比分析，结果表明本文方法仿真时间短，同时也能较好地反映大气湍流对光学成像系统的影响。

关键词：各向异性；大气湍流；湍流谱指数；图像畸变；图像模糊

中图分类号： P427.1； TN29 文献标识码： A文章编号： 1673-5048（2018）03-0058-06

0 引言

大气湍流中探测成像时，湍流退化效应会造成图像的模糊、畸变，影响成像质量和成像系统性能。通过对大气湍流中成像进行仿真，模拟大气湍流的不同退化效应，可以预测和评估光电成像探测系统性能，同时还可以为图像质量评估、图像处理算法验证、目标识别等提供测试数据。

目前大气湍流中成像仿真方法主要是针对经典的Kolmogorov大气湍流展开的，大致可分为基于随机相位屏的成像仿真方法[1-2]、图像处理成像仿真方法[3]和图像处理与湍流物理模型相结合的成像仿真方法[4-6]三类。

基于随机相位屏的成像仿真方法，从光传播的角度出发，利用随机相位屏计算光波相位畸变量，结合角谱传播理论，实现大气湍流中成像仿真。该方法从物理角度考虑了影响成像质量的各种大气湍流参量，仿真时需要根据大气湍流参量及成像条件来判断是等晕成像还是非等晕成像。在等晕成像条件下，利用随机相位屏仿真方法计算光波畸变量，进行大气湍流降质图像仿真[1-2]。在非等晕成像条件下，需要根据等晕角确定满足局部等晕的区域块大小，然后对每个区域块分别进行随机相位屏成像仿真，最终得到大气湍流降质图像。在采用随机相位屏成像仿真方法时，需要在成像路径上设置多层相位屏，并且在每层计算过程中均会涉及到傅里叶变换、反傅里叶变换等数学操作。因此，分块数量较多时，非等晕成像仿真的时间较长。

图像处理成像仿真方法通过人工设定图像畸变量和图像模糊点扩散函数，结合图像处理算法（图像反映射和双线性插值操作、图像卷积操作等），仿真得到大气湍流等晕和非等晕降质图像。该方法仿真时间较短，但是仿真过程中没有从物理角度考虑大气湍流参量的影响，不能从物理角度定量反映不同大气湍流条件下的成像情况。

成像系统重点关注大气湍流对光波相干性破坏造成的光传播相位相关的物理量，包括到达角起伏造成的图像畸变、光束扩展造成的图像模糊等。图像处理与湍流物理模型相结合的成像仿真方法利用Kolmogorov湍流理论计算得到大气湍流效应（到达角起伏、调制传递函数等），从而得到图像畸变量和图像模糊量，然后利用图像处理算法将这些湍流效应添加到图像中，仿真得到大气湍流等晕和非等晕成像条件下的降质图像。该方法从物理角度考虑了大气湍流参量的影响，并且在仿真时间方面比随机相位屏成像仿真方法具有明显优势。

长期以来，大气湍流的研究主要是在Kolmogorov大气湍流的统计均匀、各向同性假设的基础上展开的。近年来随着探测设备的不断发展和理论研究的不断深入，研究结果表明Kolmogorov大气湍流并不是湍流中的唯一存在形式，例如在高于地面2 km以上的大气层很多时候呈现各向异性和非Kolmogorov的特性。各向异性非Kolmogorov大气湍流的理论研究工作引起国内外学者的广泛关注[7-16]。而对于各向异性非Kolmogorov（nonKolmogorov， N-K）大气湍流中成像仿真，尤其是图像处理与各向异性N-K湍流物理模型相结合的成像仿真方法，国内外还没有相关研究工作。

1 各向异性N-K大气湍流中成像退化物理模型

光电成像系统重点关注与光波相位相关的物理量。因此，本文主要研究大气湍流造成的图像模糊和图像畸变这两种与光波相位密切相关的图像降质效应。各向异性N-K大气湍流造成的图像模糊效应在时域内可以用点扩散函数PSF来描述，在頻域内可以用各向异性N-K大气湍流MTF来描述。各向异性N-K这一湍流效应引起的畸变量可由各向异性N-K大气湍流到达角起伏方差来估计。通过将理论计算得到的各向异性N-K大气湍流效应模型与图像处理算法结合起来，仿真得到各向异性N-K大气湍流降质图像，可表示为

其中： O（u，v）表示仿真输出的各向异性N-K大气湍流降质图像； I（x，y）表示输入的清晰原始图像； F（·）和R（·）分别为仿真图像模糊效应和图像畸变效应的操作算子。其均可由前期推导得到的各向异性N-K大气湍流MTF理论模型和到达角起伏方差理论模型计算得到。

1.1 各向异性N-K大气湍流中图像畸变仿真原理

当光在大气湍流介质中传播成像时，到达角起伏湍流效应造成图像各像素点位置发生随机偏移。在等晕成像情况下，图像各像素点的偏移量是相同的，即图像发生整体平移运动。在等晕成像仿真时，可以不用考虑图像畸变。在非等晕成像情况下，图像各像素点在水平方向和垂直方向上的偏移量不同。各向异性N-K大气湍流中图像畸变仿真流程如图1所示。

具体实现步骤如下：

（1）产生频域内的复随机数矩阵Rw=A+iB，其中A和B为满足均值为0、方差为1的高斯函数分布的伪随机数矩阵。引入该随机数矩阵的目的是要仿真大气湍流的随机起伏特性。

（2）计算与实际成像设备相符合的能反映图像畸变的空间功率谱密度函数Sβ（f）： Sβ（f）=f-b，b=3。该函数描述了大气湍流引起的图像不同区域几何畸变量之间的空间相关情况。

（3）用步骤（1）所产生的复高斯随机数矩阵对Sβ（f）进行滤波，得到一个既能反映图像畸变量空间相关特性又能反映大气湍流随机变化特性的频域内复随机场，可以表示为[4-6]： HβRw=Sβi2+j2ΔκΔκ·Rw。与Rw相比， HβRw考虑了图像不同区域之间的相关性。与现有研究方法[4-6]不同的是，对于零频部分，采用了与随机相位屏谱反演法相似的次谐波低频补偿方法对其进行补偿。

（4）对HβRw进行反傅里叶变换并取模值，得到空域内矩阵IFFTHβRw。将IFFTHβRw的方差设定为N-K大气湍流到达角起伏方差σ2AOA，得到了N-K大气湍流引起的图像水平方向的畸变量矩阵（单位为rad）。其中， σ2AOA为各向异性N-K大气湍流中平面波到达角起伏方差理论模型[14]，其考虑了各向异性系数、湍流谱指数、传播成像距离、湍流折射率结构常数等参量的影响。

（5）重复步骤（1）～（4），得到了各向异性N-K大气湍流引起的图像垂直方向的畸变量矩阵（单位为rad）。

（6）将计算得到的图像水平和垂直方向畸变量矩阵除以角分辨率Δθ （Δθ=ΔδL，Δδ为图像每个像素代表的实际尺寸； L为传播距离），就得到各向异性N-K大气湍流引起的图像水平方向和垂直方向偏移矩阵δx和δy（单位为像素）。

（7）根据矩阵δx和δy，采用反映射和双线性插值图像处理算法，仿真得到了各向异性N-K大气湍流畸变图像。

1.2 各向异性N-K大气湍流中图像模糊仿真原理

其中： *表示卷积操作； PSF为大气湍流点扩散函数，可由前期理论推导得到的各向异性N-K大气湍流MTF模型[13]的反傅里叶变换并取模值得到。由理论研究可知，各向异性N-K大气湍流MTF所对应的PSF近似为高斯函数，而该函数可由标准差值来表征。标准差取值越大， PSF越宽，即大气湍流对图像造成的模糊效应越明显。相反，则图像模糊效应越不明显。利用式（2）进行卷积操作，得到大气湍流等晕模糊降质图像。仿真流程如图2所示。

在非等晕成像情况下，大气湍流对图像造成的模糊为空变模糊。大气湍流随机性和空间相关性造成图像上不同区域对应的PSF呈现出随机性并且满足一定空间相关性。由于大气湍流MTF所对应的PSF近似为高斯函数，因此获取PSF的关键是获取图像不同区域所对应的PSF的标准差大小。各向异性N-K大气湍流中非等晕图像模糊仿真流程如图3所示。

具体实现步骤如下：

（1）首先产生频域内的复高斯伪随机数矩阵Rw=A+iB，其中A和B为满足均值为0、方差为1的高斯函数的伪随机数矩阵。引入该随机数矩阵的目的是要仿真大气湍流的随机起伏特性。

（2）计算与实际成像设备相符合的能反映图像模糊的空间功率谱密度函数： Sα（f）=f-a， a的经验值为3。该函数表示了大气湍流引起的图像不同位置模糊效应之间的空间相关情况。 Sα（f）与Sβ（f）具有相同的表达式，但物理含义不同，分别表征了大气湍流对图像模糊和图像畸变的空间相关情况。

（3）用步骤（1）所产生的复高斯随机数矩阵对Sα（f）进行滤波，得到一个既能反映图像模糊量空间相关特性又能反映大气湍流随机变化特性的频域内复随机场，可以表示为[4-6]： HαRw=Sαi2+j2ΔκΔκ·Rw。与Rw相比， HαRw考虑了图像不同区域之间的相关性。对于零频部分，与现有研究方法[4-6]不同的是，采用了与随机相位屏谱反演法相似的次谐波低频补偿方法对零频部分进行补偿。

（4）对HαRw进行反傅里叶变换并取模值，得到空域内矩阵IFFTHαRw。将IFFTHαRw的均值设定为各向异性N-K大气湍流MTF模型所对应的点扩散函数（近似为高斯函数）的标准差σt，均方差设定为σt，得到了图像不同区域所对应的PSF标准差大小的矩阵PSF0（单位为m）。在前期工作中已经推导建立了各向异性N-K大气湍流MTF理论模型[13]，考虑了各向异性系数、湍流谱指数、传播成像距离、湍流折射率结构常数等参量的影响。

（5）将PSF0除以Δδ（图像每个像素代表的实际尺寸），就得到了图像不同区域所对应PSF（为高斯函数）标准差大小的矩阵PSF（单位为像素）。其中， PSF（i， j）表示圖像（i， j）处所对应的PSF标准差大小。

2 仿真结果与分析

对于激光成像系统来说，通常的波长为1.06 μm或1.55 μm，因此在本文仿真中，采用的参量为：波长1.55 μm， α=10/3， C^2n=1×10-14 m-1/3，各向异性因子ζ=2，成像传播距离为4 km，图像大小为512×512像素，每个像素点所代表的尺寸为10 mm，接收器孔径大小为40 mm。

2.1 各向异性N-K湍流中图像畸变量和模糊量计算

在仿真过程中，如果没有考虑空间相关性，则图像中各像素点的几何畸变量是随机起伏的，如图4（a）～（b）所示。

由于引入了能反映图像畸变空间相关性的空间功率谱密度函数Sβ（f），图像相邻像素点之间不会出现突变的几何畸变量，并且整幅图像几何畸变量表现出较为连续的变化趋势，如图4（c）～（d）所示。

按照1.2节大气湍流中图像模糊仿真流程，计算得到图像上各像素点对应的模糊PSF（i，j）的均方差值。如果没有考虑空间相关性，则图像中各点对应的PSF（i，j）是随机起伏的，如图5（a）所示。由于引入了能反映图像空间相关性的空间能量谱密度函数Sα（f），图像各像素点对应的PSF（i，j）存在较好的相关性，如图5（b）所示。

2.2 各向异性N-K湍流中成像仿真结果及分析

利用2.1节计算得到的各向异性N-K大气湍流中图像畸变量和图像模糊量，并结合图像双线性插值算法和图像模糊卷积算法，首先对图像进行几何畸变操作，然后再进行模糊操作，进而得到各向异性N-K大气湍流降质图像，仿真结果见图6。在图像仿真时，采用MATLAB R2011a版本，仿真时间为350 s。

从仿真时间和仿真精度两个方面与随机相位屏仿真方法进行对比分析。在仿真时间上，待仿真图像满足非等晕成像条件，此时需要将每个像素点（看做一个点光源）根据随机相位屏方法进行成像仿真。图像大小为512×512，每个像素点成像仿真所需要的时间为9 s，则仿真一幅图像所需要的时间为512×512×9 s，约为655 h。即使将随机相位屏的数量由13个简化为1个，即用1个相位屏来代替大气湍流的影响（此时仿真精度下降），每个像素点的仿真计算时间为0.7 s，仿真一幅图像所需要的时间为512×512×0.7 s，约为51 h。当需要生成大气湍流降质视频序列时，则需要的计算量是非常庞大的。巨大的计算量不利于获取不同湍流强度和不同传播距离下的大气湍流视频序列，不利于大气湍流中成像问题的分析。相比之下，本文仿真方法在仿真时间上具有明显优势。在仿真精度上，随机相位屏仿真结果不仅考虑了与光波相位相关联的图像模糊和图像畸变，同时还考虑了与光波振幅相关联的光强起伏等其他湍流效应。理论上，仿真精度比本文算法高。由于光强起伏并不是成像系统所要考虑的重点，本文方法的仿真精度基本满足成像系统工程的需要。

从仿真时间和仿真精度两个方面与图像处理成像仿真方法（需要人工设定图像畸变量和图像模糊量）进行比较分析。在仿真时间上，与本文仿真方法相比，图像处理成像仿真方法没有计算大气湍流效应数值这一步骤（计算时间为4.9 s），因此在仿真时间上后者稍快。在仿真精度上，图像处理成像仿真方法通过人工设定不同的图像畸变量和图像模糊量，定性仿真不同湍流强度和不同传播距离下的大气湍流降质图像。而本文仿真方法可以根据具体的成像场景，并结合各向异性N-K大气湍流理论计算出各向异性N-K大气湍流造成的图像畸变量和图像模糊量，从而能从物理角度定量仿真不同湍流强度和不同传播距离下的各向異性N-K大气湍流降质图像。

综上所述，本文采用的各向异性N-K大气湍流中成像仿真方法比随机相位屏仿真方法和图像处理成像仿真方法具有更好的工程应用价值。

3 结论

本文采用各向异性N-K大气湍流理论与图像处理算法相结合的仿真方法对各向异性N-K大气湍流中成像进行了仿真。仿真图像能较好反映大气湍流对成像造成的模糊和畸变现象。通过与随机相位屏仿真方法和图像处理成像仿真方法从仿真时间和仿真精度两个方面进行对比分析，表明本文仿真方法具有较好的工程应用价值（仿真时间短，同时也能较好反映大气湍流对光学成像系统的影响）。由于目前各向异性大气湍流的实验研究工作还不完善，即无法得到真实场景的各种湍流物理参量，包括各向异性系数、湍流谱指数等，本文仿真结果无法与真实结果进对比分析。

参考文献：

[1] Zhang Jianlin， Zhang Qiheng， He Guangming. Multiframe Blind Image Restoration Based on a Multiplicative Iterative Algorithm [J]. Optical Engineering， 2009， 48（2）： 027004-1-4.

[2] Zhang Jianlin， Zhang Qiheng， He Guangming. Blind Image Deconvolution by Means of Asymmetric Multiplicative Iterative Algorithm [J]. Journal of the Optical Society of America A， 2008， 25（3）： 710-717.

[3] Zhu X， Milanfar P. Removing Atmospheric Turbulence via SpaceInvariant Deconvolution [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence， 2013， 35（1）： 157-170.

[4] Repasi E， Weiss R. Analysis of Image Distortions by Atmospheric Turbulence and Computer Simulation of Turbulence Effects[C]∥Proceedings of SPIE， 2008， 6941： 69410S-1-13.

[5] Repasi E， Weiss R. Computer Simulation of Image Degradations by Atmospheric Turbulence for Horizontal Views[C]∥Proceedings of SPIE， 2011， 8014： 80140U-1-9.

[6] Oxford D， Espinola R L. Simulation of a Laser RangeGated SWIR Imaging System in Weak Turbulence Conditions[C]∥Proceedings of SPIE， 2011， 8014： 80140T-1-12.

[7] Zhi Dong， Tao Rumao， Zhou Pu， et al. Propagation of Ring Airy Gaussian Beams with Optical Vortices through Anisotropic NonKolmogorov Turbulence [J]. Optics Communications， 2017， 387： 157-165.

[8] Wang F， Toselli I， Li J， et al. Measuring Anisotropy Ellipse of Atmospheric Turbulence by Intensity Correlations of Laser Light [J]. Optics Letters， 2017， 42（6）： 1129-1132.

[9] Yan Xu， Guo Lixin， Cheng Mingjian， et al. Probability Density of Orbital Angular Momentum Mode of Autofocusing Airy Beam Carrying PowerExponentPhase Vortex through Weak Anisotropic Atmosphere Turbulence [J]. Optics Express， 2017， 25（13）： 15286-15298.

[10] Baykal Y. Intensity Fluctuations of Asymmetrical Optical Beams in Anisotropic Turbulence[J]. Applied Optics， 2016， 55（27）： 7462-7467.

[11] Baykal Y， Luo Y J， Ji X L. Scintillations of Higher Order Laser Beams in Anisotropic Atmospheric Turbulence [J]. Applied Optics， 2016， 55（33）： 9422-9426.

[12] Wang F， Korotkova O. Random Optical Beam Propagation in Anisotropic Turbulence along Horizontal Links [J]. Optics Express， 2016， 24（21）： 24422-24434.

[13] Cui Linyan， Xue Bindang. Influence of Anisotropic Turbulence on the LongRange Imaging System by the MTF Model [J]. Infrared Physics & Technology， 2015， 72： 229-238.

[14] Cui Linyan. Analysis of Angle of Arrival Fluctuations for Optical Waves Propagation through Weak Anisotropic NonKolmogorov Turbulence[J]. Optics Express， 2015， 23（5）： 6313-6325.

[15] Toselli I， Korotkova O. General ScaleDependent Anisotropic Turbulence and Its Impact on Free Space Optical Communication System Performance [J]. Journal of the Optical Society of America A， 2015， 32（6）： 1017-1025.

[16] Bos J P， Roggemann M C， Gudimetla V S R. Anisotropic NonKolmogorov Turbulence Phase Screens with Variable Orientation [J]. Applied Optics， 2015， 54（8）： 2039-2045.

Abstract： The atmospheric turbulence degrades the image quality， which is a bottleneck for improving the performance of long range photoelectronic imaging detection. By making the research of image simulation in atmospheric turbulence， the performance of photoelectronic imaging detection system can be evaluated， and the test data for image processing algorithm verification and target recognition can be provided. In this work， the imaging simulation in anisotropic nonKolmogorov atmospheric turbulence is performed. By combining the image processing algorithm with the theoretical imaging model in anisotropic nonKolmogorov atmospheric turbulence， the rapid imaging simulation in anisotropic atmospheric turbulence considering physical parameters including anisotropic coefficient， imaging wavelength， receiver aperture and propagation path is realized. Compared with the simulation methods of random phase screen and pure image processing from simulation time and accuracy respectively， the results show that the proposed method costs shorter simulation time and achieves better imaging effects.

Key words： anisotropy；atmospheric turbulence； turbulent spectral index； image distortion； image blur