菲涅爾衍射積分的角譜算法

韓佳成,郭文雅,趙亞楠,楊麗君,王穎,張素恒

(河北大學(xué) 物理科學(xué)與技術(shù)學(xué)院,河北 保定 071002)

菲涅爾衍射積分公式是標(biāo)量衍射公式的近似,可以計算傍軸區(qū)域的衍射光場分布[1-3],廣泛應(yīng)用于數(shù)字全息[4-6]、光學(xué)圖像加密[7-9]、光束整形[10-13]等領(lǐng)域.而在一般情況下,菲涅爾衍射積分并不存在解析解,需要借助于數(shù)值算法求解[14-15].角譜算法是一種常用的數(shù)值計算方法,它基于菲涅爾衍射的線性空間平移不變性,利用離散傅里葉變換將物光場分布與系統(tǒng)脈沖響應(yīng)函數(shù)的卷積運算轉(zhuǎn)換為物光角譜與系統(tǒng)傳遞函數(shù)的乘積運算,大大降低了計算的復(fù)雜度.得益于線性系統(tǒng)理論和快速離散傅里葉變換,角譜算法概念清晰,應(yīng)用方便,并且具有較高的計算效率,因而受到廣泛關(guān)注.

盡管已有文獻分別對數(shù)值計算菲涅爾衍射積分時二次相因子的采樣間隔[16-19]、物光矩陣的零填充[20-22]以及窗口尺寸計算[23-26]等因素進行了大量的討論.然而在角譜算法的具體實施過程中,這些因素相互依賴、彼此關(guān)聯(lián),它們共同決定了數(shù)值計算的正確性與計算精度,因此,需要從整體的角度統(tǒng)一梳理這些因素的相互關(guān)系,澄清有關(guān)問題,以確保角譜算法的正確實施.

本文首先從空間濾波的角度詳細(xì)闡明了角譜算法的具體實施過程,澄清了理想采樣條件的來源,指出理想采樣條件是角譜算法獲得正確結(jié)果的前提條件.按照理想采樣條件,采樣間隔反比于計算窗口尺寸,對采樣間隔的限制會影響到計算窗口尺寸的選擇.清晰解釋了離散傅里葉變換的循環(huán)卷積會引入卷繞誤差,指出選取較大的計算窗口尺寸相當(dāng)于進行零填充可以有效避免卷繞誤差,還明確給出了循環(huán)卷積與線性卷積相等的區(qū)域.詳細(xì)列出了影響計算窗口尺寸選取的因素,給出了計算窗口尺寸的選取依據(jù).這些方法和依據(jù)將有效地避免角譜算法應(yīng)用過程中所遇到的問題,從而保證計算的正確性與計算精度.然后,采用角譜算法數(shù)值計算了高斯光束照明余弦光柵的菲涅爾衍射場.數(shù)值解與解析解的對比表明,當(dāng)滿足理想采樣條件并正確選取計算窗口尺寸時,角譜算法可以達到非常高的計算精度.最后還采用角譜算法仿真了單透鏡相干成像實驗,進一步展示了角譜算法的應(yīng)用過程.

1 菲涅爾衍射的角譜算法

1.1 菲涅爾衍射角譜公式

在無限大不透明衍射屏上開有一長、寬分別為Wx、Wy的矩形衍射窗,物體分布于衍射窗內(nèi),以衍射窗中心為原點建立坐標(biāo)系,如圖1所示.當(dāng)單色相干光場照明衍射屏?xí)r,其后距離為z的觀察屏上的光場分布由菲涅爾衍射積分公式給出[25],見式(1).

圖1 菲涅爾衍射坐標(biāo)系示意Fig.1 Coordinate system of the Fresnel diffraction

(1)

其中,U0(x′,y′)和U(x,y)分別為衍射屏后表面物光場復(fù)振幅與觀察屏上衍射光場復(fù)振幅,λ為波長,k為波數(shù),對整個衍射屏后表面進行積分.exp(jkz)為光場縱向傳播引起的相位延遲,它在整個觀察屏上均勻分布,下文略去該相位因子.

若將菲涅爾衍射看成一個光學(xué)系統(tǒng),衍射屏后表面作為系統(tǒng)的輸入面,觀察屏作為系統(tǒng)的輸出面,并且該系統(tǒng)的輸入輸出關(guān)系由式(1)給出的卷積積分描述,于是菲涅爾衍射可看成一個線性空間平移不變系統(tǒng),其脈沖響應(yīng)函數(shù)為

(2)

這時菲涅爾衍射公式可表示為卷積形式,即

U(x,y)=U0(x,y)*h(x,y),

(3)

其中,*為卷積運算符.若從空間濾波的角度來理解菲涅爾衍射,衍射屏后表面的物光場經(jīng)濾波器濾波后即得到觀察屏上的衍射光場.

對于空間濾波系統(tǒng),還可以在頻譜域中討論輸入輸出關(guān)系.在標(biāo)量衍射理論中,光場分布的空間頻譜又稱為角譜.于是衍射屏后表面上物光場的角譜A0(ξ,η)與觀察屏上衍射光場的角譜A(ξ,η)分別為

(4)

A(ξ,η)=A0(ξ,η)H(ξ,η),

(5)

H(ξ,η)=exp[-jπλz(ξ2+η2)].

(6)

利用式(5)得到觀察屏上光場的角譜A(ξ,η)后,再進行傅里葉逆變換即可計算出觀察屏上的衍射場分布

(7)

1.2 菲涅爾衍射角譜公式的離散化

在數(shù)值計算菲涅爾衍射場時,需要對式(7)進行離散化.首先在衍射屏后表面以原點為中心選取一長、寬分別為Lx、Ly的計算窗口,使其包含衍射窗,如圖1所示.然后再以采樣間隔δx、δy沿x′、y′對計算窗口內(nèi)的物光場進行采樣,得到計算窗口內(nèi)物光的采樣矩陣

(8)

其中,各采樣點的坐標(biāo)分別為

(9)

圖2 利用ifftshift操作將坐標(biāo)原點從矩陣中心移至左上角Fig.2 Shifting the origin from the center of the matrix to the top left corner through the ifftshift operation

(10)

其中,坐標(biāo)取值分別為

(11)

對物光矩陣U0進行離散傅里葉變換得到物光角譜矩陣

A0=FFT2{U0},

(12)

其中,FFT2{·}代表二維快速離散傅里葉變換.物光角譜矩陣各取樣點的坐標(biāo)

(13)

由離散傅里葉變換的性質(zhì)可知,對于采樣間隔為δ、長度為M的離散序列,其離散頻譜序列的采樣間隔為1/Mδ.因此,物光頻譜矩陣A0在譜頻空間沿ξ與η方向的采樣間隔分別為

(14)

按照物光角譜矩陣所決定的采樣間隔(式(14))與采樣點坐標(biāo)(式(13))對衍射系統(tǒng)傳遞函數(shù)H(ξ,η)進行采樣,得到傳遞函數(shù)的采樣矩陣

(15)

然后將物光角譜矩陣A0與傳遞函數(shù)矩陣H代入式(7),可計算出觀察平面上的衍射場分布矩陣

Uw=IFFT2{A0°H},

(16)

其中,IFFT2{·}代表二維快速離散傅里葉逆變換,°代表Hadamard乘積.但直接逆變換得到的衍射場分布矩陣Uw的坐標(biāo)原點位于矩陣左上角,需要通過fftshift操作,將坐標(biāo)原點移至矩陣的中心,即

(17)

調(diào)整后的坐標(biāo)分別為

(18)

這種數(shù)值計算菲涅爾衍射積分的方法稱為角譜算法.

1.3 理想采樣條件

利用角譜算法數(shù)值計算菲涅爾衍射積分時,間接涉及對傳遞函數(shù)進行采樣.從式(6)可以看出,傳遞函數(shù)由單位振幅的二次相因子構(gòu)成,不是帶限函數(shù).在整個頻譜平面上,采樣總無法滿足采樣定理的要求.但在數(shù)值計算時,僅對有限頻率范圍內(nèi)的傳遞函數(shù)進行采樣.由式(13)、(14)可知,沿ξ與η兩個方向?qū)鬟f函數(shù)采樣的頻率范圍分別為

(19)

其中,M=Lx/δx,N=Ly/δy.2009年VOELZ等[18]采用局域空間頻率的概念證明只有當(dāng)空域采樣間隔滿足

δx=λz/Lx,δy=λz/Ly,

(20)

傳遞函數(shù)H(ξ,η)的采樣矩陣H(u,v)才與計算窗口內(nèi)脈沖響應(yīng)函數(shù)h(x,y)的采樣矩陣h(m,n)互為離散傅里葉變換對,并稱式(20)為理想采樣條件.值得一提的是,在理想采樣條件下,角譜算法與脈沖響應(yīng)函數(shù)法[18]完全等價.

理想采樣條件是確保角譜算法獲得正確結(jié)果的前提條件,下述討論均以理想采樣條件為基礎(chǔ).考慮到理想采樣條件,傳遞函數(shù)的采樣矩陣可表示為

H(u,v)=exp[-jπ(u2/M+v2/N)],

(21)

其中,各取樣點的坐標(biāo)分別為

(22)

1.4 循環(huán)卷積與卷繞誤差

利用角譜算法數(shù)值計算菲涅爾衍射積分時,會受到離散傅里葉變換帶來的循環(huán)卷積的影響.按照離散傅里葉變換的循環(huán)卷積定理,式(16)可以表示為

Uw=U0○*h,

(23)

其中,○*代表循環(huán)卷積.h=IFFT2{H}為計算窗口內(nèi)脈沖響應(yīng)函數(shù)的采樣矩陣,其中各采樣點坐標(biāo)取值由式(11)給出.由式(23)可知,計算得到的衍射光場分布矩陣Uw為計算窗口內(nèi)的物光采樣矩陣U0與脈沖響應(yīng)函數(shù)采樣矩陣h的循環(huán)卷積.然而,實際觀察屏上衍射光場分布U應(yīng)為U0與h的線性卷積,即

U=U0*h.

(24)

循環(huán)卷積與線性卷積的差常稱為卷繞誤差.

圖3 計算窗與衍射窗尺寸相同時周期延拓后的U0和hFig.3 Period extended U0and h when the sizes of the calculation window and diffraction window are the same

圖4 計算窗大于衍射窗尺寸時周期延拓后的U0和hFig.4 Period extended U0and h when the calculation window size is larger than the diffraction window size

(25)

(26)

1.5 計算窗口尺寸的選擇

在利用角譜算法計算菲涅爾衍射時,有多種因素影響計算窗口尺寸的選取.計算窗口尺寸的選取主要考慮到以下因素:

1)計算窗口尺寸應(yīng)滿足理想采樣條件(式(20)),在衍射距離確定的情況下,計算窗口尺寸與采樣間隔并不獨立.對采樣間隔的限制會影響到計算窗口尺寸的選取.對衍射窗內(nèi)物光場采樣應(yīng)滿足采樣定理,即采樣間隔應(yīng)滿足

1/δx≥Bx, 1/δy≥By,

(27)

其中,Bx、By分別為物光場沿x′、y′的帶寬.考慮到理想采樣條件,要求計算窗口尺寸滿足

Lx≥λzBx,Ly≥λzBy.

(28)

2)若觀察窗的尺寸為Gx、Gy,以z軸與觀察面的交點為中心,如圖1所示.經(jīng)采樣后,觀察窗內(nèi)采樣點坐標(biāo)集合為

(29)

Lx≥Wx+Gx,Ly≥Wy+Gy.

(30)

若要得到完整的衍射場,需要觀察窗的尺寸最小[22],

Gx=Wx+λzBx,Gy=Wy+λzBy.

(31)

這時計算窗口的尺寸應(yīng)滿足

Lx≥2Wx+λzBx,Ly≥2Wy+λzBy.

(32)

3)計算窗口的尺寸越大,采樣點越多.在理想采樣條件下,計算窗口內(nèi)沿x′、y′的采樣點數(shù)分別為

(33)

為了減小計算量,Lx和Ly應(yīng)取在式(30)的下限附近,并使式(33)得到的M、N恰好為整數(shù),這也方便計算.

2 數(shù)值仿真實驗

2.1 高斯光束照明余弦光柵的菲涅爾衍射

高斯光束照明無限大二維余弦光柵時的菲涅爾衍射積分具有解析解.這時可以將數(shù)值計算結(jié)果與采樣后的解析解進行比較,來檢驗數(shù)值算法的計算精度.

在圖1中的衍射窗內(nèi)放置一個二維余弦光柵,沿z軸正向傳播的基模高斯光束從衍射屏左側(cè)垂直入射到該余弦光柵上,且高斯光束的束腰正好位于光柵所在平面,則衍射屏后表面上的光場分布可表示為

f0(x′,y′)=e-π(x′2+y′2)/α2cos(2πξ0x′)cos(2πη0y′),

(34)

f(x,y)=[fξ0(x)+f-ξ0(x)][fη0(y)+f-η0(y)]/4,

(35)

其中,fξ(x)的具體形式為

fξ(x)=Kξexp[-π(x-λzξ)2/σ2]exp[j2π(κqx2+κlξx)],

(36)

其中,

由于照明基模高斯光束的能量幾乎全部位于-2α≤x′≤2α,-2α≤y′≤2α,因此可以認(rèn)為衍射窗尺寸為Wx=Wy=W=4α?xí)r,觀察屏上的衍射場分布基本不變,仍可由式(35)來描述.

(37)

a.數(shù)值解U的振幅分布;b.解析解f的振幅分布;c.數(shù)值解與解析解差值E的模值分布圖5 高斯光束照明余弦光柵的菲涅爾衍射圖樣Fig.5 Fresnel diffraction pattern of Gaussian beam illuminating cosine grating

為了定量描述數(shù)值計算誤差,選用均方信噪比RSN作為評價標(biāo)準(zhǔn),其定義為

RSN=10 lg(∑mn|f(m,n)|2/∑mn|E(m,n)|2),

(38)

其中求和對于矩陣中所有元素進行.根據(jù)該定義,數(shù)值計算結(jié)果U相對于解析解f的均方信噪比為RSN=116.7 dB.當(dāng)采用理想采樣間隔和適當(dāng)?shù)挠嬎愦翱跁r,角譜算法的性能十分優(yōu)異.該仿真實驗展示了解析物體經(jīng)歷單次菲涅爾衍射的計算過程.

2.2 單透鏡相干成像

通常采用菲涅爾衍射分析成像過程.下面采用角譜算法數(shù)值仿真理想薄透鏡的相干成像過程.如圖6所示,理想薄透鏡的焦距F=40 mm,其光瞳為邊長d=5 mm的正方形,物面與像面分別位于透鏡前后2倍焦距處.在物面上,物體放置于物面上邊長W=4 mm的正方形衍射窗內(nèi),衍射窗的中心與光軸重合.波長為λ=500 nm的單色平面波正入射照明物體,透過物體的光波先后經(jīng)過2次菲涅爾衍射后到達像面.

圖6 單透鏡相干成像示意Fig.6 Schematic diagram of single-lens coherent imaging

仿真物體為USAF-1951分辨率測試板(局部)以δ=4 μm的采樣間隔采樣后得到的P×Q=1 000×1 000像素的圖像,如圖7a所示.按照理想采樣條件,當(dāng)采樣間隔δ=4 μm時,計算窗口需要取邊長L=2λF/δ=10 mm的正方形,其中采樣點數(shù)為M×N=2 500×2 500.因此,首先需要將圖7a所示圖像四周零填充至2 500×2 500像素,作為物光采樣矩陣U0.利用角譜算法可以計算出透鏡前表面光場采樣矩陣Ul.由于計算窗的邊長L=10 mm,衍射窗的邊長W=4 mm,在計算出的透鏡前表面的衍射場分布矩陣Ul中未受卷繞影響的區(qū)域為L-W=6 mm的正方形.而透鏡的光瞳為邊長d=5 mm的正方形,因此,在透鏡光瞳內(nèi)是正確的衍射場.

a.離散物體圖像;b.等大成像仿真結(jié)果的振幅分布;c.輕微離焦成像仿真結(jié)果的振幅分布圖7 單透鏡相干成像中離散物體圖像與數(shù)值仿真結(jié)果Fig.7 Discrete object image and numerical simulation results in single-lens coherent imaging

透鏡的復(fù)振幅透過率為

T(x,y)=exp(-jk(x2+y2)/2F)rect(x/d)rect(y/d),

(39)

其中,rect代表矩形函數(shù).對透鏡光瞳內(nèi)的透過率函數(shù)以采樣間隔δ=4 μm進行采樣,采樣矩陣記為

(40)

其中,λF/δ2=1 250,各取樣點的坐標(biāo)為

(41)

其中,R=d/δ=1 250.

考慮到像面位于透鏡后2F處,計算窗口尺寸仍需取為L=10 mm.因此,需要將透鏡的采樣矩陣四周零填充至2 500×2 500像素,再與透鏡前表面的衍射場分布矩陣進行Hadamard乘積,得到透鏡后表面的光場分布矩陣U′l.再次利用角譜算法可以計算出像面上光場采樣矩陣Ui.對像面光場采樣矩陣Ui進行裁剪,僅保留中心1 000×1 000的矩陣元素,即為物體的像,如圖7b所示.為了便于比較,圖7b展示的是旋轉(zhuǎn)180°之后的像.

若將像面向后移動2 mm,即像面到透鏡的距離z=82 mm.這時計算窗口尺寸需要取為L=λz/δ=10.25 mm,其中采樣點數(shù)為M×N=2 563×2 563.這就需要將透鏡前表面的衍射場分布矩陣Ul和透鏡透過率函數(shù)的采樣矩陣T均四周填充至2 563×2 563,再進行Hadamard乘積得到透鏡后表面的光場分布矩陣U′l.然后利用角譜算法,計算出像面光場分布采樣矩陣Ui.最后對像面光場采樣矩陣Ui進行裁剪,僅保留中心1 000×1 000的矩陣元素,即為物體輕微離焦的像,如圖7c所示.同樣,圖7c展示的也是旋轉(zhuǎn)180°之后的像.

由圖7可以看出,利用角譜算法可以準(zhǔn)確地數(shù)值仿真不同像距下的單透鏡相干成像實驗.仿真實驗展示了離散物體經(jīng)歷級聯(lián)菲涅爾衍射的計算過程.單透鏡相干成像仿真實驗的源程序參見開源倉庫[28],以便參考.

這里通過衍射追跡的方法仿真了物體在相干照明下經(jīng)過有限孔徑的理想薄透鏡的成像情況.然而在非相干照明下,無法通過菲涅衍射直接計算物光場在成像系統(tǒng)中的傳播情況.但這時可以通過菲涅爾衍射計算出成像系統(tǒng)的相干傳遞函數(shù),進而得非相干成像系統(tǒng)的光學(xué)傳遞函數(shù),然后再利用光學(xué)傳遞函數(shù)得到像的強度分布[29].

3 結(jié)論

詳細(xì)討論了菲涅爾衍射積分的角譜算法,給出了算法的實施過程,澄清了理想采樣條件、卷繞誤差產(chǎn)生的原因和影響范圍以及計算窗口尺寸的選取依據(jù)等相關(guān)問題.高斯光束照明余弦光柵的衍射仿真實驗和單透鏡相干成像仿真實驗表明,按照本文給出的依據(jù)和過程執(zhí)行計算,可以有效地避免應(yīng)用角譜算法時出現(xiàn)錯誤,保證數(shù)值計算精度.

角譜算法和單次傅里葉變換算法[26]是2種常用的基于離散傅里葉變換的計算菲涅爾衍射積分的數(shù)值算法.在相同條件下,2種算法的計算精度基本相同,但是由于角譜算法需要進行2次離散傅里葉變換和1次矩陣Hadamard乘積運算,而單次傅里葉變換算法需要進行1次離散傅里葉變換和2次Hadamard乘積運算.因此,角譜算法花費的時間較長.以2.1節(jié)中高斯光束照明余弦光柵衍射場計算為例,角譜算法花費的時間約為單次傅里葉變換算法的1.47倍.

角譜算法基于空間濾波模型,概念比較清晰,直觀地給出了卷繞誤差產(chǎn)生的原因和影響范圍.在單次傅里葉變換算法中,卷繞誤差同樣存在.本文給出的避免卷繞誤差的方法同樣適用于單次傅里葉變換算法,為從事衍射計算相關(guān)領(lǐng)域的研究人員提供了一個簡明參考.