基于引力搜索算法的湍流相位屏生成方法*

張冬曉 陳志斌 肖程 秦夢(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ù)振幅場(chǎng)的隨機(jī)方程進(jìn)行描述,由于解析解求解困難,一般采用數(shù)值分析的方法[2?4].而湍流相位屏的模擬則是數(shù)值分析的核心.

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

FFT方法有一個(gè)顯著的缺陷,即對(duì)低頻成分的過(guò)采樣,使得生成的相位屏缺乏低頻部分,這是制約其模擬精度的主要因素[18].為解決低頻不足的問(wèn)題,Herman 和Strugala[18]提出了利用次諧波補(bǔ)償?shù)姆椒ǜ纳艶FT方法低頻不足的問(wèn)題,其主要思想是在低頻區(qū)域利用更小的采樣間隔來(lái)實(shí)現(xiàn)對(duì)功率譜更精確的模擬,從而獲得更多的低頻信息,但該方法對(duì)低頻成分的補(bǔ)償有限并且會(huì)增加大量的計(jì)算.此后,越來(lái)越多的學(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]提出了一種稀疏譜模型,即通過(guò)對(duì)低頻部分的密集采樣和高頻部分的稀疏采樣來(lái)提高相位屏的低頻模擬精度,獲得了較好的結(jié)果.2014年,蔡冬梅等[23,24]同樣提出了利用非均勻采樣的方法生成相位屏,結(jié)果表明模擬屏的相位結(jié)構(gòu)函數(shù)與理論結(jié)構(gòu)函數(shù)符合得很好.但是上述兩種方法由于采用了非均勻采樣,無(wú)法直接使用快速傅里葉變換,計(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)化過(guò)程的收斂速度,在保證低頻補(bǔ)償精度的前提下應(yīng)盡量減少低頻采樣點(diǎn)數(shù)和位置參數(shù).通過(guò)大量的試驗(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方向以及過(guò)原點(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)化算法及過(guò)程

設(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ù)一般是通過(guò)統(tǒng)計(jì)大量的模擬相位屏實(shí)現(xiàn)的,但是這勢(shì)必引入龐大的計(jì)算量.在這里選擇文獻(xiàn)[20]使用的計(jì)算方法,即通過(guò)模擬相位屏的自相關(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ù):

考慮到高頻成分過(guò)采樣引入的部分高頻誤差會(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è)體的過(guò)程,此時(shí)個(gè)體質(zhì)量用于評(píng)價(jià)個(gè)體的優(yōu)劣,質(zhì)量最大個(gè)體占據(jù)最優(yōu)位置.算法主要包括三個(gè)部分: 群體初始化、個(gè)體質(zhì)量和引力計(jì)算以及個(gè)體位置速度更新,具體步驟可參考文獻(xiàn)[26].

需要說(shuō)明的是,原始的引力搜索算法存在求解精度不高和收斂速度慢的問(wèn)題,并不適用于實(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)過(guò)大量試驗(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)化過(guò)程發(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.可見(jiàn)采用本文所提的優(yōu)化算法能將相位屏的最大相對(duì)誤差降至0.1%左右,并且算法具有很強(qiáng)的穩(wěn)定性.選取其中3個(gè)具有代表性的優(yōu)化過(guò)程,并作出其優(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é)果可見(jiàn),按照本文方法生成相位屏除了具有豐富的細(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é)開(kāi)始曾提到,相位屏生成過(guò)程中,除了兩個(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è)值得探討的問(wèn)題,為此對(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無(wú)關(guān),即當(dāng)模擬條件僅涉及大氣湍流強(qiáng)度改變時(shí),無(wú)需對(duì)低頻采樣點(diǎn)進(jìn)行調(diào)整; 對(duì)于不同的采樣點(diǎn)數(shù)N,由于低頻采樣點(diǎn)分布的不匹配導(dǎo)致誤差有增大的趨勢(shì),說(shuō)明低頻采樣點(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)過(guò)分析得到,當(dāng)兩者取得一致的相對(duì)誤差時(shí),湍流外尺度L0與相位屏尺寸L的比值L0/L是相等的,因此可以預(yù)見(jiàn),低頻采樣點(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ù)值之間差別越來(lái)越小.選取其中的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í)性要求高的場(chǎng)合,如自適應(yīng)光學(xué)系統(tǒng)的測(cè)試等,更是如此.

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

圖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)度無(wú)關(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í)性要求高的場(chǎng)合.