隨機粗糙面散射中遮蔽效應(yīng)算法的改進

宋 哲，徐 哲，龐 勇，馬思遠(yuǎn)，劉春楠

(遼寧師范大學(xué) 物理與電子技術(shù)學(xué)院，遼寧 大連 116029)

自然界中絕大部分物體表面都是隨機粗糙面，研究隨機粗糙面的電磁散射特性是獲取物體表面特征的關(guān)鍵途徑，在遙感領(lǐng)域[1]、零件檢測領(lǐng)域[2]、醫(yī)學(xué)領(lǐng)域[3]等領(lǐng)域有廣泛應(yīng)用.隨機粗糙面散射特性的研究主要有兩種方式：數(shù)值計算和理論近似.數(shù)值計算方法主要有矩量法[4]、有限元法[5]、時域有限差分法[6]等，特點是計算量大，結(jié)果較為精確.理論近似方法主要有微擾法[7]、基爾霍夫近似法[8]等，特點是計算量小，易于理解.其中，基爾霍夫近似法計算效率高，且精度與數(shù)值計算法接近，得到廣泛應(yīng)用.其核心思想是將粗糙面劃分為有限個微面元，每個微面元近似為平面，則光波在微面元上發(fā)生菲涅耳反射，由此獲得微面元上的本地場，再結(jié)合Statton-Chu積分方程得到粗糙面在空間中的散射場.

基爾霍夫近似方法的局限性在于它認(rèn)為粗糙面上所有的面元都會被照射以及被散射觀測點觀測到，但事實上，當(dāng)光波以一定角度入射到粗糙面上時會產(chǎn)生陰影區(qū)，使部分面元沒有光照射，同時，對于空間某一觀測角度也會存在部分面元的反射光因為被其他面元遮擋而無法到達觀測點，這種現(xiàn)象就是遮蔽效應(yīng).Beckmann首先提出了遮蔽效應(yīng)的概念，并給出了單站散射的遮蔽函數(shù)表達式[9].Smith研究并給出了高斯分布下隨機粗糙表面的遮蔽概率函數(shù)[10].Wagner在Beckmann單站散射遮蔽函數(shù)的基礎(chǔ)上提出了雙站散射遮蔽函數(shù)[11].王安祥等人將實驗測量和遺傳模擬退火算法相結(jié)合，對雙向反射分布函數(shù)(BRDF)模型中遮蔽函數(shù)的參數(shù)進行了反演，發(fā)現(xiàn)反演得到的遮蔽函數(shù)參數(shù)能夠反映表面粗糙度[12].柳祎等人建立了含有遮蔽函數(shù)的偏振度模型，研究了不同粗糙度下金屬與非金屬紅外自發(fā)輻射的偏振度，發(fā)現(xiàn)粗糙度對非金屬表面紅外自發(fā)輻射偏振度的影響大于金屬，可用于區(qū)分金屬和非金屬[13].郭立新等人利用基爾霍夫近似給出了考慮遮蔽效應(yīng)的分形粗糙表面散射截面近似公式，并討論了不同分形維數(shù)和空間基頻下遮蔽效應(yīng)對雙站散射截面的影響[14].馬?？频热诉\用射線追蹤法研究了一維粗糙面的遮蔽效應(yīng)和二次散射，計算了一維粗糙面的散射系數(shù)[15].陳明等人利用射線追蹤法研究了帶限分形粗糙面的遮蔽效應(yīng)，分析了影響粗糙面遮蔽效應(yīng)的因素[16].Yoshifumi等人提出一種包含多次散射的粗糙面射線追蹤(RSRT)模型，精確計算了散射分布[17].

當(dāng)前遮蔽效應(yīng)的研究主要有兩種方式：遮蔽概率函數(shù)和幾何遮蔽.遮蔽概率函數(shù)法是通過提出表面任意面元被遮蔽概率的函數(shù)來計算遮蔽效應(yīng)的，計算量小但精確度稍差.幾何遮蔽法是通過射線追蹤，根據(jù)光線與粗糙面的位置關(guān)系來逐點判斷被遮蔽情況，計算準(zhǔn)確但計算量大.本文將傳統(tǒng)的幾何遮蔽方法進行了改進，在不改變原有精確度的情況下，縮短了幾何遮蔽效應(yīng)的運算時間，并在此基礎(chǔ)上采用基爾霍夫近似法分別計算了S波和P波入射到隨機粗糙面時的散射場分布.

1 基于Kirchhoff近似法分析隨機粗糙面的散射

呈現(xiàn)高斯分布的隨機粗糙面在自然界中十分常見，可用線性濾波法來構(gòu)建隨機粗糙面.描述粗糙面的指標(biāo)主要有相關(guān)長度l(二維粗糙面x和y方向的相關(guān)長度分別為lx，ly)和均方根高度δ，則隨機粗糙面各處的高度分布可表示為[18]

(1)

其中，M、N為x和y方向離散點的數(shù)目，Lx、Ly分別為粗糙面在x和y方向的長度，(xm,xn)為離散點坐標(biāo)，當(dāng)離散點的間隔為Δx和Δy時，xm=mΔx，yn=nΔy，km、kn為x和y方向的波數(shù)，km=2πm/Lx，kn=2πn/Ly，F(xiàn)(km,kn)為傅里葉系數(shù)，可以表示為

(2)

N(0,1)為滿足高斯分布的一組隨機生成數(shù)，S(km,kn)是二維高斯粗糙面的功率譜密度，表達式為

(3)

根據(jù)式(1)～式(3)，可分別生成M=N=500、lx=ly=3 μm、δ=0.1 μm和δ=0.3 μm的隨機粗糙面，如圖1所示.

圖1 隨機粗糙面示意圖

當(dāng)平面波ki照射到隨機粗糙面上時會發(fā)生散射，根據(jù)基爾霍夫近似理論，將粗糙面劃分為若干個近似平面的微面元，光波在微面元上發(fā)生菲涅爾反射.圖2是光波在微面元上反射時入射場與反射場的示意圖，θ為微面元上的光波入射角，θr=θ為微面元上的光波反射角，n為微面元的法向量，kr為反射波矢，Ei為入射場，Er為反射場，則微面元上的本地場E(r′)為入射場與反射場的疊加，即

圖2 微面元入射場和反射場示意圖

E(r′)=Ei(r′)+Er(r′).

(4)

r′為微面元的位置矢量.借助于Stratton-Chu積分方程，可得到半球空間內(nèi)任意位置的散射場[19]:

(5)

其中,r為空間某一觀察點的位置矢量，ds′為單個微面元面積，G(r,r′)為格林函數(shù)，其遠(yuǎn)場近似形式為

(6)

其中，r為空間觀測點到原點的距離，ks為散射波矢，k為波數(shù).

2 幾何遮蔽效應(yīng)

利用式(4)得到的散射場包含了粗糙面所有面元的貢獻，忽略了有些面元可能沒有被照射到或無法被觀測到的情況，即面元的遮蔽效應(yīng)，這些面元對散射場是沒有貢獻的，在入射角或散射角較大時遮蔽效應(yīng)更加顯著，因此應(yīng)該對粗糙面散射場進行遮蔽修正.

遮蔽效應(yīng)包括入射遮蔽和散射遮蔽，入射遮蔽指由于粗糙面上的起伏使入射光被遮擋而無法照到的面元，散射遮蔽指由于面元反射的光線被其它面元遮擋而不能直接到達半球空間觀測點的面元，二者需要分開討論.

2.1 入射遮蔽效應(yīng)

傳統(tǒng)的入射幾何遮蔽主要是基于射線追蹤法來計算被遮蔽的面元，即將光線認(rèn)為是三維坐標(biāo)系下的射線，通過追蹤射線與微面元是否相交來判斷哪些面元是被遮蔽的.如判斷微面元1是否被遮蔽的思路是：假設(shè)入射光線ki能照到面元1上，計算粗糙面上除面元1外所有面元是否與ki有交點，若有交點，則面元1被遮蔽，否則面元1沒有被遮蔽，如圖3中面元1被面元2遮蔽.這種方式雖然能夠很準(zhǔn)確的判斷粗糙表面各面元被遮擋的情況，但時間復(fù)雜度較高.

圖3 面元1被遮擋情況示意圖

設(shè)入射光線ki的方位角為0°，建立坐標(biāo)系使xoz面為入射面，ki的方向指向x軸的正方向，z軸的負(fù)方向，因此二維隨機粗糙面的入射遮蔽效應(yīng)可以轉(zhuǎn)化為一維粗糙面的遮蔽效應(yīng)，如圖4所示.從圖中可以看出，入射光線在粗糙面上會產(chǎn)生陰影區(qū)，使陰影中的面元被遮蔽，每段遮蔽面元都起始于一個被照射面元和被遮蔽面元的“臨界面元”，如面元1，終止于入射到面元1的光線延長線與粗糙面的交點面元2，這段被遮蔽的面元都位于入射光線的下方，因此判斷一個面元是否被遮蔽，只要找出該面元在入射光方向最近的臨界面元，再判斷該面元是否在入射光線下方即可.

圖4 入射遮蔽效應(yīng)示意圖

首先沿x正方向計算入射光線ki與每個微面元法線n是否滿足

ki·n=0.

(7)

若式(7)成立，那么該面元即為臨界面元.設(shè)臨界面元的坐標(biāo)為(x1,z1)，則通過臨界面元的入射光線方程為

(8)

其中,kix和kiz為ki在x和z方向的分量，x和z為光線上的點.

設(shè)兩個相鄰臨界面元之間任一面元坐標(biāo)為(x2,z2)，從臨界面元起沿x方向依次代入式(8)，若

(9)

2.2 散射遮蔽效應(yīng)

散射光線的方向可以為半球空間中的任意方向，因此應(yīng)考慮二維粗糙面的遮蔽.傳統(tǒng)二維粗糙面光線追蹤法的思路是：判斷一個面元是否被遮蔽須遍歷除該面元及其相鄰面元以外的所有面元，求該面元的散射光線與其他面元是否有交點，若有交點，則該面元被遮蔽；否則沒有被遮蔽.

在實際的光線傳播過程中，散射光能夠經(jīng)過的面元是有限的，如圖5(a)所示，ks是二維粗糙面上任意面元1的散射光，虛線代表散射光可能經(jīng)過的面元.圖5(b)是二維粗糙面在xoy面的投影，其中,網(wǎng)格為粗糙面上各微面元的投影，k′s為ks的投影，陰影方框為k′s經(jīng)過的投影面元.因此在計算其他面元是否對面元1的散射光線有遮擋只需判斷陰影投影面元對應(yīng)的粗糙面面元對其是否有遮擋，若有遮擋，則面元1 被遮蔽，無需判斷粗糙面上除面元1以外的所有面元，可以大大減少計算量.

圖5 粗糙面散射遮蔽效應(yīng)示意圖

首先，若散射光線ks與面元1的法向量n滿足ks·n<0，說明面元1的散射光不能直接發(fā)射到半球空間，面元1被遮擋.若ks·n≥0，則設(shè)k′s的斜率kxy為

(10)

圖6 搜索k′s經(jīng)過的面元

設(shè)面元1的坐標(biāo)為(x1,y1,z1)，法向量為n=(nx,ny,nz)，面元1的散射光線為ks=(kx,ky,kz)，與k′s相交的投影面元對應(yīng)的粗糙面面元2坐標(biāo)為(x2,y2,z2)，法向量為n′=(n′x,n′y,n′z)，則面元1的散射光線方程為

(11)

面元2所在平面的平面方程為

(x-x2)n′x+(y-y2)n′y+(z-z2)n′z=0.

(12)

將式(11)和式(12)聯(lián)立可得到面元1的散射光線ks與面元2所在平面的交點c=(xc,yc,zc).若xc和yc滿足

(13)

則交點在面元(x2,y2,z2)內(nèi)，說明面元1被遮蔽.

2.3 改進效果

選取lx=ly=4 μm、δ=0.5 μm的隨機粗糙面，計算入射角為θi=30°，散射角為θs=70°，散射方位角為φs=45°時的遮蔽效應(yīng)，改進前與改進后的計算時間統(tǒng)計如表1所示，其中,N為x和y方向劃分面元數(shù)量，tx為原始算法的程序執(zhí)行時間，ty為改進后算法的程序執(zhí)行時間，運算平臺為i7-7500U CPU，頻率為2.90 GHz，內(nèi)存為8 GB，編程軟件為matlab R2019b.

由表1可知，相對傳統(tǒng)射線追蹤算法，本文改進算法計算遮蔽效應(yīng)的時間大大縮短，并且隨著粗糙面面元數(shù)量的增加，改進算法的優(yōu)勢更明顯.

表1 改進后的程序執(zhí)行時間

3 隨機粗糙面散射光強數(shù)值計算

選取lx=ly=4 μm、δ=0.4 μm的金屬鐵隨機粗糙表面為對象，當(dāng)入射光波長為λ=0.905 μm時，金屬鐵的折射率為n=3.12+3.87i.根據(jù)式(4)和上述的遮蔽效應(yīng)，分別模擬計算了S波、P波以30°角入射時，入射面內(nèi)散射光強的分布，如圖7所示，圖中散射光強均歸一化到S波入射時散射光強的最大值.由圖可知本文模擬結(jié)果與文獻[20]的結(jié)果符合較好，說明本文提出的遮蔽效應(yīng)計算方法是正確的.圖8給出了S波、P波30°角入射時半球空間散射光強的分布，可見半球空間中散射光強的分布是關(guān)于入射面(散射方位角為0°)對稱的，并主要集中在與入射光鏡像方向(即散射角為30°，散射方位角為0°)附近，散射角方向半寬度約40°，散射方位角方向半寬度約60°，近似在鏡像方向有最大值，P波散射光強小于S波散射光強.

圖7 入射面內(nèi)散射光強分布

圖8 全角度散射光強分布

4 結(jié) 論

粗糙面在自然界中十分常見，對粗糙面的散射場進行深入研究在目標(biāo)檢測、識別等領(lǐng)域具有重要意義.本文利用線性濾波法生成了不同相關(guān)長度和均方根高度的高斯型二維隨機粗糙面，采用基爾霍夫近似方法分析了隨機粗糙面的散射.針對傳統(tǒng)幾何遮蔽運算中時間復(fù)雜度高的問題，提出了幾何遮蔽的改進算法.對于入射遮蔽，判斷某一面元是否被遮蔽只要找出該面元在入射光方向最近的臨界面元，再判斷該面元是否在過臨界面元的入射光線下方即可.對于散射遮蔽，判斷面元是否被遮擋，只需要判斷該面元的散射光線與此散射光在粗糙面上投影所經(jīng)過的面元是否有交點即可，無需判斷粗糙面上所有面元.入射遮蔽的時間復(fù)雜度由O(N4)(O為描述時間復(fù)雜度的函數(shù))降為O(N2)，散射遮蔽的時間復(fù)雜度由O(N4)降為O(N3)，計算效率得到很大提升.模擬計算了S波和P波入射到鐵表面時入射面內(nèi)散射光強的分布，與文獻結(jié)果進行對比，基本吻合，證明了改進后方法的有效性.計算了隨機粗糙面的全角度散射光強，結(jié)果表明散射光強分布關(guān)于入射面對稱，且主要集中在與入射光鏡像方向附近，散射角方向半寬度約40°，散射方位角方向半寬度約60°，P波散射光強小于S波散射光強.研究結(jié)果能夠為目標(biāo)探測提供理論參考.