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

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

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

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

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

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

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

(1)

其中，M、N為x和y方向離散點(diǎn)的數(shù)目，Lx、Ly分別為粗糙面在x和y方向的長度，(xm,xn)為離散點(diǎn)坐標(biāo)，當(dāng)離散點(diǎn)的間隔為Δx和Δy時(shí)，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)為滿足高斯分布的一組隨機(jī)生成數(shù)，S(km,kn)是二維高斯粗糙面的功率譜密度，表達(dá)式為

(3)

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

圖1 隨機(jī)粗糙面示意圖

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

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

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

(4)

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

(5)

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

(6)

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

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

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

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

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

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

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

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

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

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

ki·n=0.

(7)

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

(8)

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

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

(9)

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

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

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

圖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相交的投影面元對(duì)應(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所在平面的交點(diǎn)c=(xc,yc,zc).若xc和yc滿足

(13)

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

2.3 改進(jìn)效果

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

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

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

3 隨機(jī)粗糙面散射光強(qiáng)數(shù)值計(jì)算

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

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

圖8 全角度散射光強(qiáng)分布

4 結(jié) 論

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