張冬曉 陳志斌 肖程 秦夢(mèng)澤 吳浩

1)(陸軍工程大學(xué)石家莊校區(qū)電子與光學(xué)工程系,石家莊 050003)

2)(32181部隊(duì),石家莊 050003)

1 引 言

大氣湍流的運(yùn)動(dòng)導(dǎo)致大氣折射率的隨機(jī)擾動(dòng),當(dāng)光波在湍流大氣中傳輸時(shí),其波面會(huì)發(fā)生畸變,嚴(yán)重影響了光的傳輸[1].大氣湍流對(duì)光傳輸?shù)挠绊懼饕蓮?fù)振幅場的隨機(jī)方程進(jìn)行描述,由于解析解求解困難,一般采用數(shù)值分析的方法[2?4].而湍流相位屏的模擬則是數(shù)值分析的核心.

目前,對(duì)于相位屏的模擬主要有基于快速傅里葉變換(FFT)的蒙特卡羅方法[5]、Zernike多項(xiàng)式法[6,7]、協(xié)方差法[8?10]、分形插值法[11?13]以及一些混合方法[14,15].這幾種相位屏生成方法中,基于FFT的蒙特卡羅方法(簡稱FFT方法)由于具有快速計(jì)算的優(yōu)勢(shì),被廣泛用于湍流相位屏的模擬,特別是空間光通信、大氣激光傳輸以及遙感成像系統(tǒng)的仿真研究等要求實(shí)時(shí)性的場合[16,17].

FFT方法有一個(gè)顯著的缺陷,即對(duì)低頻成分的過采樣,使得生成的相位屏缺乏低頻部分,這是制約其模擬精度的主要因素[18].為解決低頻不足的問題,Herman 和Strugala[18]提出了利用次諧波補(bǔ)償?shù)姆椒ǜ纳艶FT方法低頻不足的問題,其主要思想是在低頻區(qū)域利用更小的采樣間隔來實(shí)現(xiàn)對(duì)功率譜更精確的模擬,從而獲得更多的低頻信息,但該方法對(duì)低頻成分的補(bǔ)償有限并且會(huì)增加大量的計(jì)算.此后,越來越多的學(xué)者對(duì)次諧波補(bǔ)償方法進(jìn)行了改進(jìn),例如,Lane等[19]利用多級(jí)次諧波網(wǎng)格實(shí)現(xiàn)低頻補(bǔ)償,既提高了模擬精度又減小了計(jì)算量; Johansson 和Gavel[20]將低頻區(qū)域進(jìn)行了擴(kuò)展,使得低頻誤差降至5%; Sedmak[21]提出了一種加權(quán)的次諧波方法,同時(shí)補(bǔ)償了相位屏的高低頻不足,將總體誤差降低到1%,但是計(jì)算量增加了大約兩倍.從上述方法的次諧波網(wǎng)格劃分可以看出,它們對(duì)低頻成分的采樣仍是均勻等間隔的,對(duì)于功率譜衰減極快的低頻部分顯然是不合適的.2013年,Charnotskii[22]提出了一種稀疏譜模型,即通過對(duì)低頻部分的密集采樣和高頻部分的稀疏采樣來提高相位屏的低頻模擬精度,獲得了較好的結(jié)果.2014年,蔡冬梅等[23,24]同樣提出了利用非均勻采樣的方法生成相位屏,結(jié)果表明模擬屏的相位結(jié)構(gòu)函數(shù)與理論結(jié)構(gòu)函數(shù)符合得很好.但是上述兩種方法由于采用了非均勻采樣,無法直接使用快速傅里葉變換,計(jì)算速度較慢.

為提高運(yùn)算速度并兼顧非均勻采樣方法高精度的特點(diǎn),本文提出了一種混合生成方法,在高頻區(qū)域沿用經(jīng)典的FFT方法,而在低頻區(qū)域采用非等間隔的稀疏采樣方法.在低頻區(qū)域選取若干采樣點(diǎn),以模擬屏相位結(jié)構(gòu)函數(shù)和理論結(jié)構(gòu)函數(shù)的相對(duì)誤差為目標(biāo)函數(shù),采用引力搜索算法對(duì)所選采樣點(diǎn)的位置進(jìn)行優(yōu)化,從而模擬出高精度相位屏.

2 低頻補(bǔ)償屏模型

本文提出的相位屏生成方法在高頻區(qū)域與經(jīng)典的FFT方法一致,主要區(qū)別在于低頻補(bǔ)償屏的生成.低頻補(bǔ)償屏由若干采樣點(diǎn)代表的頻率成分組成,而低頻補(bǔ)償精度取決于低頻采樣點(diǎn)的位置和分布,為保證相位屏生成速度以及優(yōu)化過程的收斂速度,在保證低頻補(bǔ)償精度的前提下應(yīng)盡量減少低頻采樣點(diǎn)數(shù)和位置參數(shù).通過大量的試驗(yàn)后發(fā)現(xiàn),當(dāng)選取16個(gè)采樣點(diǎn)和兩個(gè)位置參數(shù)時(shí)能夠在保證高精度的同時(shí)提高相位屏的生成速度及優(yōu)化收斂速度.

如圖1所示,將低頻區(qū)域限定在高頻網(wǎng)格第一采樣點(diǎn)所圍成的正方形區(qū)域內(nèi),并將區(qū)域內(nèi)(包括邊緣部分)所包含的原始采樣點(diǎn)置零; 接著向低頻區(qū)域添加16個(gè)采樣點(diǎn),從圖中可看出,此時(shí)低頻區(qū)域被原始采樣點(diǎn)自然地劃分為四個(gè)子區(qū)域,分別處于頻率坐標(biāo)的四個(gè)象限,考慮到湍流功率譜的對(duì)稱性,四個(gè)象限分別包含相同個(gè)數(shù)的采樣點(diǎn)并且各象限的采樣點(diǎn)對(duì)稱分布,因此只用考慮一個(gè)象限的4個(gè)采樣點(diǎn)分布即可.

在各象限外邊緣的x,y方向以及過原點(diǎn)的對(duì)角線方向分別設(shè)置三個(gè)可活動(dòng)的采樣點(diǎn),其位置由參數(shù)f1和f2決定(f1,f2∈[0,1/L]),注意到x方向和y方向的兩個(gè)點(diǎn)由同一個(gè)參數(shù)決定.在各象限的角點(diǎn)處分別設(shè)置固定采樣點(diǎn),其坐標(biāo)為(±1/L,±1/L).至此,所有采樣點(diǎn)的空間頻率坐標(biāo)已經(jīng)確定,根據(jù)稀疏譜模型,低頻補(bǔ)償屏可表示為

其中

為Von-Karman譜的功率譜密度,fxjk和fyjk為圖1所示采樣點(diǎn)的空間頻率坐標(biāo),?fxjk和?fyjk為各采樣點(diǎn)的頻率間隔,Lx和Ly為相位屏在x,y方向的尺寸,Nx和Ny分別為相位屏在x,y方向的采樣點(diǎn)數(shù),m和n為整數(shù),f0=1/L0為與湍流外尺度L0相關(guān)的空間頻率,h(j,k)為零均值、單位方差的復(fù)高斯隨機(jī)矩陣.

圖1 低頻采樣點(diǎn)分布Fig.1.The distribution of low frequency samples.

至此,得到了低頻補(bǔ)償屏模型,最終的相位屏由低頻補(bǔ)償屏和高頻相位屏相加得到,由于高頻相位屏生成方法采用經(jīng)典的FFT方法,所以這里直接給出高頻相位屏模型:

其中 ?fx=1/Lx,?fy=1/Ly,Φ(m′?fx,n′?fy)如(2)式所示.低頻補(bǔ)償屏由參數(shù)f1和f2共同決定,不同的參數(shù)值產(chǎn)生不同的采樣點(diǎn)分布,進(jìn)而生成的補(bǔ)償屏精度也不同,利用優(yōu)化算法可對(duì)參數(shù)進(jìn)行優(yōu)化,使得生成的相位屏結(jié)構(gòu)函數(shù)與理論結(jié)構(gòu)函數(shù)之間的相對(duì)誤差達(dá)到最小,此時(shí)所生成的相位屏精度最高.

3 優(yōu)化算法及過程

設(shè)計(jì)優(yōu)化算法之前,首先要確定優(yōu)化的目標(biāo)函數(shù),在此選擇模擬相位屏結(jié)構(gòu)函數(shù)與理論結(jié)構(gòu)函數(shù)的相對(duì)誤差作為目標(biāo)函數(shù).

模擬相位屏的結(jié)構(gòu)函數(shù)一般是通過統(tǒng)計(jì)大量的模擬相位屏實(shí)現(xiàn)的,但是這勢(shì)必引入龐大的計(jì)算量.在這里選擇文獻(xiàn)[20]使用的計(jì)算方法,即通過模擬相位屏的自相關(guān)函數(shù)計(jì)算其結(jié)構(gòu)函數(shù),模擬相位屏自相關(guān)函數(shù)由高頻屏和低頻補(bǔ)償屏的自相關(guān)函數(shù)和求得.高頻屏的自相關(guān)函數(shù)求解可參考文獻(xiàn)[20]中的(19)式,并注意將第一采樣點(diǎn)以內(nèi)的頻譜置零,記為BφHigh.低頻補(bǔ)償屏的自相關(guān)函數(shù)可由下式求得:

由此可得模擬相位屏的結(jié)構(gòu)函數(shù)為[21]

其 中Dtheory(r)為理論結(jié)構(gòu)函數(shù),由下式求得[21]:

式中r0為Fried參數(shù),表示湍流的強(qiáng)弱;K5/6(·)為第三類修正分?jǐn)?shù)貝塞爾函數(shù);G(·)為Gamma函數(shù).由于相對(duì)誤差為距離的函數(shù),以此作為目標(biāo)函數(shù)是不合適的,因此選擇相位屏尺寸范圍內(nèi)的最大相對(duì)誤差作為優(yōu)化目標(biāo)函數(shù):

考慮到高頻成分過采樣引入的部分高頻誤差會(huì)對(duì)小尺度結(jié)構(gòu)函數(shù)產(chǎn)生影響,且不能被低頻屏補(bǔ)償,因此須將這部分點(diǎn)去除.根據(jù)文獻(xiàn)[25]有關(guān)高頻誤差的補(bǔ)償結(jié)果可知,高頻誤差影響的這部分點(diǎn)占整個(gè)尺度的5%左右[25],因此剔除這部分點(diǎn)后的最大相對(duì)誤差為

采用引力搜索算法(gravitational search algorithm,GSA)對(duì)低頻采樣點(diǎn)的兩個(gè)未知參數(shù)進(jìn)行優(yōu)化,該算法由Rashedi等[26]于2009年首次提出,是利用個(gè)體間的引力作用尋找質(zhì)量最大個(gè)體的過程,此時(shí)個(gè)體質(zhì)量用于評(píng)價(jià)個(gè)體的優(yōu)劣,質(zhì)量最大個(gè)體占據(jù)最優(yōu)位置.算法主要包括三個(gè)部分: 群體初始化、個(gè)體質(zhì)量和引力計(jì)算以及個(gè)體位置速度更新,具體步驟可參考文獻(xiàn)[26].

需要說明的是,原始的引力搜索算法存在求解精度不高和收斂速度慢的問題,并不適用于實(shí)際應(yīng)用,為此需要對(duì)基本引力搜索算法進(jìn)行改進(jìn),受微粒群優(yōu)化算法中記憶功能的啟發(fā)[27],為基本引力搜索算法的速度和位置更新方程添加記憶項(xiàng),使得算法能夠在最優(yōu)解附近不產(chǎn)生劇烈的振蕩.此時(shí)速度和位置更新方程變?yōu)閇28]:

其中c1s,c2s表示c1和c2的迭代初始值,c1e,c2e表示c1和c2的迭代終值,tmax為最大迭代次數(shù).

4 優(yōu)化結(jié)果與討論

根據(jù)上節(jié)的分析,采用(9)式作為引力搜索算法的適應(yīng)度函數(shù),經(jīng)過大量試驗(yàn),最終確定算法各參數(shù)設(shè)置為: 樣本數(shù)M=50 ; 迭代次數(shù)T=100 ;引力常數(shù)G(t)=G0e?αt/T隨迭代次數(shù)變化,其中G0=100,α=50 ;c1s=1 ,c2s=5 ;c1e=2 ,c2e=2.由(6)式和(7)式可知,適應(yīng)度函數(shù)中的固定參數(shù)有r0,L0以及r,它們與大氣湍流狀態(tài)以及需要生成的相位屏尺寸以及采樣點(diǎn)數(shù)有關(guān),不會(huì)隨優(yōu)化過程發(fā)生變化.而(4)式中的fxjk,?fxjk和fyjk,?fyjk由第二節(jié)所述的參數(shù)f1和f2決定,兩個(gè)參數(shù)的值即(11)式中的二維位置參數(shù)將隨迭代次數(shù)逐漸趨于最優(yōu)值.因此在進(jìn)行優(yōu)化前,首先根據(jù)需要設(shè)置固定參數(shù)r0,L0,L與N,其中L為相位屏尺寸,N為采樣點(diǎn)數(shù)(假設(shè)相位屏為方形屏),為與傳統(tǒng)的次諧波方法作比較,采用與文獻(xiàn)[20]一致的參數(shù):r0=0.1m ,L0=10m ,L=2m ,N=256.

為驗(yàn)證算法的穩(wěn)定性,重復(fù)進(jìn)行了50次優(yōu)化,并統(tǒng)計(jì)最優(yōu)適應(yīng)值的均值和方差,其中均值μ=1.06×10?3,方差σ2=1.88×10?8.可見采用本文所提的優(yōu)化算法能將相位屏的最大相對(duì)誤差降至0.1%左右,并且算法具有很強(qiáng)的穩(wěn)定性.選取其中3個(gè)具有代表性的優(yōu)化過程,并作出其優(yōu)化曲線,如圖2所示.

圖2 引力搜索算法優(yōu)化曲線Fig.2.The optimization curve of GSA.

從圖2可以看出,算法前期能夠快速地進(jìn)入局部搜索,到后期在局部搜索空間中逐步逼近最優(yōu)值,其中最佳適應(yīng)值(最大相對(duì)誤差)僅為6.34 ×10?4,此時(shí)參數(shù)c1和c2的取值分別為6.438和9.048,生成的相位屏如圖3所示.

從圖3所示的相位屏模擬結(jié)果可見,按照本文方法生成相位屏除了具有豐富的細(xì)節(jié)外還包含有充足的低頻部分,體現(xiàn)為相位屏的整體傾斜.為與經(jīng)典的低頻補(bǔ)償方法作比較,我們?cè)谙嗤膮?shù)條件下對(duì)Johansson 和Gavel[20]所提方法和本文方法分別進(jìn)行了模擬,各生成1000幅相位屏,計(jì)算相位結(jié)構(gòu)函數(shù),并與(7)式所示的理論結(jié)構(gòu)函數(shù)相比,得到了兩種方法生成相位屏的結(jié)構(gòu)函數(shù)曲線以及相對(duì)誤差曲線,如圖4所示.

圖3 湍流相位屏模擬結(jié)果(a)相位屏二維分布;(b)相位屏三維分布Fig.3.A realization of turbulence phase screen:(a)Two dimensional distribution;(b)three dimensional distribution.

圖4 兩種方法的相位屏結(jié)構(gòu)函數(shù)對(duì)比(a)結(jié)構(gòu)函數(shù)曲線;(b)相對(duì)誤差曲線Fig.4.Expected structure functions generated by Johansson’s method and our method,where the theoretical structure function is shown for reference:(a)Phase structure functions(b)relative errors.

從圖4(a)中可以看出,利用本文所提方法生成的相位屏其相位結(jié)構(gòu)函數(shù)與理論結(jié)構(gòu)函數(shù)曲線幾乎一致.從圖4(b)中可以看出,在低頻部分本文所提方法可將低頻區(qū)域的最大相對(duì)誤差降低至0.063%,明顯優(yōu)于經(jīng)典次諧波補(bǔ)償方法的 5%.而對(duì)于高頻部分,由于本文采用經(jīng)典的FFT方法,因此對(duì)于這部分誤差并未做出補(bǔ)償,從圖中可知高頻部分誤差與經(jīng)典次諧波方法是一致的.對(duì)高頻誤差的補(bǔ)償,Xiang[25]曾提出一種有效的方法,能夠?qū)⒏哳l誤差降至0.1%,可直接利用該方法對(duì)相位屏高頻誤差進(jìn)行補(bǔ)償.

本節(jié)開始曾提到,相位屏生成過程中,除了兩個(gè)優(yōu)化參數(shù)外還有三個(gè)固定參數(shù)r0,L0以及r,它們與大氣湍流狀態(tài)(r0表示大氣湍流強(qiáng)度,L0表示湍流外尺度)、相位屏尺寸以及采樣點(diǎn)個(gè)數(shù)相關(guān).而低頻采樣點(diǎn)分布是否與三個(gè)固定參數(shù)相關(guān),也即是否與不同的湍流強(qiáng)度r0、外尺度L0以及相位屏尺寸L和采樣點(diǎn)個(gè)數(shù)N相關(guān)是一個(gè)值得探討的問題,為此對(duì)每個(gè)具體參數(shù)取7個(gè)不同的值,考察低頻采樣點(diǎn)分布與每個(gè)參數(shù)之間的相關(guān)性時(shí)保持其余參數(shù)不變,將各參數(shù)值分別代入相位屏生成模型中,計(jì)算各取值下的結(jié)構(gòu)函數(shù)最大相對(duì)誤差,結(jié)果如表1所示.

由表1可知,低頻采樣點(diǎn)的分布與參數(shù)r0無關(guān),即當(dāng)模擬條件僅涉及大氣湍流強(qiáng)度改變時(shí),無需對(duì)低頻采樣點(diǎn)進(jìn)行調(diào)整; 對(duì)于不同的采樣點(diǎn)數(shù)N,由于低頻采樣點(diǎn)分布的不匹配導(dǎo)致誤差有增大的趨勢(shì),說明低頻采樣點(diǎn)分布與采樣點(diǎn)數(shù)N有關(guān),但從表中可知,除去小采樣點(diǎn)數(shù)(32和64點(diǎn))下誤差有明顯變化外,其他采樣點(diǎn)數(shù)下誤差增加并不明顯,在實(shí)際應(yīng)用中可忽略采樣點(diǎn)數(shù)的影響.而誤差在小采樣點(diǎn)數(shù)時(shí)的明顯變化,是因?yàn)?6個(gè)低頻采樣點(diǎn)在小采樣點(diǎn)數(shù)情況下占的比重很大,使得在大采樣點(diǎn)數(shù)中的少許偏差被放大了; 而低頻采樣點(diǎn)分布與相位屏尺寸和湍流外尺度的關(guān)系較為復(fù)雜,采用同樣的分布模擬不同的尺寸、不同湍流外尺度的相位屏?xí)r,其誤差明顯增大.但是注意到表中L0=2,4,5,10時(shí)的相對(duì)誤差與L=10,5,4,2時(shí)的相對(duì)誤差是一致的,經(jīng)過分析得到,當(dāng)兩者取得一致的相對(duì)誤差時(shí),湍流外尺度L0與相位屏尺寸L的比值L0/L是相等的,因此可以預(yù)見,低頻采樣點(diǎn)分布與L0/L相關(guān),相同的比值采用相同的低頻采樣點(diǎn)分布.

表1 不同參數(shù)及參數(shù)值下的最大相對(duì)誤差Table 1.The maximum relative errors with different parameters.

利用本文提出的優(yōu)化算法對(duì)不同比值下的相位屏進(jìn)行模擬,并計(jì)算其相位結(jié)構(gòu)函數(shù).其中在不同比值下優(yōu)化得到的參數(shù)值c1和c2如表2所示.

表2 不同L0/L下的最優(yōu)參數(shù)Table 2.The optimization parameters with different L0/L.

從表2中數(shù)據(jù)可以看出,對(duì)于不同的比值,低頻采樣點(diǎn)分布是不同的,并且還可看到隨著比值的增加,參數(shù)值之間差別越來越小.選取其中的5項(xiàng)作出其結(jié)構(gòu)函數(shù)曲線以及理論結(jié)構(gòu)函數(shù)曲線,如圖5所示.

由圖5所示曲線可看出,本文所提的優(yōu)化算法能夠適用于不同情況下的相位屏模擬,并且模擬精度很高.此外,相位屏的生成速度也是評(píng)價(jià)相位屏生成算法的一項(xiàng)重要指標(biāo),特別是對(duì)于實(shí)時(shí)性要求高的場合,如自適應(yīng)光學(xué)系統(tǒng)的測(cè)試等,更是如此.

對(duì)于本文所提方法,當(dāng)優(yōu)化參數(shù)確定后,相位屏的生成速度由兩部分影響,一是二維快速傅里葉變換速度,二是低頻補(bǔ)償屏的離散傅立葉變換速度.對(duì)比經(jīng)典的FFT方法,同時(shí)也是最快的生成方法來說,本文方法引入的額外時(shí)間主要由低頻補(bǔ)償屏的離散傅里葉變換決定,對(duì)于本文設(shè)置的16個(gè)采樣點(diǎn)來講,其引入的計(jì)算復(fù)雜度為O(16·N2),作為對(duì)比,Johansson 和Gavel[20]所提方法至少需要2層次諧波甚至當(dāng)模擬無限外尺度湍流時(shí),需達(dá)到10層,引入的計(jì)算復(fù)雜度為O(36·p·N2),其中p為次諧波層數(shù).因此利用本文所提方法生成相位屏的速度至少是次諧波方法的4.5倍,可應(yīng)用于實(shí)時(shí)性要求高的場合.

圖5 不同L0/L下的相位屏結(jié)構(gòu)函數(shù)曲線與理論結(jié)構(gòu)函數(shù)曲線Fig.5.The expected structure functions vs.theoretical structure functions with different L0/L.

5 結(jié) 論

本文提出了一種全新的相位屏生成方法,借助改進(jìn)的引力搜索算法,在相位屏的低頻采樣區(qū)尋找最優(yōu)的采樣分布,彌補(bǔ)了經(jīng)典FFT方法低頻采樣不足的缺點(diǎn).與Johansson和Gavel[20]所提的次諧波補(bǔ)償方法相比,本文提出的方法能將相位屏在低頻區(qū)域的最大相對(duì)誤差降低至 0.063%.除此之外,本文還對(duì)相位屏的各項(xiàng)參數(shù)與低頻采樣點(diǎn)分布之間的關(guān)系進(jìn)行了分析,結(jié)果表明低頻采樣點(diǎn)分布與湍流強(qiáng)度無關(guān),因此在遇到湍流強(qiáng)度變化的情形時(shí),可直接模擬出相應(yīng)的相位屏而不需要重新進(jìn)行優(yōu)化; 相位屏的采樣點(diǎn)數(shù)與低頻采樣點(diǎn)分布相關(guān)性很小,不同采樣點(diǎn)數(shù)間的誤差變化并不明顯,因此可利用小采樣點(diǎn)進(jìn)行優(yōu)化,減小優(yōu)化時(shí)間; 低頻采樣點(diǎn)分布與湍流外尺度和相位屏尺寸之間的比值直接相關(guān),需要根據(jù)不同的比值對(duì)低頻采樣點(diǎn)分布進(jìn)行優(yōu)化.

對(duì)相位屏的生成速度進(jìn)行了討論,結(jié)果表明,利用本文所提方法生成相位屏的速度至少是次諧波方法的4.5倍,可用于實(shí)時(shí)性要求高的場合.