基于多面循環Gerchberg-Saxton算法相位復原研究

彭金錳，李志鵬

(中國人民解放軍69213部隊21分隊，喀什844900)

引 言

光在傳輸過程中，受傳播介質、光學儀器誤差、光學元件熱或力影響下的變形等因素影響，光波會產生畸變，這些因素綜合作用的結果就是成像模糊不清、分辨率降低［1］。高能激光系統中，受大氣湍流帶來的波前畸變影響，激光光束在遠處無法有效聚焦，不能得到接近衍射極限的輸出激光光束，激光在湍流大氣中傳輸時，湍流將對它產生各種效應，如光強閃爍、相位畸變、光斑擴展和漂移等，這些效應均會導致激光光束質量的嚴重退化［2］，因此在實際的光學系統中，常需要測量畸變波前，為自適應校正等應用提供依據。

Gerchberg-Saxton(GS)算法在20世紀70年代由GERCHBERG和 SAXTON共同提出［3-4］，通過測量光瞳面和成像面的光強分布，利用迭代算法進行相位恢復。在實際應用中，GS算法可能會面臨著光瞳面光強無法測量的情況，必須利用其它兩個或多個位置處的光強信息，首先恢復這些位置的相位信息，然后利用光傳輸理論得到光瞳處的相位和振幅分布。在利用得到的相位信息計算光瞳的光場時，往往會增加誤差，因此，必須提高焦面或離焦面的相位恢復精度。

作為一種局部搜索算法，GS算法極易陷入局部極值，而不能得到理想的精度［5-7］，尤其是在被恢復的相位比較復雜時。為了得到更精確的解，GU等人利用中間某兩次的迭代得到的相位值進行加權作為下次迭代的初始條件，得到了更為準確的結果［8］。HUANG等人利用了3個面(1個輸入面和2個輸出面)的光強信息，采用與GS算法類似的迭代算法，并將當次的相位迭代值與上次的迭代值之差作為梯度方向，在下次循環中引入了該方向的一個附加值，大幅提升了收斂的速度［9］。本文中提供了一種擴展的循環GS算法，增加一個面的光強信息，使算法有機會跳出極值，達到比較高的精度。

1 光場傳輸理論

假定平面j的相位是 θj(j=1，2，…)，振幅分布為Fj且已知，則其光場為:

式中，r=1，2，3…，M;s=1，2，3…，N;M 和 N 分別為平面j兩個方向的采樣點數。使用角譜傳輸的方法，第 k個面的光場［10］為:

式中，DFT(discrete Fourier transform)和 IDFT(inverse discrete Fourier transform)分別表示正反離散傅里葉變換，λ為波長，t和u為傅里葉變換系數，dr和ds為第j個面的橫向和縱向采樣間距，與第k個面一致。Δzkj為第k個面到第j個面的傳輸距離。由第j個面計算第k個面光場，公式為:

透鏡后的光場可以由焦平面或離焦面光場通過(3)式計算，而通過透鏡后，振幅可以認為不變，只是增加了一個相位因子θ:

式中，f表示透鏡焦距，p和q表示采樣點坐標。

2 算法描述

光路中任意位置的兩個或多個面光強分布已知的情況下，都可以利用GS算法得出這些面的相位分布，進而計算出其它位置的光場。本文中選擇了焦平面和兩個離焦面光強，如圖1所示，對光瞳處的光場進行恢復。算法流程如下:(1)給定平面1的初始嘗試相位φ0，與平面1的振幅結合;(2)利用(2)式計算平面2的光場，保留相位信息，與平面2的振幅結合;(3)利用(3)式，計算平面1的光場，保留相位信息，與平面1的振幅結合;(4)循環第2步至第3步至給定步數或評價函數達到給定值。

Fig.1 Schematic diagram of GS algorithm

只利用兩個面的光強信息時，往往迭代幾十次后收斂，得到的精度比較低，增加迭代次數并不能解決問題。為了能跳出局部極值，達到理想精度，這里提供了一種簡單的循環方式，利用了3個面的光強信息進行迭代，如圖2所示。首先利用平面1和平面2進行迭代，收斂后利用得到平面2的相位作為初始值，與平面2的振幅結合，進入平面2與平面3之間的迭代，得到新的平面2的相位信息。循環上述過程幾次。最后以計算得到的平面2的光場，利用(3)式和(4)式逆向計算光瞳處的光場。

Fig.2 Schematic diagram of circulatory GS algorithm

3 改進算法仿真結果與分析

3.1 高斯光束恢復

透鏡前的光場振幅為高斯分布，波前由4階至10階ZERNIKE多項式［11］生成，系數分別為 2.9，0.72，-1.2，0.72，-0.96，-0.72，0.72 和 -0.24，光瞳區域內波長畸變約為0.3λ，以后的例子中，光瞳處波前分布一致。波長為1.064μm，透鏡焦距500mm，光瞳直徑6mm，采樣點數為512×512。定義振幅誤差R為:

焦平面和距離焦平面150mm和300mm處得到的振幅分布，如圖4所示。

Fig.3 Amplitude and phase of round pupil plane a—amplitude b—phase

Fig.4 Amplitude of focal plane and two defocused planes a—focal plane b—150mm c—300mm

以零相位面作為算法的初始值，采用改進GS算法，最終恢復出來的透鏡前的振幅和分布(見圖5)，圖5c表示復原出的波前殘差，單位為λ，后文中同此表示。

Fig.5 Retrieved amplitude，phase and residual phase error of round pupil planea—amplitude b—phase c—phase error

恢復出的波前誤差為 0.007λ，如果不考慮0.4mm的邊緣區域，光強在峰值的0.035以內區域精度達到了0.002λ。

盡管透鏡前的光場未知，但在仿真中仍然可以監視光瞳平面的恢復出的振幅與實際振幅分布的誤差，如圖6所示。

Fig.6 Amplitude error of plane 2 and round pupil plane vs.iteration numbera—plane2 b—pupil plane

前100次迭代即為GS算法，由圖6a可以看出，GS算法迭代幾十次后左右陷入極值，復原出的振幅誤差較大，復原出的波前分布幾乎無參考價值。這種情況在進行迭代的光場波前相對復雜時更容易發生。改進后，進入下次循環，算法跳出了局部極值，最初平面2振幅誤差大幅增加，隨后迅速下降，收斂到更小的值。而透鏡前的振幅誤差一直下降，表明復原出的光場更接近實際分布。

3.2 圓形平頂光束恢復

對圓形平頂光束進行相位恢復時，GS算法更容易陷入局部極值，在經過10次循環后，最終波前誤差為0.0049λ?；謴徒Y果如圖7所示。

與高斯光束波前恢復相比，需要更多循環次數，而且振幅和相位恢復結果弱于高斯光束恢復效果。

Fig.7 a—real amplitude of pupil plane b—retrieved amplitude of pupil plane c—phase error of pupil plane d—amplitude error of plane2 vs.iteration number e—amplitude error of pupil plane vs.iteration number

Fig.8 Amplitude and phase of annulus pupil plane a—amplitude b—phase

Fig.9 Retrieved amplitude，phase and phase error of annulus pupil plane a—amplitude b—phase c—phase error

3.3 圓環光束恢復

光瞳處的圓環內徑1mm，外徑4mm。光束如圖8所示?；謴统龅恼穹筒ㄇ皻埐钊鐖D9所示，最終波前誤差為0.0064λ。

在對內徑為2mm、外徑4mm的圓環光束進行恢復時，算法經過多次循環后精度未有提升，算法在多循環中，只在兩個局部極值間變動，不能進一步提高精度，這與首次循環迭代得到的結果精度不高有關。

4 結論

在光瞳處振幅未知的情況下，采用用GS算法，利用焦平面和附近一個平面上的光強信息對光瞳平面上的光場進行復原，誤差較大，以至于無參考價值。增加一個平面后，可以比較準確地重建光瞳光場，事實上，每一處的光場都可以得到準確復原，優于只利用兩個面光強信息的迭代得到的結果。信息量的增加，使得算法能夠有效地跳出局部極值，達到更準確的結果。

［1］WANG X，ZHAO D Z，MAO H，et al.The developmentof phase diversity［J］.Optical Technique，2009，35(3):454-460(in Chinese).

［2］GE X L，FENG X X，FAN Ch Y.Progress of the study of phase discontinuity of laser propagation through atmosphere［J］.Laser Technology，2012，36(4):485-489(in Chinese).

［3］GERCHBERG RW，SAXTON W O.Phase determination from image and diffraction plane pictures in electron microscope［J］.Optik，1971，34(3):237-246.

［4］GERCHBERG RW，SAXTONW O.A practical algorithm for determination of phase from image and diffraction plane pictures［J］.Optik，1972，35(2):275-284.

［5］FIENUP JR.Phase retrieval algorithms:a comparison［J］.Applied Optics，1982，21(15):2758-2769.

［6］FIENUP JR.Phase-retrieval algorithms for a complicated Optical system［J］.Applied Optics，1993，32(10):1737-1746.

［7］RODDIER C，RODDIER F.Combined approach to the hubble space telescopewave-front distortion analysis［J］.Applied Optics，1993，32(15):2992-3008.

［8］GU X，XU K Sh.Some improvements on the GS algorithm in phase retrieval problem［J］.Journal of Fudan University，2000，33(2):205-211(in Chinese).

［9］HUANG L X，YAO X，CAID M，et al.A high accuracy and fast iterative algorithm for phase retrieval［J］.Chinese Journal of Lasers，2010，37(5):1218-1221(in Chinese).

［10］BRADY G R，FIENUP JR.Nonlinear optimization algorithm for retrieving the full complex pupil function［J］.Optics Express，2006，14(2):474-486.

［11］NOLL R J.Zernike polynomial and atmospheric turbulence［J］.Journal of the Optical Society of America，1976，A66(3):207-211.