大氣湍流等效相位屏的仿真研究*

徐瑞超，高 明

(西安工業(yè)大學(xué) 光電工程學(xué)院， 西安 710021)

光在大氣環(huán)境中傳播時，由于大氣溫度和壓強(qiáng)的變化會產(chǎn)生大氣折射率的隨機(jī)擾動，從而發(fā)生光束波前畸變.在激光通信領(lǐng)域，大氣湍流對激光傳輸?shù)挠绊憰鹜鈭龉怆娫O(shè)備的非正常工作，因此了解湍流的變化，掌握激光在大氣湍流中的傳輸特性對于激光通信中提高光束的光通量大有作用，所以需要對大氣湍流進(jìn)行相關(guān)方面的研究.在大氣湍流的相關(guān)研究中，通常采用數(shù)值方法對大氣湍流相位屏進(jìn)行模擬仿真，其具有清晰直觀的優(yōu)點.模擬大氣湍流等效相位屏的數(shù)值方法可以分為兩大類：① 功率譜反演法，由Mcglamery[1]提出，其是通過大氣湍流的功率譜密度函數(shù)得到擾動的大氣湍流相位分布.文獻(xiàn)[2]利用功率譜反演研究了一個完整的自適應(yīng)光學(xué)激光大氣傳輸?shù)南辔谎a(bǔ)償系統(tǒng)的數(shù)值模擬.文獻(xiàn)[3-4]從產(chǎn)生大氣湍流隨機(jī)相位屏的功率譜反演法原理出發(fā)，分析了均勻采樣造成的隨機(jī)相位屏大量低頻信息泄漏的不足，并且對大氣湍流功率譜非均勻采樣進(jìn)行了研究，認(rèn)為其可以有效改善傳統(tǒng)功率譜反演法低頻采樣嚴(yán)重不足的缺陷，實現(xiàn)高精度的大氣湍流相位屏的模擬；② 采用正交的澤尼克多項式法來模擬大氣湍流，這種方法對于大多形式的湍流功率譜研究較少，而在Kolmogrov 譜的模擬上效果比較顯著.文獻(xiàn)[5]利用澤尼克多項式對大氣湍流相位屏波前補(bǔ)償系統(tǒng)需要獨立的改正數(shù)進(jìn)行了解析描述.文獻(xiàn)[6]基于澤尼克多項式對大氣湍流波前相位畸變進(jìn)行了模擬.文獻(xiàn)[7]通過澤尼克多項式法仿真了大氣湍流相位屏并利用相位結(jié)構(gòu)函數(shù)進(jìn)行了驗證.以上的研究均使用單一方法對大氣湍流相位屏進(jìn)行了模擬，對于整個空間頻率大氣湍流適用哪種方法的研究較少.

因此，本文分別使用功率譜反演法和澤尼克多項式展開法模擬仿真了大氣湍流畸變波前相位屏，數(shù)值模擬了兩種方法產(chǎn)生的相位屏在不同參數(shù)下的相位結(jié)構(gòu)函數(shù)曲線，通過比較相位結(jié)構(gòu)函數(shù)與理論結(jié)果的差異來對比兩種方法在模擬大氣湍流相位屏不同空間頻率部分的優(yōu)缺點.

1 大氣湍流相位屏模擬方法

1.1 譜反演法

采用譜反演法模擬大氣湍流相位屏，首先假設(shè)一個復(fù)高斯隨機(jī)數(shù)矩陣，然后使用大氣湍流功率譜對其進(jìn)行濾波處理，最后通過傅里葉逆變換(Inverse Fast Fourier Transform,IFFT)得到大氣湍流中大氣擾動的相位分布[8].

假設(shè)有一零均值，其單位方差復(fù)值高斯隨機(jī)過程a(k)滿足

〈a(k)a(k′)〉=δ(k-k′)

(1)

式中：k,k′均為空間頻率；δ為Dirac函數(shù).進(jìn)一步假定a(k)是Hermitian對稱的，即

a(-k)=a*(k)

(2)

其中a*(k)為Hermitian共軛隨機(jī)過程.

用濾波函數(shù)G(k)來得到符合湍流大氣統(tǒng)計特性的等效相位屏S(r)，G(k)反映了相位屏上空間頻率為k的起伏分量幅度的期望值.S(r)的自相關(guān)函數(shù)為

exp(ik·r-ik′·r′)dkdk′

(3)

式中：r,r′均為空間任意點；S*(r′)為等效相位屏的共軛；G*(k′)為相位屏上空間頻率為k′的起伏分量幅度的共軛期望.

利用式(1)，將式(3)化簡為

(4)

由于〈S(r)S*(r′)〉是大氣湍流等效相位屏起伏的相關(guān)函數(shù)，所以根據(jù)Rytov理論可得復(fù)相位變換為

(5)

式中：z為湍流中任意點；L為空間傳輸距離；Φn(k)為大氣折射率功率譜密度函數(shù);ξ為比例系數(shù)；Δz為湍流薄層厚度.比較式(4)和式(5)可得

(6)

由于G(k)為一個實值非負(fù)函數(shù)，且和Φn(k)一樣各向同性，則

(7)

大氣折射率譜采用vonkarman譜，大氣折射率功率譜密度為

(8)

式中：Cn(h)為在斜程路徑傳輸?shù)耐牧鞔髿庵须S傳輸高度起伏的大氣結(jié)構(gòu)常數(shù)模型；h為接收機(jī)高度；k0,km均為譜參數(shù).

采用簡化模型，內(nèi)、外尺度對模擬結(jié)果的影響并未展開討論，所以內(nèi)、外尺度統(tǒng)一取為內(nèi)尺度l0=1 cm,外尺度L0=10 m.利用相關(guān)公式轉(zhuǎn)化為

(9)

其中r0為大氣相干長度.

應(yīng)用譜反演法生成隨機(jī)等效相位屏的公式為

(10)

利用譜反演法來模擬大氣湍流等效相位屏，算法步驟為

① 首先生成一個復(fù)隨機(jī)數(shù)矩陣k，其滿足高斯分布；

② 依據(jù)大氣折射率功率譜密度函數(shù)生成二維的功率譜密度函數(shù)矩陣Φ(k)；

(11)

④ 對矩陣Y(k)進(jìn)行傅里葉逆變換，從而得到其空間域形式，即

Y(l)=IFFTY(k)

(12)

其中l(wèi)為空間頻率.

⑤ 對于式(12)得到的矩陣Y(l)，可以將其分解為實部和虛部.其實部和虛部均可以代表一種隨機(jī)的大氣湍流相位屏S1和S2,表達(dá)式為

S1=ReY(l)，S2=ImY(l)

(13)

1.2 澤尼克多項式法

澤尼克多項式由無窮數(shù)量的多項式完全集合組成，具有兩個變量ρ和θ，其在單位圓內(nèi)部連續(xù)正交.澤尼克多項式只有在單位圓的內(nèi)部連續(xù)區(qū)域為正交，而在單位圓內(nèi)部的離散的坐標(biāo)上不具有正交性質(zhì).定義在單位圓上的澤尼克多項式極坐標(biāo)表達(dá)式[9]為

(14)

(15)

(16)

式中：m為多項式的角向級次；n為多項式的徑向級次.

澤尼克系數(shù)aj一般情況下被看作是具有零均值的Gauss 隨機(jī)變量，設(shè)其系數(shù)向量為

A=[a2,a3,…,ap]T

(17)

協(xié)方差陣為

C=E[AAT]

=E(a2,a2)E(a2,a3),…,E(a2,ap)

E(a3,a2)E(a3,a3),…,E(a3,ap)

E(ap,a2)E(ap,a3),…,E(ap,ap)

(18)

其中協(xié)方差為

(19)

式中：D為光學(xué)系統(tǒng)的口徑；Г[r]為伽馬函數(shù);δ為協(xié)方差參數(shù)；K為頻率特征因子；n,n′分別為澤尼克數(shù)Zi,Zj的徑向級次和角向級次，頻率特征因子K依賴于Zi,Zj.由式(19)可見，在統(tǒng)計上澤尼克系數(shù)間的關(guān)系是相關(guān)的.所以很有必要引入Karhunen-Loeve(K-L)函數(shù)，即可通過有確定方差的隨機(jī)量組合來表示大氣湍流的隨機(jī)波前，即

(20)

式中：Φ(r)為大氣湍流的隨機(jī)波前；bj為在統(tǒng)計獨立中的隨機(jī)系數(shù)；Kj(r)為K-L函數(shù).

矩陣B為

B=[b1b2…bj]T

(21)

K-L函數(shù)目前為止還沒有統(tǒng)一的解析式，但可以展開為Zernike多項式的形式，表達(dá)式為

(22)

式中：Zj(r)為澤尼克多項式；Vij為一般多項式.

將式(22)代入式(20)可得

(23)

比較式(19)和式(23)可見，矩陣A和矩陣B滿足關(guān)系式：

A=VB

(24)

其中V為酉陣.

由于矩陣C為厄米陣，所以存在酉陣U從而使得UCUT為對角陣.同理矩陣C可分解為

C=VSVT

(25)

其中矩陣S為對角陣.令B=UA，故

E[BBT]=E[UAATUT]

=UE[AAT]UT=S

(26)

A=U-1B=UTB=VB

(27)

綜上所述，應(yīng)用澤尼克多項式展開法來模擬大氣湍流等效相位屏，其重點在于系數(shù)aj的生成，生成步驟主要分為3步：

① 明確澤尼克多項式的階數(shù)，得到協(xié)方差矩陣C.

② 分解協(xié)方差矩陣C的對角陣S，得到C=VSVT，其中V為酉陣.

③ 生成一個零均值隨機(jī)向量，且其協(xié)方差矩陣C為對角陣S，從而得到多項式系數(shù)A=VB.

波前相位函數(shù)W(r,θ)可以應(yīng)用正交的澤尼克多項式表示為

(28)

2 兩種算法的模擬仿真比較

2.1 譜反演和澤尼克多項式的仿真實現(xiàn)

為了比較兩種算法生成相位屏，特別選取了一定的大氣湍流強(qiáng)度(使用大氣相干長度r0表示)、不同口徑D以及不同澤尼克多項式階數(shù)N的數(shù)值，用兩種算法分別模擬生成等效相位屏.將兩種方法生成的等效相位屏作初步比較，然后再用相位結(jié)構(gòu)函數(shù)進(jìn)行分析比較.

圖1(a)是譜反演法r0=0.01 m,口徑寬度為0.3 m的波前相位圖；圖1(b)是澤尼克多項式法r0=0.01 m，N=100，D=0.3 m的波前相位圖；圖1(c)是澤尼克多項式法r0=0.01 m，N=200，D=0.3 m的波前相位圖；圖1(d)是澤尼克多項式法r0=0.01 m，N=400，D=0.3 m的波前相位圖.

圖2(a)是譜反演法r0=0.01 m,口徑寬度為0.5 m的波前相位圖；圖2(b)是澤尼克多項式法r0=0.01 m,D=0.5 m,N=100的波前相位圖；圖2(c)是澤尼克多項式法r0=0.01 m,D=0.5 m,N=200的波前相位圖；圖2(d)是澤尼克多項式法r0=0.01 m,D=0.5 m,N=400的波前相位圖.

圖3(a)是譜反演法r0=0.1 m，口徑寬度為0.3 m的波前相位圖；圖3(b)是澤尼克多項式法r0=0.1 m,D=0.3 m,N=100的波前相位圖；圖3(c)是澤尼克多項式法r0=0.1 m,D=0.3 m,N=200的波前相位圖；圖3(d)是澤尼克多項式法r0=0.1 m,D=0.3 m,N=400的波前相位圖.

圖1 不同參數(shù)下兩種算法模擬生成的大氣湍流相位屏

圖2 不同參數(shù)下兩種算法模擬生成的大氣湍流相位屏

圖4(a)是譜反演法r0=0.1 m，口徑寬度為0.5 m的波前相位圖；圖4(b)是澤尼克多項式法r0=0.1 m,D=0.5 m,N=100的波前相位圖；圖4(c)是澤尼克多項式法r0=0.1 m,D=0.5 m,N=200的波前相位圖；圖4(d)是澤尼克多項式法r0=0.1 m,D=0.5 m,N=400的波前相位圖.相位結(jié)構(gòu)函數(shù)為

DΦ(r)=6.88(ρ/r0)5/3

(29)

式中：ρ為空間相干長度；r0為大氣相干長度.

圖3 不同參數(shù)下兩種算法模擬生成的大氣湍流相位屏

圖4不同參數(shù)下兩種算法模擬生成的大氣湍流相位屏

Fig.4 Phase screen of atmospheric turbulence generated by two algorithms with different parameters

2.2 相位結(jié)構(gòu)函數(shù)的比較評價

采用相位結(jié)構(gòu)函數(shù)DΦ(r)和大氣相干長度r0可以描述大氣湍流中大氣擾動的相位分布統(tǒng)計特性，對于 Kolgmgorov 譜而言，其相位結(jié)構(gòu)函數(shù)的確定方法見文獻(xiàn)[10].

圖5為在波長為 6.328×10-7m ，相位屏大小為 256×256 pixel，采樣點間隔為 0.003 m，相位屏間距為 500 m條件下，應(yīng)用譜反演法模擬大氣湍流相位屏?xí)r的相干長度r0=0.1 m，對應(yīng)大氣湍流結(jié)構(gòu)常數(shù)為2×10-15的相位屏，并對相位屏疊加不同級次(1，2，3，4級)的次諧波后得到的相位屏相位結(jié)構(gòu)函數(shù).

圖 6 為在波長為 6.328×10-7m，相位屏大小為 256×256 pixel，口徑為 0.6 m 的條件下,采用不同階次(5,100,200,400階)Zernike多項式模擬相干長度r0= 0.1 m(中等強(qiáng)度)的大氣湍流相位屏得到的相位結(jié)構(gòu)函數(shù).

a-理論；b-4次諧波；c-3次諧波 d-2次諧波； e-1次諧波；f-無諧波

a-理論；b-5階澤尼克多項式；c-100階澤尼克多項式；d-200階澤尼克多項式；e-400階澤尼克多項式

通過橫向比較圖1～4以及縱向比較圖5和圖6，可以看出：澤尼克多項式法產(chǎn)生的相位屏在低頻空間部分比較吻合，但是在高頻空間部分明顯不如譜反演法更貼近理論曲線；通過增加澤尼克多項式的系數(shù)可以改善高頻不足的現(xiàn)象，但不能完全消除影響；模擬仿真中應(yīng)用澤尼克多項式法得到仿真結(jié)果明顯慢于應(yīng)用譜反演法；應(yīng)用譜反演法產(chǎn)生大氣湍流等效相位屏，低頻部分會出現(xiàn)與理論值不吻合的現(xiàn)象，進(jìn)行次諧波補(bǔ)償產(chǎn)生較大計算量，而澤尼克多項式法相對來說具有計算量小的優(yōu)點.

3 結(jié) 論

通過譜反演法和澤尼克多項式法對大氣湍流相位屏進(jìn)行了模擬仿真比較，并利用相位結(jié)構(gòu)函數(shù)加以驗證，得出結(jié)論為

1) 通過譜反演法生成的相位屏的結(jié)構(gòu)函數(shù)高頻部分與理論值吻合，但存在低頻不足的現(xiàn)象，需要進(jìn)行低頻補(bǔ)償來消除反演誤差；

2) 澤尼克多項式法產(chǎn)生的相位屏的結(jié)構(gòu)函數(shù)與譜反演法相反，在低頻率空間部分其與理論值基本吻合，但在高頻率空間部分則誤差較大，可通過增加澤尼克多項式的階數(shù)來改善高頻不足的情況，隨著多項式階數(shù)的增加會引起計算量過大的問題.

3) 對于湍流的研究是為了能夠及時了解湍流相位在大氣中的分布，在實際研究工作中可以將兩種方法結(jié)合考慮，減小了使用單一方式而造成的數(shù)值仿真誤差.

參 考 文 獻(xiàn)：

[1] MCGLAMERY B L.Restoration of Turbulence-degraded Images[J].Journal of the Optical Society of America,1967,57(3):293.

[2] YAN H X,LI S S,ZHANG D L,et al.Numerical Simulation of an Adaptive Optics System with Laser Propagation in the Atmosphere[J].Applied Optics,2000,39(39):3023.

[3] 蔡冬梅,王昆,賈鵬,等.功率譜反演大氣湍流隨機(jī)相位屏采樣方法的研究[J].物理學(xué)報,2014,63(10):227.

CAI Dongmei,WANG Kun,JIA Peng,et al.Sampling Methods of Power Spectral Density Method Simulating Atmospheric Turbulence Phase Screen[J].Physics Journal,2014,63(10):227.(in Chinese)

[4] 蔡冬梅,遆培培,賈鵬,等.非均勻采樣的功率譜反演大氣湍流相位屏的快速模擬[J].物理學(xué)報,2015,64(22):248.

CAI Dongmei,TI Peipei,JIA Peng,et al.Fast Simulation of Atmospheric Turbulence Phase Screen Based on Non-uniform Sampling[J].Physics Journal,2015,64(22):248.(in Chinese)

[5] NOLL R J.Zernike Polynomials and Atmospheric Turbulence[J].Journal of the Optical Society of America,1976,66(3):207.

[6] RODDIER N A.Atmospheric Wavefront Simulation Using Zernike Polynomials[J].Optical Engineering,1990,29(10):1174.

[7] 王奇濤,佟首峰,徐友會.采用Zernike多項式對大氣湍流相位屏的仿真和驗證[J].紅外與激光工程,2013,42(7):1907.

WANG Qitao,TONG Shoufeng,XU Huiyou.On Simulation and Verification of the Atmospheric Turbulent Phase Screen with Zernike Polynomials[J].Infrared and Laser Engineering,2013,42(7):1907.

(in Chinese)

[8] BAHR G V.Investigations into the Spherical and Chromatic Aberrations of the Eye and Their Influence on Its Refraction[J].Acta Ophthalmologica,1945,23(1):1.

[9] 劉良清.Matlab 輔助激光分析與應(yīng)用[M].武漢:武漢凌云光電科技有限公司,2008.

LIU Liangqing.Matlab Assisted Laser Analysis and Application[M].Wuhan:Wuhan Lingyun Optoelecronic Technology Cooperation,2008.(in Chinese)

[10] 王立瑾,李強(qiáng),魏宏剛,等.大氣湍流隨機(jī)相位屏的數(shù)值模擬和驗證[J].光電工程,2007,34(3):1.

WANG Lijin,LI Qiang,WEI Honggang,et al.Numerical Simulation and Validation of Phase Screen Distorted by Atmospheric Turbulence[J].Opto-electronic Engineering,2007,34(3):1.(in Chinese)