各向異性大氣湍流中成像仿真研究

崔林艷

摘 要： 大氣湍流降低成像質量， 一直以來是制約遠距離光電成像探測的瓶頸問題。 通過對大氣湍流中成像仿真進行研究， 一方面可以評估光電成像探測系統性能， 同時還可以為圖像處理算法驗證、 目標識別等提供測試數據。 本文主要針對各向異性非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