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

張冬曉 陳志斌 肖程 秦夢澤 吳浩

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

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

1 引 言

大氣湍流的運動導致大氣折射率的隨機擾動,當光波在湍流大氣中傳輸時,其波面會發生畸變,嚴重影響了光的傳輸[1].大氣湍流對光傳輸的影響主要由復振幅場的隨機方程進行描述,由于解析解求解困難,一般采用數值分析的方法[2?4].而湍流相位屏的模擬則是數值分析的核心.

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

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

為提高運算速度并兼顧非均勻采樣方法高精度的特點,本文提出了一種混合生成方法,在高頻區域沿用經典的FFT方法,而在低頻區域采用非等間隔的稀疏采樣方法.在低頻區域選取若干采樣點,以模擬屏相位結構函數和理論結構函數的相對誤差為目標函數,采用引力搜索算法對所選采樣點的位置進行優化,從而模擬出高精度相位屏.

2 低頻補償屏模型

本文提出的相位屏生成方法在高頻區域與經典的FFT方法一致,主要區別在于低頻補償屏的生成.低頻補償屏由若干采樣點代表的頻率成分組成,而低頻補償精度取決于低頻采樣點的位置和分布,為保證相位屏生成速度以及優化過程的收斂速度,在保證低頻補償精度的前提下應盡量減少低頻采樣點數和位置參數.通過大量的試驗后發現,當選取16個采樣點和兩個位置參數時能夠在保證高精度的同時提高相位屏的生成速度及優化收斂速度.

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

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

其中

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

圖1 低頻采樣點分布Fig.1.The distribution of low frequency samples.

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

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

3 優化算法及過程

設計優化算法之前,首先要確定優化的目標函數,在此選擇模擬相位屏結構函數與理論結構函數的相對誤差作為目標函數.

模擬相位屏的結構函數一般是通過統計大量的模擬相位屏實現的,但是這勢必引入龐大的計算量.在這里選擇文獻[20]使用的計算方法,即通過模擬相位屏的自相關函數計算其結構函數,模擬相位屏自相關函數由高頻屏和低頻補償屏的自相關函數和求得.高頻屏的自相關函數求解可參考文獻[20]中的(19)式,并注意將第一采樣點以內的頻譜置零,記為BφHigh.低頻補償屏的自相關函數可由下式求得:

由此可得模擬相位屏的結構函數為[21]

其 中Dtheory(r)為理論結構函數,由下式求得[21]:

式中r0為Fried參數,表示湍流的強弱;K5/6(·)為第三類修正分數貝塞爾函數;G(·)為Gamma函數.由于相對誤差為距離的函數,以此作為目標函數是不合適的,因此選擇相位屏尺寸范圍內的最大相對誤差作為優化目標函數:

考慮到高頻成分過采樣引入的部分高頻誤差會對小尺度結構函數產生影響,且不能被低頻屏補償,因此須將這部分點去除.根據文獻[25]有關高頻誤差的補償結果可知,高頻誤差影響的這部分點占整個尺度的5%左右[25],因此剔除這部分點后的最大相對誤差為

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

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

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

4 優化結果與討論

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

為驗證算法的穩定性,重復進行了50次優化,并統計最優適應值的均值和方差,其中均值μ=1.06×10?3,方差σ2=1.88×10?8.可見采用本文所提的優化算法能將相位屏的最大相對誤差降至0.1%左右,并且算法具有很強的穩定性.選取其中3個具有代表性的優化過程,并作出其優化曲線,如圖2所示.

圖2 引力搜索算法優化曲線Fig.2.The optimization curve of GSA.

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

從圖3所示的相位屏模擬結果可見,按照本文方法生成相位屏除了具有豐富的細節外還包含有充足的低頻部分,體現為相位屏的整體傾斜.為與經典的低頻補償方法作比較,我們在相同的參數條件下對Johansson 和Gavel[20]所提方法和本文方法分別進行了模擬,各生成1000幅相位屏,計算相位結構函數,并與(7)式所示的理論結構函數相比,得到了兩種方法生成相位屏的結構函數曲線以及相對誤差曲線,如圖4所示.

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

圖4 兩種方法的相位屏結構函數對比(a)結構函數曲線;(b)相對誤差曲線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)中可以看出,利用本文所提方法生成的相位屏其相位結構函數與理論結構函數曲線幾乎一致.從圖4(b)中可以看出,在低頻部分本文所提方法可將低頻區域的最大相對誤差降低至0.063%,明顯優于經典次諧波補償方法的 5%.而對于高頻部分,由于本文采用經典的FFT方法,因此對于這部分誤差并未做出補償,從圖中可知高頻部分誤差與經典次諧波方法是一致的.對高頻誤差的補償,Xiang[25]曾提出一種有效的方法,能夠將高頻誤差降至0.1%,可直接利用該方法對相位屏高頻誤差進行補償.

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

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

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

利用本文提出的優化算法對不同比值下的相位屏進行模擬,并計算其相位結構函數.其中在不同比值下優化得到的參數值c1和c2如表2所示.

表2 不同L0/L下的最優參數Table 2.The optimization parameters with different L0/L.

從表2中數據可以看出,對于不同的比值,低頻采樣點分布是不同的,并且還可看到隨著比值的增加,參數值之間差別越來越小.選取其中的5項作出其結構函數曲線以及理論結構函數曲線,如圖5所示.

由圖5所示曲線可看出,本文所提的優化算法能夠適用于不同情況下的相位屏模擬,并且模擬精度很高.此外,相位屏的生成速度也是評價相位屏生成算法的一項重要指標,特別是對于實時性要求高的場合,如自適應光學系統的測試等,更是如此.

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

圖5 不同L0/L下的相位屏結構函數曲線與理論結構函數曲線Fig.5.The expected structure functions vs.theoretical structure functions with different L0/L.

5 結 論

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

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