趙旺 董理治 楊平 王帥? 許冰

1)(中國科學院自適應(yīng)光學重點實驗室,成都 610209)

2)(中國科學院光電技術(shù)研究所,成都 610209)

3)(中國科學院大學,北京 100049)

1 引 言

激光在大氣中傳輸時,由于強湍流或長傳輸距離的影響,會出現(xiàn)強閃爍效應(yīng),此時波前中出現(xiàn)相位不連續(xù)點[1,2].常規(guī)波前復原算法不能準確復原波前中的不連續(xù)結(jié)構(gòu),使得自適應(yīng)光學系統(tǒng)校正效果下降甚至失效[3-5].為了解決該問題,Le等[6]開展了復指數(shù)因子波前復原算法(complex exponential reconstructor,CER)研究,其將相位差用復指數(shù)表示,利用迭代計算實現(xiàn)了對單個相位奇點的復原,該方法需要上萬次迭代才能收斂,實用性差.Fried[7]和Barchers等[8]將CER算法中迭代計算簡化為降采樣、最小二乘求解、數(shù)據(jù)重構(gòu)三個過程,實現(xiàn)了對CER算法的加速.但是該方法要求子孔徑數(shù)目滿足 2N×2N,子孔徑數(shù)目不滿足要求時,需要對子孔徑進行擴充,使得波前復原結(jié)果存在較大誤差.同時,通光口徑不是正方形時,由于引入過多的權(quán)值為零數(shù)據(jù)點,波前復原結(jié)果將存在較大殘差.此外,Aubailly 和Vorontsov[9]以及Yazdani和Fallah[10]將相位恢復算法用于相位奇點復原,此方法需要測量微透鏡后不同位置的光強數(shù)據(jù),利用光傳輸原理復原波前.其光強測量系統(tǒng)復雜,波前復原算法計算量大,難以應(yīng)用于自適應(yīng)光學系統(tǒng).

為了解決復指數(shù)波前復原算法計算量大的問題,本文提出了基于瀑布型多重網(wǎng)格(cascadic multigrid method,CMG)加速的復指數(shù)波前復原算法，該算法利用夏克-哈特曼波前傳感器測量相位差給出不同網(wǎng)格層上光場間的關(guān)系,在最粗網(wǎng)格上計算滿足收斂條件的解,并將其插值到細網(wǎng)格層,作為該層迭代計算的初值,減少了迭代計算次數(shù),加速了波前復原過程.本文利用CMG算法和直接迭代復原了不連續(xù)相位和隨機連續(xù)相位，比較了兩種波前復原算法的波前復原精度和算法復雜度，對比了采用CMG算法和直接斜率法的自適應(yīng)光學系統(tǒng)對大氣湍流像差的校正效果,并給出了相應(yīng)結(jié)論.

2 復指數(shù)波前復原算法

由傅里葉光學可知第i個子孔徑焦平面上光場復振幅[11]為

式中u(x0,y0) 為入射光場復振幅,λ為激光波長,f為微透鏡焦距,k=2π/λ.

焦平面上光強分布為

式中

子孔徑聚焦光斑x方向質(zhì)心位置為

由傅里葉變換微分性質(zhì)和卷積定理可得：

其中

利用(4)式可得

將(6)式代入(3)式可得x方向質(zhì)心位置為：

同理,子孔徑聚焦光斑y方向質(zhì)心位置為：

入射光波x,y方向波前斜率為：

包含相位奇點的波前[1]中,枝切線(branch cut)兩側(cè)存在 2nπ 跳變.假設(shè)某個子孔徑中測量相位φ(x)存在枝切線,波前被枝切線分割成兩部分,即

其中,φt(x)為不包含 2nπ 跳變的連續(xù)相位,S1和S2表示被枝切線分割的兩部分.依據(jù)復指數(shù)性質(zhì)有eiφ(x)=eiφt(x),由(7)式—(9)式可知,夏克-哈特曼波前傳感器不能探測波前中的 2nπ 相位跳變.

Hudgin模型[12]中,重建相位點在柵格點上,測量數(shù)據(jù)是柵格點間的相位差,如圖1所示,紅色實心點表示重建相位點,箭頭表示測量相位差.子孔徑大小為d時,第(i,j)個子孔徑x,y方向測量相位差為：

式中,A為 + 1,0,—1組成的稀疏矩陣,用于表示測量相位差和真實相位的關(guān)系.ΔφHS為向量形式表示的測量相位差,φr為向量形式表示的真實相位.

波前中不存在枝切線時,可以利用迭代計算求解(12)式的方程組,迭代計算公式[13]為：

圖1 測量相位差和重建相位點的關(guān)系Fig.1.The relationship between phase differences and phase.

存在枝切線的相位中,測量相位差和真實相位差滿足 ΔφHS=Δφr+2nπ,Δφr表示真實相位差.此時,(13)式不滿足等式關(guān)系,不能用其重建相位.由復指數(shù)性質(zhì),有 eiΔφHS=ei(Δφr+2nπ)=eiΔφr,測量相位差和真實相位差復指數(shù)相等.將(13)式中相位和相位差用復指數(shù)代替,相位和相位差間的加運算變成復指數(shù)乘運算,有

(14)式給出了復指數(shù)波前復原算法迭代計算過程,對(14)式所得結(jié)果取矢量輻角可得重建相位.復指數(shù)波前復原算法更為通用的迭代計算公式為：

迭代計算終止條件為

3 瀑布型多重網(wǎng)格法

瀑布型多重網(wǎng)格法[14](cascadic multigrid method,CMG)網(wǎng)格結(jié)構(gòu)和計算流程如圖2所示,相鄰兩層網(wǎng)格間距滿足hi=2hiˉ1.圖2(b)中紅色方框表示最粗網(wǎng)格上的迭代計算,圓點表示細網(wǎng)格上的迭代計算,箭頭表示插值過程.CMG算法計算流程為：首先在最粗網(wǎng)格上計算滿足收斂條件的解,將其插值到細網(wǎng)格,此數(shù)據(jù)作為細網(wǎng)格迭代計算的初值,重復迭代計算和插值過程直到在最細網(wǎng)格上得到滿足收斂條件的解.

圖2 瀑布型多重網(wǎng)格法示意圖 (a)網(wǎng)格結(jié)構(gòu)；(b)CMG算法計算流程Fig.2.Schematic of the CMG method：(a)Structure of network layers；(b)calculation process.

圖3 CMG算法降采樣過程 (a)細網(wǎng)格上光場；(b)粗網(wǎng)格上光場Fig.3.Downsampling process of the CMG method：(a)Data on the fine network；(b)data on the coarse network.

利用CMG算法加速(15)式時,夏克-哈特曼波前傳感器只給出了最細網(wǎng)格上測量相位差,粗網(wǎng)格相位差需通過細網(wǎng)格數(shù)據(jù)降采樣得到,降采樣過程如圖3所示.圖3(a)中實心點表示細網(wǎng)格上光場,圖3(b)中實心點表示粗網(wǎng)格上光場,圓圈表示剔除的細網(wǎng)格光場.粗網(wǎng)格上相位差復指數(shù)表示為細網(wǎng)格上相鄰復指數(shù)相乘,即

CMG算法插值過程如圖4所示,圖4(a)中實心點表示粗網(wǎng)格光場,空心點表示插值得到的細網(wǎng)格光場.圖4(b)中,u1,u2,u3,u4四個粗網(wǎng)格數(shù)據(jù)組成一個正方形,正方形中心點處光場通過多路加權(quán)平均得到,即

其中

圖4 CMG算法插值過程 (a)細網(wǎng)格光場和粗網(wǎng)格光場的關(guān)系；(b)待插值數(shù)據(jù)位于正方形中心；(c),(d)待插值數(shù)據(jù)位于正方形四邊上Fig.4.Interpolation process of the CMG method：(a)The relationship between grid points on coarse network and fine network；(b)the new grid point located at the center of the unit square；(c),(d)the new grid point located on the edge of the unit square.

式中,x1,x2,…,x6和y1,y2,…,y6為細網(wǎng)格上相位差復指數(shù),ωx1,ωx2,…,ωx6和ωy1,ωy2,…,ωy6為測量相位差權(quán)值,最細網(wǎng)格上相位差權(quán)值等于子孔徑光斑峰值信噪比,其他網(wǎng)格層上相位差權(quán)值由(18)式得到.

正方形四邊上待插值數(shù)據(jù)和已知光場關(guān)系如圖4(c)和圖4(d)所示,圖中實心環(huán)表示利用(19)式得到的光場,插值過程表示為：

其中,u1,u2,u3,u4為已知光場,x1,x2,y1,y2為細網(wǎng)格上相位差復指數(shù),ωx1,ωx2,ωy1,ωy2為細網(wǎng)格上相位差權(quán)值.

夏克-哈特曼波前傳感器子孔徑數(shù)目為N×N時,降采樣過程所需浮點乘數(shù)為

插值過程所需浮點乘數(shù)為

忽略網(wǎng)格邊界數(shù)據(jù)和內(nèi)部數(shù)據(jù)迭代計算的差異,單次迭代過程需要的浮點乘運算數(shù)目為

4 CMG算法波前復原性能

圖5(a)—圖5(d)給出了四種包含不同相位奇點的波前分布,分別用Phase1,Phase2,Phase3和Phase4表示.Phase1,Phase2中包含一個正相位奇點,Phase3中有兩個正相位奇點,Phase4包含一正一負兩個相位奇點.波前傳感器子孔徑數(shù)目等于20×20時,圖5(e)—圖5(h)給出了最小二乘法波前復原結(jié)果,圖5(i)—圖5(l)給出了復指數(shù)波前復原算法結(jié)果,仿真中波前傳感器探測噪聲為零,有效通光口徑內(nèi)的相位差權(quán)值為1,通光口徑外的相位差權(quán)值為0.由圖5中波前分布可知,最小二乘法不能重建相位奇點,利用復指數(shù)波前復原算法重建的相位中包含相位奇點和枝切線,但枝切線位置和輸入波前存在偏差.如文獻[15]和文獻[16]所述,枝切線為正負相位奇點的連線或相位奇點和邊線的連線,piston像差會使枝切線位置發(fā)生變化,實際系統(tǒng)中不考慮piston像差對復原結(jié)果的影響,去除piston像差后,復指數(shù)波前復原算法波前復原殘差為0.019λ,0.030λ,0.018λ和0.016λ.

圖5 (a)-(d)Phase1,Phase2,Phase3和Phase4二維分布；(e)-(h)最小二乘法波前復原結(jié)果；(i)-(l)復指數(shù)波前復原算法結(jié)果Fig.5.(a)-(d)Two-dimensional distribution of Phase1,Phase2,Phase3 and Phase4；(e)-(h)wavefront reconstructed by the leastsquares reconstruction algorithm；(i)-(l)wavefront reconstructed by the CER algorithm.

不同子孔徑數(shù)目時,直接迭代和CMG算法波前復原殘差如圖6(a)—圖6(c)所示.子孔徑數(shù)目等于20×20,40×40,80×80時,CMG算法分別選用3層、4層、5層網(wǎng)格加速迭代計算.如圖6所示,兩種算法波前復原殘差RMS值最大相差0.005λ,兩種波前復原算法都能有效復原帶有相位奇點的波前,且波前復原精度相近.

直接迭代和CMG算法復原Phase1,Phase2,Phase3和Phase4所需浮點乘運算數(shù)目如圖7(a)—圖7(c)所示.子孔徑數(shù)目為20×20時,直接迭代過程重建Phase1,Phase2,Phase3和Phase4所需浮點乘運算數(shù)約為106.此時,CMG算法大約需要2×104次浮點乘運算即可得到相近復原精度的結(jié)果,其所需浮點乘數(shù)相比直接迭代下降近2個數(shù)量級.子孔徑數(shù)目為80×80時,直接迭代過程所需浮點乘數(shù)的數(shù)量級為108,而CMG算法所需浮點乘數(shù)的數(shù)量級約為105,CMG算法所需浮點乘數(shù)相比直接迭代過程下降近3個數(shù)量級.

硬件條件為3.20 GHz Intel(R)Xeon(R)Gold 6134 CPU,128 G內(nèi)存的平臺上,直接迭代和CMG算法波前復原過程所需時間見表1.子孔徑數(shù)目為20×20,40×40,80×80時,CMG算法復原Phase1所需時間是直接迭代計算的1.9%,0.4%,0.6‰.相比直接迭代過程,CMG算法波前復原所需計算時間大大減少.隨著子孔徑數(shù)目增多,CMG算法引入更多的網(wǎng)格層用于加速波前復原過程,其在計算時間上更具優(yōu)勢,仿真測試結(jié)果和利用浮點乘數(shù)得到的結(jié)論一致.

利用前35階zernike多項式生成服從Kolmogorov統(tǒng)計規(guī)律的多組隨機像差檢驗CMG算法復原連續(xù)相位的性能[17].生成隨機像差時,D/r0等于10,D表示望遠鏡口徑,r0表示大氣相干長度.子孔徑數(shù)目等于20×20,40×40和80×80時,CMG算法和超松弛迭代法[18](successive over relaxation,SOR)波前復原殘差統(tǒng)計結(jié)果如圖8所示,圖中誤差線為標準差.不同子孔徑數(shù)目時,兩種波前復原算法復原精度沒有明顯差異.子孔徑數(shù)目等于20×20,40×40和80×80時,CMG算法所需浮點乘數(shù)目是SOR算法的30%,10%和3%.可見,復原連續(xù)相位時,同等復原精度下CMG算法所需計算量比SOR算法少.

圖6 直接迭代和CMG算法波前復原殘差 (a)子孔徑數(shù)目為20×20；(b)子孔徑數(shù)目為40×40；(c)子孔徑數(shù)目為80×80Fig.6.Wavefront residual error of the direct iteration method and the CMG method,the number of subapertures is (a)20×20；(b)40×40；(c)80×80.

圖7 CMG算法和直接迭代過程所需浮點乘運算數(shù)目(a)子孔徑數(shù)目為20×20；(b)子孔徑數(shù)目為40×40；(c)子孔徑數(shù)目為80×80Fig.7.Float point multiplications required by the CMG method and the process of the direct iteration the number of subapertures is (a)20×20；(b)40×40；(c)80×80.

綜上,CMG算法既可用于重建帶有相位奇點的波前,又可復原連續(xù)相位.重建相位奇點時,CMG算法所需浮點乘數(shù)相比直接迭代過程減少近2個數(shù)量級,且隨著子孔徑數(shù)目增加,CMG算法在計算量上更具優(yōu)勢.復原連續(xù)相位時,CMG算法復原精度與現(xiàn)有SOR算法相當,但所需浮點乘運算量更少.

表1 直接迭代和CMG算法波前復原時間(單位：s)Table 1.Time required by the direct iteration and CMG method (in s).

圖8 CMG算法和SOR算法波前復原殘差統(tǒng)計結(jié)果Fig.8.Wavefront residual statistics of the CMG method and SOR method.

5 采用CMG算法的自適應(yīng)光學系統(tǒng)校正能力

5.1 激光大氣傳輸補償系統(tǒng)數(shù)值仿真

多層相位屏法模擬激光在湍流大氣中傳輸時[19],將大氣湍流像差等效成相位屏,激光在兩張相位屏之間的傳輸過程用菲涅爾衍射描述.相位屏采用功率譜反演法[20]生成,大氣湍流功率譜選用von-Karman譜,湍流外尺度l0= 100 m.

數(shù)值計算中,望遠鏡口徑D= 600 mm,激光波長λ= 1064 nm,大氣相干長度r0= 6 cm,調(diào)節(jié)激光傳輸距離和大氣折射率結(jié)構(gòu)常數(shù)保證r0不變.自適應(yīng)光學系統(tǒng)中變形鏡驅(qū)動器和波前傳感器子孔徑的匹配關(guān)系如圖9所示,紅色圓圈表示驅(qū)動器位置,方格表示子孔徑位置,藍色圓圈表示有效通光口徑.變形鏡驅(qū)動器個數(shù)為20×20,有效驅(qū)動器個數(shù)為384.波前傳感器子孔徑數(shù)目為40×40,有效子孔徑數(shù)目為1240.變形鏡響應(yīng)函數(shù)交聯(lián)值ω= 0.08,高斯指數(shù)α= 2.2[21].

圖9 變形鏡驅(qū)動器和哈特曼波前傳感器子孔徑匹配關(guān)系Fig.9.Matching relation between actuators of deformable mirror and subapertures of Shack-Hartmann sensor.

CMG算法復原波前和驅(qū)動器控制電壓的關(guān)系為

其中Φ+表示變形鏡面形響應(yīng)函數(shù)廣義逆矩陣,φ為CMG算法復原波前.

激光光束質(zhì)量用峰值Strehl比表示,峰值Strehl比定義為實際光束遠場峰值光強與理想無像差光束遠場峰值光強之比[22],即

5.2 自適應(yīng)光學系統(tǒng)校正效果

不同Rytov方差時,自適應(yīng)光學系統(tǒng)校正前后遠場光強分布如圖10所示,圖上數(shù)字為峰值Strehl比.比較圖中數(shù)據(jù)可以發(fā)現(xiàn),經(jīng)自適應(yīng)光學系統(tǒng)校正后,遠場光強能量更加集中,采用CMG算法的自適應(yīng)光學系統(tǒng)校正效果優(yōu)于采用直接斜率法的系統(tǒng).圖10中不同Rytov方差的波前中包含相位奇點數(shù)目為50,98,120,126,136,隨著Rytov方差增大,波前中相位奇點數(shù)目增多,自適應(yīng)光學系統(tǒng)校正效果下降.

圖10 不同Rytov方差時,自適應(yīng)光學系統(tǒng)校正前后遠場光強分布及其峰值Strehl比Fig.10.Far field intensity and Strehl ratio of laser beam before and after corrected by the adaptive optics system.

圖11 不同Rytov方差時,自適應(yīng)光學系統(tǒng)校正光束Strehl比Fig.11.Strehl ratio of laser beam after corrected by the adaptive optics system in different Rytov number.

不同Rytov方差時,自適應(yīng)光學系統(tǒng)校正前后Strehl比平均值如圖11所示,圖中結(jié)果由20次仿真數(shù)據(jù)得到.由圖11數(shù)據(jù)可知CMG算法波前復原效果優(yōu)于直接斜率法,Rytov方差大于0.4時,自適應(yīng)光學系統(tǒng)采用CMG算法后校正光束Strehl比相比采用直接斜率法的系統(tǒng)提升1倍.出現(xiàn)這種差異的原因在于直接斜率法基于測量斜率和控制電壓滿足線性方程這一假設(shè)計算控制電壓,不能復原相位奇點,而CMG算法能夠重建相位奇點,從而提升自適應(yīng)光學系統(tǒng)校正效果.

6 結(jié) 論

本文提出了基于瀑布型多重網(wǎng)格加速的復指數(shù)波前復原算法,分析了最小二乘法不能復原相位奇點的原因.測量相位差和待復原相位點滿足Hudgin模型時,給出了CMG算法中降采樣、插值計算過程.同等復原精度下,相比直接迭代,CMG算法復原相位奇點所需浮點乘數(shù)下降近2個數(shù)量級,隨著子孔徑數(shù)目增加,其在計算量上的優(yōu)勢更加明顯.仿真結(jié)果表明,相比直接斜率法,自適應(yīng)光學系統(tǒng)采用CMG算法后校正光束Strehl比提升1倍.本文所述方法在近地面激光大氣傳輸校正、天文望遠鏡低仰角觀測等領(lǐng)域具有潛在應(yīng)用價值,后續(xù)將開展相關(guān)實驗研究.