基于瀑布型多重网格加速的复指数波前复原算法*

赵旺 董理治 杨平 王帅† 许冰

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倍.本文所述方法在近地面激光大气传输校正、天文望远镜低仰角观测等领域具有潜在应用价值,后续将开展相关实验研究.