一種基于角譜理論的改進型相位恢復迭代算法*

劉宏展 紀越峰

1)(北京郵電大學，信息光子學與光通信國家重點實驗室，北京 100876)

2)(華南師范大學，廣東省微納光子功能材料與器件重點實驗室，廣州 510006)

(2012年12月20日收到;2013年1月29日收到修改稿)

1 引言

相位恢復是指利用光的衍射理論，對輸入光場進行衍射計算，得到輸出面光場的場強分布，將其與實測(或理想)的輸出場強數據進行比較，以能量轉換效率最大、誤差最小為準則，通過迭代或者搜索找到最符合實測(或理想)場強數據的相位分布.相位恢復是物理及工程中的一個基礎性問題[1]，由于其在信號恢復、空間光通信、光學衍射元件的設計等場合有廣泛的應用[2-4]，它已成為一個很重要的研究領域，其核心是要找到合適的相位恢復算法.而采用迭代算法進行相位恢復是當前主要的研究思路之一，并被運用到實踐中[5].早在1972年，Gerehberg等[6]提出了G-S迭代算法，G-S算法簡單而實用，但它的誤差并不隨著迭代次數增加而遞減，相對于其他算法而言最小誤差偏大，所得結果是相對最優.在此基礎上，發展了許多改進算法，例如輸入輸出算法(IO)、混合輸入輸出算法(HIO)[7]等，可它們并不能保證迭代過程總收斂到正確解，有時甚至會停滯在某個局部極小值附近;另一方面，混合遺傳-模擬退火算法[8]、免疫遺傳算法[9]和蟻群算法[10]等也相繼產生，但共同的缺點就是原理相對比較復雜，編程難度較大，且收斂速度相對較慢.以上算法都只適用于么正變換系統，基于此，楊國禎和顧本源提出任意線性變換系統中振幅-相位恢復的一般理論和楊-顧(Y-G)算法[11]，它已成功應用到各種實際問題和各種變換系統中.

本文結合星間衛星通信光學發射天線光束整形的實際需求，在傳統G-S算法的基礎上，提出了一種基于角譜傳播理論的、能對復雜光場進行快速高精度相位恢復的幅度梯度加成迭代算法(下文簡稱加成算法)，即利用角譜傳播理論，構建輸入和輸出面光場之間的正、逆向衍射過程，把整形要求達到的理想(或實際)輸出光場振幅通過幅度參數加入到迭代過程的中間環節，并引入梯度方法，綜合運用輸入、輸出光場信息，找到能量集中度高、誤差小的相位分布.由于角譜理論嚴格滿足亥姆霍茲方程，用它來處理迭代過程中輸入、輸出面光場間的衍射計算，可以得到更精確、可靠的結果[12];理想輸出光強加入到迭代中間環節，使得迭代一次時便綜合了輸入、輸出面的光強信息，不像傳統的G-S那樣只進行簡單的替換，從而加速了迭代的速度;梯度法的引入，通過自適應的加大迭代步長，確保了對復雜光場相位恢復的良好收斂性.文章首先對算法中將要用到的角譜理論進行簡短介紹，然后著重對新算法進行了描述與分析，最后對算法的性能，通過數值仿真進行了詳細的比較.

2 算法描述與分析

要進行相位恢復迭代運算，就必須進行光場衍射及逆衍射運算，因此我們將首先簡單介紹文中要用到的角譜理論，然后重點對加成算法進行描述與分析.

2.1 角譜理論

在直角坐標系中，令輸入面平面坐標為xi，yi，輸入光場振幅為U(xi,yi)、相位為φ(xi,yi);經某一空間距離d衍射后，在輸出面的坐標為xo，yo，輸出光場振幅為U(xo,yo)、相位為φ(xo,yo)，輸入、輸出光場的復振幅為Ei，Eo.引用傅里葉變換及逆變換符號：F{}，F-1{}，則由輸入光場，獲得輸出光場的角譜衍射積分為[13]

H?(fx,fy)是 H(fx,fy)的共軛，k=2π/λ，λ 為光波長，ΔLh為衍射光場的計算寬度，m，n均為采樣點數，N為采樣總數.由(1)，(2)式進行離散數值計算時，可借助快速傅里葉變換(FFT)完成.衍射光場的角譜理論嚴格滿足亥姆霍茲方程，用它來處理相位恢復迭代算法過程中輸入輸出面光場間的衍射計算，可以得到更精確、可靠的結果，這為算法的設計提供了理論支撐.

2.2 幅度梯度加成迭代算法

傳統G-S算法具有收斂速度較慢，對初始相位敏感等缺點，為了盡量消除這些不利因素，必須對傳統G-S算法進行改進，本文把改進的算法稱為幅度梯度加成迭代算法.其流程如圖1所示.

其中W 為輸出面整形光場的范圍，βk，gk，γk的表達式是參考文獻[14]得到的.

圖1 幅度梯度加成迭代算法流程圖

與傳統G-S算法相比，文中提出的加成算法做了三方面的改進：一方面，通過變量參數α把實際(或理想)輸出光場振幅與角譜衍射得到的輸出光場振幅共同組成新的輸出光場振幅，而不是像傳統G-S算法那樣，僅僅把角譜衍射得到的輸出光場振幅用實際(或理想)輸出光場振幅來替換.通過這樣的處理，讓每次的迭代衍射結果疊加在輸出面光場信息中，并通過逆衍射過程，反饋回輸入光場，從而構成了反饋環路，加速了迭代速度，具體見算法流程第4)步.另一方面，通過MATLAB中的自帶angle()函數，直接對光場求對應的相位，避免了某些情況下，因出現分母取零奇點而迭代算法無法進行的情況，具體見算法流程第4)，6)步.第三方面，利用文獻[14]的思想，把一種簡單的梯度方法引用到迭代過程，在保證算法向目標最優解迭代的過程中，通過加大每次迭代的步長，加速了整個算法的收斂速度，且進行迭代時，不需要預先知道迭代值之外的其他信息，計算過程簡單，具體見算法流程的第8)步.

3 加成算法仿真與結果分析

在星間激光通信系統中，因其通信距離達到幾萬公里，要求光束以近衍射極限發散角進行發射，加上有效星載的限制，光學天線多采用同軸兩鏡式反射望遠鏡結構[2].由于次鏡對光束中心部分的遮擋，嚴重減損光束能量，因此，在光束進入發射天線前，必須對其進行有效的整形，使其變換成空心環狀激光束，保證光束能量能夠盡可能多地從光學天線發射出去.下邊就以設計合適的星間激光通信星載光學天線光束衍射整形器為例子，來檢驗新提出的加成算法的功能.考慮到功率受限對星載系統的影響，光束整形應盡量保持光束的能量，因此具體迭代過程以能量集中度η為判斷依據.

3.1 參數設置

具體參數如下，光波長λ=830 nm，整形后的高斯光束要求整形成外徑R=3 mm，內徑r=1.5 mm，即遮擋比為R/r=2的中空光束，從而減少光學天線中的次鏡對光束能量的遮擋.高斯光束的束腰為0.75 mm，假設光束整形器放置在距入射光束束腰d=450 mm處，入射、出射光場所在平面的范圍都設定為半徑為5 mm的圓形區域，取樣點數N=256.有了這些參數，利用(1)—(9)式，就上文中提出的加成算法，即可對其性能進行深入探討.

3.2 仿真與分析

為了突出加成算法的優點，本節將從以下三方面著手：首先，從迭代速度以及誤差下降速度與G-S算法進行比較;其次，通過初始相位取不同隨機值，證明加成算法對初始相位是否具有適應性;最后，通過變換能量集中度η指標，證明算法是否具有收斂一致性.

先考查加成算法的迭代速度.從算法的流程圖1可以看到，必須先確定幅度變量參數α，這里采用參考文獻[15]中的試探法來實現，根據文中給出的數據，得到α=0.8.令輸入光束的初始相位φ(xi,yi)=0，則新算法與G-S算法的相關數據如表1所示，表中列出了η分別為90%，94%時，兩種算法的迭代次數與最小誤差的取值.從表1可見，當η=90%，G-S算法的迭代次數為22次，加成算法為16次，加成算法的迭代速度略快，但隨著能量集中度指標的提升，當η=94%時，GS算法要355次，而加成算法為213次，可見，加成算法的收斂速度明顯優于G-S算法.另外，從最小誤差來看，見表中第3、5列，似乎G-S算法能達到的最小誤差的絕對值比加成算法小，但其實，它們都處在同一數量級上，因此從這個角度來看，它們能達到的最小誤差是基本一致的.同樣還是考察最小誤差，若從單位迭代次數引起誤差的下降速率來看，加成算法明顯比G-S算法快.比如，利用表1數據，加成算法的誤差下降速率=(0.1520-0.0950)/(213-16)=3.12×10-4，同理，G-S的下降速率=1.83×10-4，兩者相比，得到加成算法單位迭代次數引起誤差的下降速率是G-S的1.7倍.因此，在達到同樣級別誤差的約束條件下，加成算法收斂速度明顯優于G-S算法.

表1 G-S算法與幅度梯度加成算法比較

由于G-S算法對輸入光場的初始相位敏感，在進行迭代運算時，通常需要經過多次試探后，才能找到合適的初始相位值，進行迭代運算，這樣費時費力.接下來，我們研究加成算法對初始相位的適應性.首先選定能量集中度η指標為90%，然后令輸入初始相位φ(xi,yi為取值在[0 π]之間的隨機數，具體由MATLAB中的rand()函數產生，進行一次完整的加成算法運算，得到當前隨機初始相位下的迭代總次數M，并記錄;重新產生新的隨機初始相位φ(xi,yi)，進行新一輪加成算法運算，得到新的M值;依次重復10次，所得結果如2所示.

圖2 不同隨機初始相位下迭代次數的分布

圖2中，橫坐標表示初始相位取10次不同值時，對應的加成算法運行次序n，縱坐標表示加成算法運行一次的迭代總次數M.對于確定的能量集中度，從圖2可以看到，當輸入光場的初始相位隨機變化時，加成算法的最小迭代次數為10次，最大為12次，絕大多數都是11次.因此，對于不同的隨機初始相位值，加成算法表現出極好的適應性，它可以克服G-S算法對初始相位敏感的缺陷.不僅如此，與零初始相位相比，見表1中迭代次數的數據，初始相位取隨機取值時，加成算法具有更快的收斂速度.

對于相位恢復迭代算法，都希望對于不同的隨機初始相位值，算法最后都能收斂于最小的誤差點.為此，令能量集中度η指標分別為90%，94%，分10次，對隨機初始相位取不同的值，研究加成算法的收斂一致性，其仿真結果如圖3所示.

圖3 不同隨機初始相位下迭代誤差的分布

圖3中，有兩條圖線，橫坐標與圖2的含義一致，縱坐標表示每次加成算法能達到的最小誤差e.上邊那條對應η=90%，10次運算的結果是，其迭代誤差e密集分布在0.14兩側，即迭代誤差趨向于0.14，可見，加成算法的迭代收斂具有一致性.同時，與初始相位為0時e=0.1520(見表1)相比，其誤差有所減小.下邊的圖線對應η=94%，其迭代誤差e更加密集地分布于0.09附近，與η=90%相比，此時迭代誤差更加減小，且收斂一致性更優.這也從一個側面說明，加成算法中用能量集中度η或者迭代誤差e作為判斷依據是統一的，這是因為，能量更集中了，自然迭代所能達到的最小誤差便變小;而迭代誤差變小了，則光束能量就更集中.

3.3 仿真實例

最后，設定能量集中度η以達到92%為標準，運用加成算法得到的星間激光通信光學天線光束整形器如圖4所示.從圖4(a)看到，迭代過程迭代誤差e是單調遞減的，能量集中度η是單調遞增的，沒有出現吉布斯現象[16]，而理想輸出光場圖4(b)與迭代所得光場圖4(c)的誤差e=0.1150.考慮到衍射光場邊界的限制，在文中給定數據的約束下，當迭代誤差進步一下降后，能量集中度η可以達到95%，因此，從能量利用率的角度來看，基于加成迭代算法進行光學天線光束整形器設計是可行的.圖4(d)是所設計光束整形器的相位分布圖，它決定了光束整形的最終效果，根據其分布，通過相關的運算、處理，運用衍射光學技術，就可以制作出來，其具體的器件實施過程，這里將不作討論.

圖4 仿真實例 (a)迭代誤差、能量集中度與迭代次數關系圖;(b)理想的整形光束幅度;(c)迭代所得整形幅度;(d)衍射整形元件的相位分布

4 結論

基于光場的角譜傳播理論，在傳統G-S相位恢復迭代算法的基礎上，本文提出了一種幅度梯度加成相位恢復迭代算法.該算法通過引入幅度變量參數α，把輸入光場的幅度信息通過正向衍射過程，加入到輸出光場幅度的表征當中，并通過逆衍射過程，饋送至輸入光場，形成反饋回路，從而加速迭代的收斂速度，而不像傳統G-S算法那樣只是單純的把正向衍射得到的光場振幅用實際(或理想)的輸出光場振幅置換，在整個迭代過程中，輸入或輸出光場的振幅都是單向傳播，不能形成閉合回路，不利于加快迭代過程的收斂.另外，相對于傳統的G-S，文中提出的加成算法在進行新一輪迭代過程時，引入了梯度運算，自適應地加大了迭代步長，從而也加快了收斂速度.從仿真結果看，在約束能量集中度η后，對于不同的隨機初始相位，加成算法的最終迭代次數M趨同，可見，加成算法的收斂性并不依賴于初始相位分布，具有優良的適應性;與零初始相位相比，隨機初始相位的迭代次數和最小誤差都更小，這更加驗證了加成算法對隨機初始相位的適應性，而且收斂速度更快，這些都得益于幅度反饋與梯度的聯合作用.不僅如此，對于設定的不同能量集中度標準，輸入不同隨機初始相位，新算法每次得到的最終迭代最小誤差都趨于同一數值，比如，當η=94%時，最小誤差e都密集分布在0.09附近，加成算法表現出優良的收斂一致性.綜上所述，幅度梯度加成相位恢復迭代算法，具有收斂速度快，對初始相位適應性好，收斂一致性好等優點，它為復雜光場的相位恢復提供了一種新的方法，對設計衍射光學元件、光學整形器件，尤其是對能量集中度要求高的整形器件，提供了技術支持.

[1]Huang L X，YaoX，CaiD M，GuoY K，YaoJ，GaoF H 2010 Chin.J.Laser.37 1218(in Chinese)[黃利新，姚新，蔡冬梅，郭永康，姚軍，高福華2010中國激光37 1218]

[2]Yu J J，Tan L Y，Ma J，Han Q Q，Yang Y Q，LiM 2009 Chin.J.Laser.36 581(in Chinese)[俞建杰，譚立英，馬晶，韓琦琦，楊玉強，李密2009中國激光36 581]

[3]Wu R，Hua N，Zhang X B，CaoG W，ZhaoD F，Zhou S L 2012 Acta Phys.Sin.61 224202(in Chinese)[鄔融，華能，張曉波，曹國威，趙東峰，周申蕾2012物理學報61 224202]

[4]Yu B，Peng X，Tian J D，Niu H B 2005 Acta Phys.Sin.54 2034(in Chinese)[于斌，彭翔，田勁東，牛憨笨2005物理學報54 2034

[5]Deng X P，ZhaoD M 2011 Appl.Opt.50 6019

[6]Gerchberg R W，Saxton W O 1972 Optik 35 237

[7]Fienup J R 1982 Appl.Opt.21 2758

[8]Lu Y，LiQ，Dong Y H，GaoH D，Ma Z G 2001 J.Opt.Laser 12 365(in Chinese)[魯建業，李琦，董蘊華，高惠德，馬祖光2001光電子激光12 365]

[9]Fang L，Ye Y T，Wu Y F，Lu J J，Yang X M，Cheng Z Q Opt 2006 Opt.El.Eng.33 42(in Chinese)[方亮，葉玉堂，吳云峰，陸佳佳，楊先明，成志強2006光電工程33 42]

[10]LiS L，LiH T，Yang X J 2008 J.Appl.Opt.29 758(in Chinese)[李社蕾，李海濤，楊喜娟2008應用光學29 758]

[11]Yang G Z，Gu B Y 1981 Acta Phys.Sin.30 410(in Chinese)[楊國楨，顧本源1981物理學報30 410]

[12]Goodman J W 1996 Introduction toFourier Optics(2nd EdN.)(New York：McGraw-Hill)p55

[13]LiJ C 2009 Acta Opt.Sin.29 1163(in Chinese)[李俊昌2009光學學報29 1163]

[14]Biggs D S C，Andrews M 1997 Appl.Opt.36 1766

[15]Wen C L，JiJ R，Dou W H，Song Y S 2010 Acta Opt.Sin.30 2473(in Chinese)[溫昌禮，季家镕，竇文華，宋艷生2010光學學報30 2473]

[16]Deng X G，LiY P，Qiu Y，Fan D Y 1995 Chin.J.Laser B 4 447