基于瀑布型多重網格加速的復指數波前復原算法*

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

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

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

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

1 引 言

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

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

2 復指數波前復原算法

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

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

焦平面上光強分布為

式中

子孔徑聚焦光斑x方向質心位置為

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

其中

利用(4)式可得

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

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

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

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

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

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

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

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

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

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

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

迭代計算終止條件為

3 瀑布型多重網格法

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

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

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

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

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

其中

圖4 CMG算法插值過程 (a)細網格光場和粗網格光場的關系；(b)待插值數據位于正方形中心；(c),(d)待插值數據位于正方形四邊上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為細網格上相位差復指數,ωx1,ωx2,…,ωx6和ωy1,ωy2,…,ωy6為測量相位差權值,最細網格上相位差權值等于子孔徑光斑峰值信噪比,其他網格層上相位差權值由(18)式得到.

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

其中,u1,u2,u3,u4為已知光場,x1,x2,y1,y2為細網格上相位差復指數,ωx1,ωx2,ωy1,ωy2為細網格上相位差權值.

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

插值過程所需浮點乘數為

忽略網格邊界數據和內部數據迭代計算的差異,單次迭代過程需要的浮點乘運算數目為

4 CMG算法波前復原性能

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

圖5 (a)-(d)Phase1,Phase2,Phase3和Phase4二維分布；(e)-(h)最小二乘法波前復原結果；(i)-(l)復指數波前復原算法結果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.

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

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

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

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

圖6 直接迭代和CMG算法波前復原殘差 (a)子孔徑數目為20×20；(b)子孔徑數目為40×40；(c)子孔徑數目為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算法和直接迭代過程所需浮點乘運算數目(a)子孔徑數目為20×20；(b)子孔徑數目為40×40；(c)子孔徑數目為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算法既可用于重建帶有相位奇點的波前,又可復原連續相位.重建相位奇點時,CMG算法所需浮點乘數相比直接迭代過程減少近2個數量級,且隨著子孔徑數目增加,CMG算法在計算量上更具優勢.復原連續相位時,CMG算法復原精度與現有SOR算法相當,但所需浮點乘運算量更少.

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

圖8 CMG算法和SOR算法波前復原殘差統計結果Fig.8.Wavefront residual statistics of the CMG method and SOR method.

5 采用CMG算法的自適應光學系統校正能力

5.1 激光大氣傳輸補償系統數值仿真

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

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

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

CMG算法復原波前和驅動器控制電壓的關系為

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

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

5.2 自適應光學系統校正效果

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

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

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

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

6 結 論

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