李 淼 唐建波 劉 儉

（1.海軍駐431廠軍事代表室 葫蘆島 125004）（2.哈爾濱工業(yè)大學超精密光電儀器工程研究所 哈爾濱 150080）

1 引言

二元光學器件（Binary Optical Element，BOE）是利用計算機輔助設計而成的微浮雕結構衍射光學器件[1]，在光束的整形和變換[2]、光束耦合、聚焦、分束[3]、消除光學系統(tǒng)色差[4]及超分辨濾波器設計[5]等方面有廣泛的應用，二元光學元件計算機輔助設計已成為當今光學領域的一大研究熱點。通常的BOE設計算法大致可以分為兩類：基于線性變換方法如 GS算法[6]、楊－顧算法（Y－G算法）[7]等，基于搜索極值的優(yōu)化法如模擬退火算法（Simulated Annealing，SA算法）[8]、遺傳算法（GA 算法）[9]等，以及如參考文獻[1、3、10、13]等所述的眾多改進算法。BOE設計具有結構復雜、計算量大等特點，實現計算量小、收斂快、對初始值選取敏感性低、抗局部極值能力強、易于求得問題的全局最優(yōu)解是BOE計算機優(yōu)化設計算法研究的重要任務[1]。

傳統(tǒng)的GS算法收斂速度快，但其存在著對初始值敏感、易陷入局部極值等問題[10]；SA算法以一定的概率接受使評價函數變差的自變量取值，從而使算法能跳出局部極值，但其實質是以足夠高的初始溫度，緩慢的退火速度，大量的迭代次數及同一溫度下足夠多的擾動次數為保證的，從而導致了收斂效率十分低下，算法的計算量非常大[12]。本文在深入分析了模擬退火原理和傳統(tǒng)GS算法優(yōu)點的基礎上，結合模擬退火原理抗局部極值的思想和GS算法收斂目標明確、效率高的特點，提出了基于模擬退火原理的改進GS算法，并用改進算法和傳統(tǒng)GS算法對不同振幅分布目標進行了計算機仿真實驗。

2 BOE設計理論及基本設計算法

2.1 BOE設計理論

BOE設計理論基于圖1所示的傅里葉光學變換系統(tǒng)，輸入面P1和輸出面P2分別位于傅里葉透鏡的前焦面和后焦面，二元光學元件BOE置于P1面上，單色平面波fipt（x，y）垂直入射到BOE平面，經過BOE的衍射作用和透鏡L的傅里葉變換，在P2面得到所需的光場振幅分布Ftgt（ξ，η）。

BOE設計的實質是[11]：入射光波透過BOE后波前相位受到調制，然后經過透鏡的傅里葉變換作用在輸出平面上呈現出某種光強分布，并且BOE的相位分布與輸出平面上的光強分布是一一對應的，即一種相位分布對應一種輸出平面上的光強分布。這個過程在物理上是正過程，而BOE的設計正好是這個過程的逆過程，即從任意給定的輸出平面上的目標分布函數Ftgt（ξ，η）和入射光波分布函數fipt（x，y）去求解引起此光強分布的BOE的相位分布函數φ（x，y），使最終求得的相位分布函數φ（x，y）、入射光波分布函數fipt（x，y）及目標振幅分布函數Ftgt（ξ，η）滿足公式：

圖1 傅里葉光學變換系統(tǒng)

2.2 模擬退火原理

模擬退火算法（SA）作為局部搜索算法的擴展，其核心思想是在每一次修改模型的過程中，在先前最優(yōu)狀態(tài)的鄰域內隨機產生一個新的狀態(tài)模型，然后以一定的概率選擇鄰域中評價函數更好的狀態(tài)，這種接受新模型的方式使其成為一種全局最優(yōu)算法[3]。模擬退火原理基本思想描述如下[12]：

1）設初始狀態(tài)i為當前狀態(tài)，該狀態(tài)的能量為Ei，然后對該狀態(tài)做隨機微擾，得到一個新狀態(tài)j，新狀態(tài)的能量為Ej；

2）如果Ej＜Ei，則取狀態(tài)j為當前狀態(tài)；如果Ej＞Ei，則考慮到熱運動的影響，狀態(tài)j是否為重要狀態(tài)根據固體處于該狀態(tài)的概率來判斷，若狀態(tài)i和j的概率比值exp[（Ei－Ej）／（kb·T）]＞r，r為區(qū)間[0，1]上的隨機數，則以狀態(tài)j取代狀態(tài)i成為當前狀態(tài)，否則仍以狀態(tài)i為當前狀態(tài)，再重復以上新狀態(tài)的產生過程；

3）在大量固體狀態(tài)變換后，系統(tǒng)趨于能量較低的平衡狀態(tài)，固體狀態(tài)的概率分布趨于Boltzm－ann概率分布。

模擬退火原理的實質是在以100%的概率接受使評價函數變好的鄰域擾動自變量取值的同時以exp（（Ei－Ej）／T）＞rand（1）的概率接受使評價函數變壞的自變量改變，使其能跳出局部極值，模擬退火算法應用于二元光學計算機輔助設計時計算量很大，效率十分低下，這也是一直影響其真正走向實用的主要障礙[1]。

2.3 GS算法

GS算法應用于二元光學元件計算機輔助設計的基本思路是[6]：以初始相位和給定的入射光場分布，通過做正向傅里葉變換，得到輸出平面光場分布；在輸出平面用期望的光場振幅分布取代正向傅里葉變換得到的光場振幅分布，并且保持相位不變；然后做逆向傅里葉變換，得到輸入平面光場分布；在輸入平面用給定的光場振幅分布取代逆向傅里葉變換得到的振幅分布，并且保持相位不變；接著再次做正向傅里葉變換，如此循環(huán)下去，直到得到滿意結果或達到足夠多的循環(huán)次數為止。

GS算法是一種收斂速度很快的算法，但值得注意的是GS算法對初始參數很敏感，容易陷入局部極值[10]。

3 改進算法原理

GS算法是反復采用正傅里葉變換和逆傅里葉變換逐步收斂到目標函數的線性變換算法，具有目標明確、計算量小、收斂速度很快的特點，但用這種算法搜索函數的極小值時，如果函數在其定義域范圍內單調性不好，波形起伏頻繁，則收斂的效果就不會很好，不能跳出局部最優(yōu)解，收斂結果嚴重依賴于初始相位的選取[10]；模擬退火算法的核心思想是逐步的退火計劃和以一定的概率接受惡化解以跳出局部最優(yōu)解，但其計算量大，效率低下[1]。針對這兩種算法應用于二元光學元件設計的缺點，人們提出了多種改進算法，如參考文獻[1、3、10、13]所述。本文基于對模擬退火原理及GS算法特點的分析，設計了基于模擬退火原理的改進GS算法，并通過仿真實驗驗證了改進算法的優(yōu)越性，算法設計如下：

1）對初始溫度T0、目標溫度Tend、降溫速度dT、回火升溫閾值溫度Tlow、初始位相解φ0（x，y）等初始化。

2）對位相解φbest（x，y）做鄰域擾動，得到位相新解φnew（x，y）：

其中S（x，y）為鄰域擾動因子，當T＞Tlow時，采用指數遞減因子Se（x，y）：

其中k為迭代次數，a、b為自設系數；

當T≤Tlow時，采用依賴于溫度的似Cauchy遞減因子Sc（x，y）：

其中T為當前溫度。

3）以位相新解φnew（x，y）和給定的入射光場分布fipt（x，y）做正傅里葉變換得到輸出面上的復振幅分布F（ξ，η）：

由｜F（ξ，η）｜計算評價函數SSEnew和評價函數改變量ΔSSE：

4）如果評價函數改變量ΔSSE滿足條件：

（當T≤Tlow時取下降速度更快的T′，以保證最后只接受使評價函數變好的結果）

則以目標振幅分布和在步驟（3）中得到的位相解φ進行傅里葉逆變換得到f（x，y）：

并接受φangle（x，y）為φbest（x，y），SSEnew為SSEbest。

5）否則結束本次循環(huán)，φbest（x，y）、SSEbest保持不變，重復執(zhí)行2）～5），直到滿足當前溫度下內循環(huán)條件。

6）降溫，進行下一組迭代直到溫度下降到目標溫度Tend。

在步驟2）中當T≤Tlow時算法進行了回火升溫處理并且改變了鄰域擾動方式，依賴于溫度的似Cauchy分布產生新鄰域的優(yōu)點是：在高溫情況下搜索范圍大，在低溫時搜索僅在當前模型附近，即具有從高溫到低溫擾動范圍逐步平滑下降的特性，從而使其易于迅速跳出局部極值，并能很好保持搜索到的最優(yōu)解[12]，通過這一改進可以使模型具有再次跳出局部極值的機會，提高了模型的全局尋優(yōu)能力，并且逐步降低接受惡化解的可能性，使模型最終能很好保持搜索到的最優(yōu)解；在步驟4）中改進算法進行有選擇的傅里葉逆變換，即當滿足式（7）所示條件時進行傅里葉逆變換，減少了運算時間，提高了迭代效率。基于模擬退火原理的改進GS算法流程如圖2所示。

圖2 改進GS算法流程圖

圖3 目標振幅分布

通過對目標a所示圖形進行仿真實驗可以驗證算法對呈中心對稱分布圖形的設計性能，對目標b所示圖形進行仿真實驗則可驗證算法對任意分布圖形的設計性能。

均方誤差SSE和擬合系數η是衡量再現圖樣與目標圖樣差異的主要指標，能夠反映出再現圖樣的均勻性以及再現圖樣與目標圖樣的接近程度，在實際設計計算中，一般將均方誤差SSE或擬合系數η作為優(yōu)化設計算法的目標函數[10]，本文采用均方誤差SSE作為評價函數，擬合系數η作為輔助評價手段。

傳統(tǒng)的模擬退火算法采用對BOE元件各個相位單元進行逐個尋優(yōu)，這是導致SA算法效率低下的一個重要原因[13]，本文所探討的改進算法采用的是對整個BOE元件的相位面進行鄰域擾動尋優(yōu)，其核心思想是以GS迭代為基本框架，但在每一次GS迭代中都采用了模擬退火原理思想，在保證算法收斂效率的同時增加了算法的抗局部極值能力，能更好的全局尋優(yōu)，并且通過引入判斷條件減少了進行傅里葉逆變換的次數，收斂效率相對于GS算法得到了進一步提高。

4 仿真實驗與分析

本文采用圖1所示的傅里葉變換光路進行仿真實驗，設定入射高斯光束波長632.8nm，束腰半徑1.5mm，BOE元件邊長為x＝y(tǒng)＝5mm，輸入與輸出面采樣點數均為128×128。

為了驗證改進算法對二元光學元件計算機輔助設計的有效性，目標圖形采用如圖3（a）所示呈中心對稱分布的方形光斑，其中方形光斑邊長為2mm×2mm，以及如圖3（b）所示模擬任意形狀的北京奧運會會徽圖案：

φ（x，y）即為所求二元光學元件的位相分布。

通過取不同的初始相位可以判斷算法對初始相位的敏感性，因此當以圖3（a）為目標分布時，初始相位分別取以下四種情況：

1）φ0（x，y）取128×128維全0矩陣；

2）φ0（x，y）取128×128維全1矩陣；

3）φ0（x，y）取通過目標振幅分布Ftgt（ξ，η）傅里葉逆變換獲得的相位值；

4）φ0（x，y）取元素為[－π，π]之間隨機值的128×128維矩陣；

用GS算法和改進算法在以上四種情況下取相同的初始相位得到的仿真結果如表1所示，計算時間是在計算機CPU為AMD Athlon64，主頻2.31GHZ下得到的。

由表1可以看出在以上四種初始相位情況下，改進算法相對于GS算法均方誤差SSE平均減小了14.83%，擬合系數η平均提高了0.04%，迭代4000次所用時間平均減少了15.63%，證明改進算法相對于傳統(tǒng)GS算法具有計算量小、收斂快、全局尋優(yōu)能力強的特點；同時我們也發(fā)現當初始相位取前三種情況時，在每一種情況下GS算法迭代4000次都得到相同的結果，說明GS算法在很大程度上受制于初始相位的選取，抗局部極值能力弱；而改進算法因為采用了模擬退火原理的思想，通過鄰域擾動及以一定概率接受使評價函數變壞的結果，降低了對初始值選取的敏感性，提高了全局尋優(yōu)能力。

其中，Ftgt（ξ，η）為目標復振幅分布，F（ξ，η）為入射波光束垂直入射到BOE后經過透鏡的傅里葉變換在成像面上得到的復振幅分布：

表1 不同初始相位下傳統(tǒng)GS算法和改進GS算法仿真結果

圖4（a）所示為采用隨機初始相位時改進算法和GS算法的評價函數曲線圖，由圖可知GS算法在迭代2453次時達到局部極值，其后的迭代結果沒有明顯的改進，并且最終不能保持搜索到的最優(yōu)解；改進算法在第一階段收斂到局部極值后，經過再次升溫擾動，模型具有了跳出局部極值的能力，經過第二階段的擾動尋優(yōu)后，取得了更好的尋優(yōu)結果，并且在其后的迭代過程中逐步降低接受評價函數變壞的概率直到只接受使評價函數變好的擾動，使得改進算法的尋優(yōu)能力得到了明顯的提高。

圖4（b）所示為采用改進算法和GS算法取第四種初始相位情況下進行25000次迭代獲得的評價函數曲線圖，仿真結果顯示：GS算法獲得的最小均方誤差SSE＝4.5%，最優(yōu)擬合系數η＝97.77%，耗時755.2s；改進算法獲得的最小均方誤差SSE＝3.68%，最優(yōu)擬合系數η＝97.82%，耗時611.59s，并且自第二次擾動后評價函數值一直優(yōu)于GS算法并在迭代結束時獲得最優(yōu)解。進一步驗證了改進算法隨著迭代次數的增加可進一步搜索全局最優(yōu)解，具有很強的全局尋優(yōu)能力，而GS算法已陷入局部極值不能繼續(xù)尋優(yōu)。

圖4 GS算法和改進算法仿真實驗評價函數曲線圖

圖5（a）、（b）分別為在第4種情況下取相同的初始相位時GS算法和改進算法迭代4000次仿真得到的輸出面振幅分布截面圖，由圖可知改進算法獲得的結果中心方形光束部分能量更集中，并且方形光束之外區(qū)域光場分布更少，所得結果更接近目標振幅分布。

當以圖3（b）所示圖形（北京奧運會會徽圖案）為目標振幅分布時，采用改進算法及GS算法取相同的隨機初始相位分別迭代4000次得到的評價函數曲線如圖6所示。

圖5 輸出面歸一化振幅分布截面圖

圖6 GS算法和改進算法仿真實驗評價函數曲線圖

仿真結果顯示：GS算法獲得的最小均方誤差SSE＝6.577%，最優(yōu)擬合系數η＝98.04%，耗時131.3s；改進算法獲得的最小均方誤差SSE＝5.739%，最優(yōu)擬合系數η＝98.08%，耗時115.4s，并且自第二次擾動后評價函數值一直優(yōu)于GS算法并在迭代結束時獲得最優(yōu)解；證明改進算法對GS算法的優(yōu)越性同樣適用于對任意不對稱分布目標圖形的二元光學元件計算機輔助設計。

灰度直方圖反映了圖像中各個灰度級與該灰度級像素出現的頻率之間的關系，明亮圖像的灰度直方圖傾向于灰度級高的一側，對比度高的圖像直方圖覆蓋的灰度級很寬而且像素分布不均勻，集中分布在高灰度級和低灰度級上[14]，因此通過灰度直方圖可以評斷圖像的明亮程度及對比度。圖3（b）具有對比度高、目標圖像明亮的特點，因此若仿真圖形的直方圖覆蓋的灰度級寬、并且趨近于灰度級高的一側，則獲得的效果更好。圖7（a）、（b）分別為GS算法和改進算法仿真得到的輸出面圖形的灰度直方圖，由圖可知GS算法獲得的直方圖明暗波峰覆蓋范圍為182、明亮成分中心灰度級為196，而改進算法獲得的直方圖明暗波峰覆蓋范圍為194、明亮成分中心灰度級為208，說明改進算法獲得的圖形的直方圖覆蓋的灰度級寬、并且趨近于灰度級高的一側，因此改進算法的優(yōu)化設計效果更好。

圖7 輸出面圖像灰度直方圖

圖8（a）、（b）分別為GS算法和改進算法仿真得到的效果圖。

圖8 輸出面振幅分布

5 結語

改進算法以GS迭代為基本框架，在每一次迭代過程中又引入了模擬退火的思想，首先確保了收斂的速度，具有很高的效率，然后對每一次迭代的自變量都引入了鄰域擾動，并且以一定的概率接收惡化解，使算法具有跳出局部最優(yōu)解的能力，在溫度下降到一定閾值后進行升溫處理并且改變鄰域擾動因子，給模型再次跳出局部極值的機會，在很大程度上保證了全局尋優(yōu)。本文通過采用兩種不同的目標圖像進行仿真，驗證了改進算法的優(yōu)越性，其運行效率高，全局尋優(yōu)能力強，并且降低了對位相初值的敏感性，能對任意目標圖像進行計算機仿真設計，是進行二元光學元件計算機輔助設計的有效算法，在二元光學計算機優(yōu)化設計中具有很好的應用前景。

[1]康果果，謝景輝，莫曉麗，等.用改進的兩步模擬退火法進行二元光學元件的設計[J].光子學報，2008，37（7）：1416－1419.

[2]謝夢，曾曉東，安毓英.大功率半導體激光器二元光學消像散器件設計[J].紅外與激光工程，2007，36（1）：102－105.

[3]尹霄麗，余重秀.用改進的模擬退火算法設計二元光學陣列器件[J].光電子·激光，2001，12（2）：154－157.

[4]孫強，唐同斌，董國才，等.紅外場景產生器折射／衍射準直光學系統(tǒng)設計[J].紅外與激光工程，2007，36（6）：881－883.

[5]Jian Liu，Jiubin Tan，Yanchao Wang.Synthetic Complex Super－resolving Pupil Filter Based on Double－beam Phase Modulation[J].Applied Optics，2008，47（21）：3803－3807.

[6]Gerchberg RW，Saxton WO.A practical algorithm for the determination of phase from image and diffraction plane pictures[J].Optic，1972，35：237－246.

[7]楊國禎，顧本源，董碧珍.衍射光學元件的設計新方法[C]／／中國光學學會95年會特邀報告專輯，1994，39（5）：371－376.

[8]KIRKPATRICK S，GEL ATT C D，VECCHI M P.Optimization by simulated annealing[J].Science，1983，220（4598）：671－680.

[9]Gabreal Cormier，Roger Boudreau，Sylvain Theriault.Real－codid genetic algorithm for Bragg grating parameter synthesis[J].Opt.Soc.Am.B，2001，18（12）：1771－1776.

[10]鄒杰宇，盧亞雄，黃子強，等.基于改進GS算法的衍射光學光束整形元件的設計[J].紅外與激光工程，2006，35（S）：48－52.

[11]孔令彬，葉敦范，易新建，等.多相位256×256衍射微透鏡的設計及其光學性能研究[J].紅外與激光工程，2002，31（3）：253－256.

[12]陳華根，吳健生，王家林，等.模擬退火算法機理研究[J].同濟大學學報（自然版），2006，34（18）：1121－1125.

[13]謝敬輝，劉錫宇.基于DPDV算法的二元光學元件設計[J].光學技術，2000，26（3）：225－227.

[14]Rafael C.Gonzalez，Richard E.Woods.數字圖像處理[M].阮秋琦，阮宇智，譯.北京：電子工業(yè)出版社，2008，11.