随机粗糙面散射中遮蔽效应算法的改进

宋 哲，徐 哲，庞 勇，马思远，刘春楠

(辽宁师范大学 物理与电子技术学院，辽宁 大连 116029)

自然界中绝大部分物体表面都是随机粗糙面，研究随机粗糙面的电磁散射特性是获取物体表面特征的关键途径，在遥感领域[1]、零件检测领域[2]、医学领域[3]等领域有广泛应用.随机粗糙面散射特性的研究主要有两种方式：数值计算和理论近似.数值计算方法主要有矩量法[4]、有限元法[5]、时域有限差分法[6]等，特点是计算量大，结果较为精确.理论近似方法主要有微扰法[7]、基尔霍夫近似法[8]等，特点是计算量小，易于理解.其中，基尔霍夫近似法计算效率高，且精度与数值计算法接近，得到广泛应用.其核心思想是将粗糙面划分为有限个微面元，每个微面元近似为平面，则光波在微面元上发生菲涅耳反射，由此获得微面元上的本地场，再结合Statton-Chu积分方程得到粗糙面在空间中的散射场.

基尔霍夫近似方法的局限性在于它认为粗糙面上所有的面元都会被照射以及被散射观测点观测到，但事实上，当光波以一定角度入射到粗糙面上时会产生阴影区，使部分面元没有光照射，同时，对于空间某一观测角度也会存在部分面元的反射光因为被其他面元遮挡而无法到达观测点，这种现象就是遮蔽效应.Beckmann首先提出了遮蔽效应的概念，并给出了单站散射的遮蔽函数表达式[9].Smith研究并给出了高斯分布下随机粗糙表面的遮蔽概率函数[10].Wagner在Beckmann单站散射遮蔽函数的基础上提出了双站散射遮蔽函数[11].王安祥等人将实验测量和遗传模拟退火算法相结合，对双向反射分布函数(BRDF)模型中遮蔽函数的参数进行了反演，发现反演得到的遮蔽函数参数能够反映表面粗糙度[12].柳祎等人建立了含有遮蔽函数的偏振度模型，研究了不同粗糙度下金属与非金属红外自发辐射的偏振度，发现粗糙度对非金属表面红外自发辐射偏振度的影响大于金属，可用于区分金属和非金属[13].郭立新等人利用基尔霍夫近似给出了考虑遮蔽效应的分形粗糙表面散射截面近似公式，并讨论了不同分形维数和空间基频下遮蔽效应对双站散射截面的影响[14].马保科等人运用射线追踪法研究了一维粗糙面的遮蔽效应和二次散射，计算了一维粗糙面的散射系数[15].陈明等人利用射线追踪法研究了带限分形粗糙面的遮蔽效应，分析了影响粗糙面遮蔽效应的因素[16].Yoshifumi等人提出一种包含多次散射的粗糙面射线追踪(RSRT)模型，精确计算了散射分布[17].

当前遮蔽效应的研究主要有两种方式：遮蔽概率函数和几何遮蔽.遮蔽概率函数法是通过提出表面任意面元被遮蔽概率的函数来计算遮蔽效应的，计算量小但精确度稍差.几何遮蔽法是通过射线追踪，根据光线与粗糙面的位置关系来逐点判断被遮蔽情况，计算准确但计算量大.本文将传统的几何遮蔽方法进行了改进，在不改变原有精确度的情况下，缩短了几何遮蔽效应的运算时间，并在此基础上采用基尔霍夫近似法分别计算了S波和P波入射到随机粗糙面时的散射场分布.

1 基于Kirchhoff近似法分析随机粗糙面的散射

呈现高斯分布的随机粗糙面在自然界中十分常见，可用线性滤波法来构建随机粗糙面.描述粗糙面的指标主要有相关长度l(二维粗糙面x和y方向的相关长度分别为lx，ly)和均方根高度δ，则随机粗糙面各处的高度分布可表示为[18]

(1)

其中，M、N为x和y方向离散点的数目，Lx、Ly分别为粗糙面在x和y方向的长度，(xm,xn)为离散点坐标，当离散点的间隔为Δx和Δy时，xm=mΔx，yn=nΔy，km、kn为x和y方向的波数，km=2πm/Lx，kn=2πn/Ly，F(km,kn)为傅里叶系数，可以表示为

(2)

N(0,1)为满足高斯分布的一组随机生成数，S(km,kn)是二维高斯粗糙面的功率谱密度，表达式为

(3)

根据式(1)～式(3)，可分别生成M=N=500、lx=ly=3 μm、δ=0.1 μm和δ=0.3 μm的随机粗糙面，如图1所示.

图1 随机粗糙面示意图

当平面波ki照射到随机粗糙面上时会发生散射，根据基尔霍夫近似理论，将粗糙面划分为若干个近似平面的微面元，光波在微面元上发生菲涅尔反射.图2是光波在微面元上反射时入射场与反射场的示意图，θ为微面元上的光波入射角，θr=θ为微面元上的光波反射角，n为微面元的法向量，kr为反射波矢，Ei为入射场，Er为反射场，则微面元上的本地场E(r′)为入射场与反射场的叠加，即

图2 微面元入射场和反射场示意图

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

(4)

r′为微面元的位置矢量.借助于Stratton-Chu积分方程，可得到半球空间内任意位置的散射场[19]:

(5)

其中,r为空间某一观察点的位置矢量，ds′为单个微面元面积，G(r,r′)为格林函数，其远场近似形式为

(6)

其中，r为空间观测点到原点的距离，ks为散射波矢，k为波数.

2 几何遮蔽效应

利用式(4)得到的散射场包含了粗糙面所有面元的贡献，忽略了有些面元可能没有被照射到或无法被观测到的情况，即面元的遮蔽效应，这些面元对散射场是没有贡献的，在入射角或散射角较大时遮蔽效应更加显著，因此应该对粗糙面散射场进行遮蔽修正.

遮蔽效应包括入射遮蔽和散射遮蔽，入射遮蔽指由于粗糙面上的起伏使入射光被遮挡而无法照到的面元，散射遮蔽指由于面元反射的光线被其它面元遮挡而不能直接到达半球空间观测点的面元，二者需要分开讨论.

2.1 入射遮蔽效应

传统的入射几何遮蔽主要是基于射线追踪法来计算被遮蔽的面元，即将光线认为是三维坐标系下的射线，通过追踪射线与微面元是否相交来判断哪些面元是被遮蔽的.如判断微面元1是否被遮蔽的思路是：假设入射光线ki能照到面元1上，计算粗糙面上除面元1外所有面元是否与ki有交点，若有交点，则面元1被遮蔽，否则面元1没有被遮蔽，如图3中面元1被面元2遮蔽.这种方式虽然能够很准确的判断粗糙表面各面元被遮挡的情况，但时间复杂度较高.

图3 面元1被遮挡情况示意图

设入射光线ki的方位角为0°，建立坐标系使xoz面为入射面，ki的方向指向x轴的正方向，z轴的负方向，因此二维随机粗糙面的入射遮蔽效应可以转化为一维粗糙面的遮蔽效应，如图4所示.从图中可以看出，入射光线在粗糙面上会产生阴影区，使阴影中的面元被遮蔽，每段遮蔽面元都起始于一个被照射面元和被遮蔽面元的“临界面元”，如面元1，终止于入射到面元1的光线延长线与粗糙面的交点面元2，这段被遮蔽的面元都位于入射光线的下方，因此判断一个面元是否被遮蔽，只要找出该面元在入射光方向最近的临界面元，再判断该面元是否在入射光线下方即可.

图4 入射遮蔽效应示意图

首先沿x正方向计算入射光线ki与每个微面元法线n是否满足

ki·n=0.

(7)

若式(7)成立，那么该面元即为临界面元.设临界面元的坐标为(x1,z1)，则通过临界面元的入射光线方程为

(8)

其中,kix和kiz为ki在x和z方向的分量，x和z为光线上的点.

设两个相邻临界面元之间任一面元坐标为(x2,z2)，从临界面元起沿x方向依次代入式(8)，若

(9)

2.2 散射遮蔽效应

散射光线的方向可以为半球空间中的任意方向，因此应考虑二维粗糙面的遮蔽.传统二维粗糙面光线追踪法的思路是：判断一个面元是否被遮蔽须遍历除该面元及其相邻面元以外的所有面元，求该面元的散射光线与其他面元是否有交点，若有交点，则该面元被遮蔽；否则没有被遮蔽.

在实际的光线传播过程中，散射光能够经过的面元是有限的，如图5(a)所示，ks是二维粗糙面上任意面元1的散射光，虚线代表散射光可能经过的面元.图5(b)是二维粗糙面在xoy面的投影，其中,网格为粗糙面上各微面元的投影，k′s为ks的投影，阴影方框为k′s经过的投影面元.因此在计算其他面元是否对面元1的散射光线有遮挡只需判断阴影投影面元对应的粗糙面面元对其是否有遮挡，若有遮挡，则面元1 被遮蔽，无需判断粗糙面上除面元1以外的所有面元，可以大大减少计算量.

图5 粗糙面散射遮蔽效应示意图

首先，若散射光线ks与面元1的法向量n满足ks·n<0，说明面元1的散射光不能直接发射到半球空间，面元1被遮挡.若ks·n≥0，则设k′s的斜率kxy为

(10)

图6 搜索k′s经过的面元

设面元1的坐标为(x1,y1,z1)，法向量为n=(nx,ny,nz)，面元1的散射光线为ks=(kx,ky,kz)，与k′s相交的投影面元对应的粗糙面面元2坐标为(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)联立可得到面元1的散射光线ks与面元2所在平面的交点c=(xc,yc,zc).若xc和yc满足

(13)

则交点在面元(x2,y2,z2)内，说明面元1被遮蔽.

2.3 改进效果

选取lx=ly=4 μm、δ=0.5 μm的随机粗糙面，计算入射角为θi=30°，散射角为θs=70°，散射方位角为φs=45°时的遮蔽效应，改进前与改进后的计算时间统计如表1所示，其中,N为x和y方向划分面元数量，tx为原始算法的程序执行时间，ty为改进后算法的程序执行时间，运算平台为i7-7500U CPU，频率为2.90 GHz，内存为8 GB，编程软件为matlab R2019b.

由表1可知，相对传统射线追踪算法，本文改进算法计算遮蔽效应的时间大大缩短，并且随着粗糙面面元数量的增加，改进算法的优势更明显.

表1 改进后的程序执行时间

3 随机粗糙面散射光强数值计算

选取lx=ly=4 μm、δ=0.4 μm的金属铁随机粗糙表面为对象，当入射光波长为λ=0.905 μm时，金属铁的折射率为n=3.12+3.87i.根据式(4)和上述的遮蔽效应，分别模拟计算了S波、P波以30°角入射时，入射面内散射光强的分布，如图7所示，图中散射光强均归一化到S波入射时散射光强的最大值.由图可知本文模拟结果与文献[20]的结果符合较好，说明本文提出的遮蔽效应计算方法是正确的.图8给出了S波、P波30°角入射时半球空间散射光强的分布，可见半球空间中散射光强的分布是关于入射面(散射方位角为0°)对称的，并主要集中在与入射光镜像方向(即散射角为30°，散射方位角为0°)附近，散射角方向半宽度约40°，散射方位角方向半宽度约60°，近似在镜像方向有最大值，P波散射光强小于S波散射光强.

图7 入射面内散射光强分布

图8 全角度散射光强分布

4 结 论

粗糙面在自然界中十分常见，对粗糙面的散射场进行深入研究在目标检测、识别等领域具有重要意义.本文利用线性滤波法生成了不同相关长度和均方根高度的高斯型二维随机粗糙面，采用基尔霍夫近似方法分析了随机粗糙面的散射.针对传统几何遮蔽运算中时间复杂度高的问题，提出了几何遮蔽的改进算法.对于入射遮蔽，判断某一面元是否被遮蔽只要找出该面元在入射光方向最近的临界面元，再判断该面元是否在过临界面元的入射光线下方即可.对于散射遮蔽，判断面元是否被遮挡，只需要判断该面元的散射光线与此散射光在粗糙面上投影所经过的面元是否有交点即可，无需判断粗糙面上所有面元.入射遮蔽的时间复杂度由O(N4)(O为描述时间复杂度的函数)降为O(N2)，散射遮蔽的时间复杂度由O(N4)降为O(N3)，计算效率得到很大提升.模拟计算了S波和P波入射到铁表面时入射面内散射光强的分布，与文献结果进行对比，基本吻合，证明了改进后方法的有效性.计算了随机粗糙面的全角度散射光强，结果表明散射光强分布关于入射面对称，且主要集中在与入射光镜像方向附近，散射角方向半宽度约40°，散射方位角方向半宽度约60°，P波散射光强小于S波散射光强.研究结果能够为目标探测提供理论参考.