劉宏展 紀(jì)越峰

1)(北京郵電大學(xué)，信息光子學(xué)與光通信國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100876)

2)(華南師范大學(xué)，廣東省微納光子功能材料與器件重點(diǎn)實(shí)驗(yàn)室，廣州 510006)

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

1 引言

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

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

2 算法描述與分析

要進(jìn)行相位恢復(fù)迭代運(yùn)算，就必須進(jìn)行光場(chǎng)衍射及逆衍射運(yùn)算，因此我們將首先簡(jiǎn)單介紹文中要用到的角譜理論，然后重點(diǎn)對(duì)加成算法進(jìn)行描述與分析.

2.1 角譜理論

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

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

2.2 幅度梯度加成迭代算法

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

其中W 為輸出面整形光場(chǎng)的范圍，βk，gk，γk的表達(dá)式是參考文獻(xiàn)[14]得到的.

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

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

3 加成算法仿真與結(jié)果分析

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

3.1 參數(shù)設(shè)置

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

3.2 仿真與分析

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

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

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

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

圖2 不同隨機(jī)初始相位下迭代次數(shù)的分布

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

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

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

圖3中，有兩條圖線，橫坐標(biāo)與圖2的含義一致，縱坐標(biāo)表示每次加成算法能達(dá)到的最小誤差e.上邊那條對(duì)應(yīng)η=90%，10次運(yùn)算的結(jié)果是，其迭代誤差e密集分布在0.14兩側(cè)，即迭代誤差趨向于0.14，可見，加成算法的迭代收斂具有一致性.同時(shí)，與初始相位為0時(shí)e=0.1520(見表1)相比，其誤差有所減小.下邊的圖線對(duì)應(yīng)η=94%，其迭代誤差e更加密集地分布于0.09附近，與η=90%相比，此時(shí)迭代誤差更加減小，且收斂一致性更優(yōu).這也從一個(gè)側(cè)面說明，加成算法中用能量集中度η或者迭代誤差e作為判斷依據(jù)是統(tǒng)一的，這是因?yàn)椋芰扛辛?，自然迭代所能達(dá)到的最小誤差便變小;而迭代誤差變小了，則光束能量就更集中.

3.3 仿真實(shí)例

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

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

4 結(jié)論

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

[1]Huang L X，YaoX，CaiD M，GuoY K，YaoJ，GaoF H 2010 Chin.J.Laser.37 1218(in Chinese)[黃利新，姚新，蔡冬梅，郭永康，姚軍，高福華2010中國(guó)激光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)[俞建杰，譚立英，馬晶，韓琦琦，楊玉強(qiáng)，李密2009中國(guó)激光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)[鄔融，華能，張曉波，曹國(guó)威，趙東峰，周申蕾2012物理學(xué)報(bào)61 224202]

[4]Yu B，Peng X，Tian J D，Niu H B 2005 Acta Phys.Sin.54 2034(in Chinese)[于斌，彭翔，田勁東，牛憨笨2005物理學(xué)報(bào)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)[魯建業(yè)，李琦，董蘊(yùn)華，高惠德，馬祖光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)[方亮，葉玉堂，吳云峰，陸佳佳，楊先明，成志強(qiáng)2006光電工程33 42]

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

[11]Yang G Z，Gu B Y 1981 Acta Phys.Sin.30 410(in Chinese)[楊國(guó)楨，顧本源1981物理學(xué)報(bào)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光學(xué)學(xué)報(bào)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光學(xué)學(xué)報(bào)30 2473]

[16]Deng X G，LiY P，Qiu Y，F(xiàn)an D Y 1995 Chin.J.Laser B 4 447