张 弛,吴 鑫,谢 建
(西安电子科技大学 物理与光电工程学院,陕西 西安 710071)
随着CCD技术和探测技术的发展,针对目标与背景的光学散射特性的研究被广泛地应用到目标识别[1]、医学检测[2]和计算机图形学[3]等领域中。海面的光学散射特性是目标与背景光学特性研究的重要内容,受到了国内外学者的广泛关注。受到海面背景的影响,目标的光强信息常常淹没于背景中,降低了海上目标探测与识别的正确率与效率。例如,海面对太阳或天空辐射的反射会产生强烈的耀斑,当耀斑进入光学探测系统中时,产生的强烈杂波会严重影响成像效果,造成目标丢失。与光强探测相比,基于偏振信息的探测技术能够获得更多关于目标和背景的轮廓与细节信息。研究海面的偏振特性,可为海背景下的目标探测与识别、目标隐身与反隐身以及海上搜救等应用提供新的手段。
在红外波段,海面对太阳和天空散射所产生的太阳耀光以及海面的自发辐射存在着特定的偏振特性。太阳耀光和海面自发辐射的偏振方向和偏振强度有明显差异。Cooper[4]和Gregoris[5]通过Cox-Munk 波浪斜率模型分析了太阳耀光的偏振特性,并通过实验测量发现在红外偏振成像中,在中波红外波段,探测器接收的辐射主要来自海面耀光;在长波红外波段,接收的辐射主要来自海面的自发辐射。Shaw[6]等的研究表明,海面的偏振发射率会随着观测角度、风速和波长等因素的变化而发生改变,但文中未涉及海面散射光的偏振特性研究。Chang[7]等认为海面自发辐射的偏振度与观测角度以及粗糙度等因素有关,海面越粗糙,偏振度越低。但同样缺乏海面散射光的偏振特性建模和分析。
以往的研究表明,海面的红外偏振辐射特性与风速、海面的几何分布情况以及海水的介电特性等多种因素有关。由于实时监测与现场实验测量等研究手段会耗费大量的人力物力,为了减少研究成本和缩短研究周期,红外偏振成像建模与仿真技术是一种研究海面红外偏振辐射特性的有效手段。国外早在2002年就开始了对红外偏振成像仿真的研究[8-10],国内在偏振双向反射函数(polarized Bidirectional Reflectance Distribution Function,pBRDF)理论[11-12]与材料的偏振特性方面[13-14]取得了一定的进展,但对海面的红外偏振特性研究较少,尤其是红外偏振成像仿真领域的研究则鲜有报道。
在研究海面光学性质的过程中,通常采用Cox-Munk模型[10,15]对海面的斜率分布进行预测。Cox-Munk模型中假设海面为高斯粗糙面,斜率方差、峰度和偏度等参数是基于海拔高度为12.5 m处的风速计算得出,所以当风速过大或过小时,Cox-Munk模型中参数的不确定度增加[16]。因此,在海面的红外偏振成像仿真研究中,建立一种符合实际海情的海面模型具有重要意义。与Cox-Munk模型相比,基于海浪谱和快速傅里叶变换的海面生成方法不仅可以获得海面的斜率信息,还可以获得海面的高度信息。目前,海浪谱主要有Neumann谱、PM谱、Jonswap谱和Elfouhaily谱[17]等。其中,Elfouhaily谱是基于实测数据得到的,并且是一个全波数范围内的海浪谱,因此本文选用Elfouhaily谱。在红外波段,与海面几何尺寸相比,波长较小,需将海面分成许多微小面元。在偏振仿真成像过程中,大多数文献[15,18]均将小面元近似看成光滑平面,仅考虑了镜面反射,忽略了其他方向上的散射,在相同面元数的情况下,这种方法较为粗略。
本文提出将pBRDF模型和偏振自发辐射模型结合起来,建立了一种红外偏振辐射模型。首先,在海面光学特性表征时,本文将每个微小面元近似为有一定斜率分布的微粗糙面,这样更加符合实际情况,也提高了偏振仿真成像的效率,计算更加便捷。然后,通过数值计算得到了不同风速下偏振度随观测角的变化曲线,对比分析了不同风速对计算结果的影响。最后,建立中波红外偏振成像的仿真模型,生成了海面和船舰的红外偏振图像,并结合提出的红外偏振辐射模型对仿真图像进行了理论分析。该研究结果可以为海上目标探测、识别与分类等应用提供理论依据和数据支撑。
本文利用Monte Carlo方法模拟生成海面的高度场。Monte Carlo法的基本思想是对白噪声进行傅里叶变换,在频域对海面的功率谱进行滤波,再进行傅里叶逆变换得到海面的高度场[20]。设海面在x和y方向上的长度分别为Lx和Ly,等间隔离散点数分别为Nx和Ny,相邻点的距离分别为Δx=Lx/Nx和Δy=Ly/Ny。海面上每点的坐标可表示为xm=m·Δx和yn=n·Δy,其中m=-Nx/2+1,…,Nx/2且n=-Ny/2+1,…,Ny/2。每点的海面高度z(xm,yn)可表示为[23]:
(1)
(2)
式中:kmk=2πmk/Lx,knk=2πnk/Ly,N(0,1)表示均值为0,方差为1的高斯分布。方向海浪谱Ψ(kmk,knk)的表达式为[17]:
Ψ(kmk,knk)=
(3)
斜率方差σ2是描述海面粗糙程度的一个重要参数,具体计算公式如下[24]:
(4)
其中:
(5)
(6)
风速分别为3,5,10 m/s时,通过FFT海面生成方法以及Cox-Munk模型计算得到的斜率方差σ2以及绝对误差δ如表1所示。从表1中可以看出,本文FFT海面生成方法计算得到的σ2与经验值吻合得较好,当风速为5 m/s时,绝对误差仅为0.000 2。
表1 不同风速下的斜率方差σ2和绝对误差δ
Tab.1 Slope varianceσ2and absolute errorδwith different wind speeds
Wind speed/(m·s-1)σ2Cox-MunkFFTδ30.018 360.0210.002 6450.028 60.028 80.000 2100.054 20.0520.002 2
图1 海面坐标系和微面元坐标系Fig.1 Coordinate systems of sea surface and single microfacet
BRDF可以表示为散射辐亮度与入射辐照度的比值[25]:
(7)
式中:θi和φi分别是入射光的天顶角和方位角;θr和φr分别是散射光的天顶角和方位角;Lr是散射辐亮度;Ei是入射辐照度。真实的海面是大量微小面元的集合体,受到风速的影响,海表面存在着高低起伏不定的波浪,因此这些微小面元的法线方向各异。菲涅尔定律和斯涅尔定律不能直接应用到海面的BRDF计算中,必须考虑不同微面元对散射光的不同贡献。图1为海面坐标系和微面元坐标系。其中,n是粗糙面的法线,nμ是微面元的法线,nx是粗糙面本地坐标系的x轴;θ是微面元法线与粗糙面表面法线的夹角;β是微面元上法线与入射光的夹角;θ和β可由球面三角几何得到[21]:
cos 2β=cosθicosθr+sinθisinθrcos(φi-φr),
(8)
(9)
偏振的pBRDF在标量的BRDF基础上发展而来,描述了材质对不同偏振方向的光的反射特性,可以表示为[25]:
(10)
式中:Lr和Ei分别是反射辐亮度和入射辐照度的矢量形式。根据微面元理论,通过将标量BRDF函数与4×4 Mueller矩阵M作用可以得到F的具体表达式[25]:
Sh(θi,θr,σ)·M(θi,φi;θr,φr;λ),
(11)
式中Mueller矩阵的各分量Mij可由琼斯矩阵各分量推出,文献[21]给出了Mij的表达式。Sh为遮挡函数,其表达式[26]为:
(12)
(13)
其中:θi和θr分别为入射光和反射光的天顶角,erfc(x)是互补误差函数。考虑到红外波段海水对辐射的衰减很大,经过多次反射并被探测器接收到的漫反射光所占的比重极少,故F中未包含漫反射部分。
来自海面的红外辐射由两部分组成:海面对太阳以及天空的反射辐射Lr和海面的自发辐射Le,不考虑大气对辐射传输的影响,探测器接收到的红外辐射Ld表示[15]为:
Ld=Lr+Le,
(14)
Lr=Lsun,r+Lsky,r,
(15)
其中,Lr由海面对太阳的反射辐射Lsun,r和海面对天空的反射辐射Lsky,r两部分组成。不同波段下,海面对太阳的反射和对天空的反射所占的比重不同。在中波波段,主要是对太阳的反射[22]。因此,式(14)可以简化为:
Ld=Lsun,r+Le,
(16)
其中[18-19]:
Lsun,r=F(θi,φi;θr,φr;λ)·E(θin,φin,λ)·Sin,
(17)
Le=εsea·LBB(λ,T),
(18)
式中:Lsun,r是海面反射的太阳光谱辐亮度,Sin是入射的Stokes矢量,E(θin,φin,λ)为入射到海面上的太阳光谱辐照度,LBB(λ,T) 是与海面温度T相同的黑体辐射亮度,其表达式如下[27]:
(19)
式中:第一辐射常数C1=3.741 6×108(W·μm4·m-2) ,第二辐射常数C2=1.438 8×104(μm·K)。
由于海水对红外辐射的衰减很大,可以将海水看成是不透明的,因此根据基尔霍夫定律,发射与反射是相关的。当海水处于局部热平衡状态时,能量达到守恒状态,吸收和发射相等。海水的发射率εsea与理想黑体发射率εbb的关系[19,28]如下:
εsea=
(20)
式中:dΩi是单位立体角。理想黑体辐射是非偏振的,黑体发射率可以表示为εbb=[1,0,0,0]T。因此可得出海面的红外发射率的表达式[19,28]为:
(21)
式中:F00,F10,F20和F30分别为矩阵F中对应分量。
在传统的光线跟踪方法中,需将二维随机粗糙海面划分为不同的单个微小面元,并将每个微面元近似当作镜面。为了有效地描述海面的纹理特征,需要采用更精确的几何划分方式,这势必会增大计算时间。为了提高计算效率,本文对每个微小面元采用pBRDF模型进行计算,将单个面元散射贡献进行叠加,得出总的接收辐射为:
(22)
式中:A为生成的二维随机海面的面积,M和N分别为在x轴和y轴上的三角面元个数。
由式(14)可知,在红外波段下,探测器接收到的偏振信息受反射的辐射和自发辐射共同影响。当入射方位和探测方位确定后,可分别计算出特定入射几何和观测几何条件下的反射辐射亮度向量Lsun,r和自发辐射亮度向量Le,进一步可计算海面自发辐射的偏振度DOLPe和海面反射辐射的偏振度DOLPr。为了验证理论模型中自发辐射模型和pBRDF模型的有效性,分别计算了海面的DOLPe和DOLPr,结果如图2~图4所示。
图2为通过pBRDF模型和镜面反射模型[6]计算出来的海面自发辐射DOLPe随波长λ(3.0~4.8 μm)的变化曲线。其中,风速为10 m/s,微小面元上的斜率方差σ2=0.052,观测天顶角分别为10°,45°,70°。虚线为文献[6]中的不同观测天顶角下的自发辐射DOLPe随波长变化的曲线,实线为本文通过pBRDF模型计算得到的不同观测天顶角下自发辐射DOLPe随波长变化的曲线。
图2 不同观测天顶角下DOLPe随波长变化的关系Fig.2 DOLPe versus wavelength with different observation zenith angles
由图2可知,相同观测天顶角下镜面反射和pBRDF模型计算得出的DOLPe具有相同的变化趋势,但是数值上略有差异,并且观测天顶角越大,差异越大。产生这种现象的原因是:由于镜面反射模型仅在垂直和水平方向产生偏振,而在pBRDF模型中Mueller矩阵的作用下,海面的自发辐射在垂直和水平方向以及45°和135°方向同时具有偏振分量。在光谱段3.0~4.8 μm,海面自发辐射在垂直和水平方向线偏振光强的差异高于45°和135°方向线偏振光强的差异,从而导致基于pBRDF模型得到的海面自发辐射DOLPe在数值上低于镜面反射模型。
图3为风速分别为3,5,10 m/s下自发辐射DOLPe随观测天顶角变化的关系,仿真波长为4 μm,微小面元上的斜率方差σ2分别为0.021,0.028 8,0.052。由图3可知,通过pBRDF模型仿真得到的DOLPe与文献[6]中结果较为接近,两者的DOLPe都随着观测角的增大而增大。在0°~30°时,自发辐射的DOLPe较低,与风速的相关性较低;当观测天顶角大于30°时,DOLPe随着风速的增大而减小。并且在数值上,通过镜面反射模型得到的DOLPe高于由pBRDF模型得到的DOLPe。这是由于风速增大时,海面粗糙度增大,对自发辐射的起偏效果减弱,海面为镜面时,偏振效果最明显。
图3 不同风速下DOLPe随观测天顶角变化的曲线Fig.3 DOLPe versus observation zenith angle with different wind speeds
图4 不同入射天顶角下DOLPr随观测天顶角变化的曲线Fig.4 DOLPr versus observation zenith angle with different incident zenith angles
图4是入射天顶角分别为10°,30°,45°,70°时,通过pBRDF模型计算得到的DOLPr随观测天顶角变化的曲线。观测方向在入射平面内(即|φr-φi|=180°),风速为10 m/s,微小面元上的斜率方差σ2=0.052。
如图4所示,当入射天顶角分别为30°,45°,70°时,DOLPr峰值分别出现在观测天顶角为76.5°,61.5°,37.5°处,此时反射光达到了完全偏振;而入射天顶角为10°时,DOLPr峰值处的观测天顶角接近90°,但是在峰值处达不到百分百偏振。光在无风海面上的反射可以看作是镜面反射,只有当入射天顶角为海水的布儒斯特角(约为53.5°)时,反射光在θr=53.5°处达到完全偏振。然而当风速为10 m/s时,由于风力的影响,海面可以看作是由倾斜的具有一定粗糙度的微粗糙面元组成的粗糙面,所以入射光在每个微面元上朝着不同方向散射。尽管每个微面元的倾斜角度不同,入射光在不同微面元上的入射天顶角不同,但整体上还是遵循入射天顶角和观测天顶角之和约等于两倍布儒斯特角。
基于本文提出的红外偏振辐射模型,利用Modtran软件计算太阳的辐照度,对海上探测器接收的红外辐射进行仿真,并结合基于海浪谱生成的三维随机海面的高度及斜率信息,模拟生成海面红外偏振图像,如图5~图7所示。仿真参数如表2所示,仿真平台为MATLAB,海面的粗糙度根据风速以及式(4)计算得到。太阳以30°的天顶角以及0°的方位角入射,经过海面各个微小面元反射至位于(200,200,500)处的探测器中,通过每个面元与探测器的相对位置关系即可计算出反射天顶角和方位角。此处,方位角为沿着x轴正方向逆时针方向旋转的角度。
表2 红外海面偏振成像仿真参数
Tab.2 Simulation parameters for infrared polarization imaging of sea surface
参 数设 置海面区域128 m×128 m微面元斜率方差0.052风速10 m/s风向0°海水温度298.15 K仿真波长4 μm太阳天顶角30°太阳方位角0°探测器位置(200,200,500)舰船折射率6.77舰船粗糙度1.134仿真平台MATLAB
图5 海面自发辐射偏振仿真图像Fig.5 Simulated polarimetric infrared images of sea surface without Sun
图6 包含反射辐射和自发辐射的海面偏振仿真图像Fig.6 Simulated polarimetric infrared images of sea surface with Sun
图5为海面自发辐射偏振仿真图像。由于海面自发辐射的圆偏振分量占比极少,所以S3图像不予显示。从图中可以看出海面自发辐射的S0图像较为平滑,由于S0图像是非偏振的红外强度图像,整体对比度较低,由于海面高低起伏所造成的方向非偏振发射率差异而产生的光强分布变化不太明显。偏振分量S1图像和S2图像相对于S0图像拥有更强的对比度,可以看到海面更多的细节,并且S1图像比S2图像的整体灰度更亮,说明在当前观测方向下,海面自发辐射的水平和垂直线偏振光强的差异高于45°和135°线偏振光强的差异。此外,DOLP图像的对比度和细节最为明显,表明利用偏振探测获取海面自发辐射的偏振度图可以提高图像的质量,能够为识别海上目标提供依据。
图6为包含太阳反射和自发辐射的海面红外偏振仿真图像。图6(a)中,太阳耀光现象比较明显。产生这种现象的原因是在当前太阳入射几何条件下,探测器从坐标为(200, 200, 500)处朝向海面中心探测时,太阳光的入射方向和探测器的接收方向关于部分面片的法向量对称,形成镜面反射,此方向反射率较高,从而产生了太阳耀光。在进行海上目标检测时,太阳耀光会干扰海上目标的检测与识别。由于遮挡与阴影函数的影响,从某些面片上反射的偏振光会被附近的面元遮挡,探测器接收到此面元的光强贡献为0,因此这些面元在图像上表现为黑色。偏振图像S1和S2在一定程度上增强了图像的对比度,但是仍然无法消除太阳耀光的影响。从图6(d)中可以看出,DOLP图像有效滤除了太阳耀光,使图像变得更加均匀。通过对图6中的S0,S1,S2和DOLP图像进行方差、信息熵及对比度等图像质量评价的计算,如表3所示。偏振图像的上述3个评价指标均高于非偏振图像S0,其中S1图像的3个指标最大,且对比度最高,图像所包含信息最为丰富,这是由于海面上S1分量即垂直和水平方向线偏振光差异较大,偏振效应较为明显。
表3 图像质量评价结果
为了更加直观地对比海面与人造目标的偏振特性差异,本文对海面以及舰船进行了红外偏振成像仿真。其中,海面使用pBRDF和偏振自发辐射模型进行仿真,而舰船使用典型的P-G模型[21]进行仿真,仿真结果如图7所示。S0为非偏振图像,图像的灰度与红外辐射亮度成正比,从图中可以看出,舰船温度高于海面温度,细节信息较不明显。图7(b)、7(c)和7(d)为偏振图像,其灰度亮暗与温度无关,是偏振度的体现。灰度越高,偏振度越大,灰度越低则反之。从偏振图像中可以看出,与S0相比,舰船和海面的对比度显著提升。并且,偏振图像下,舰船的纹理细节信息更加明显,不同几何朝向的面元在偏振度上具有明显的差异,说明本文提出的偏振模型可以用于仿真海场景,可为海上目标的检测、识别及跟踪等应用提供数据支持。
图7 包含反射辐射和自发辐射的海面目标偏振仿真图像Fig.7 Simulated polarimetric infrared images of sea surface and target with Sun
本文综合考虑了海面的自发辐射与反射效应对探测器接收辐射的影响,通过建立pBRDF模型和偏振自发辐射模型,分别对海面反射辐射和自发辐射的线性偏振度DOLP进行了分析。结合海浪谱和快速傅里叶变换对二维海面进行了几何建模。在此基础上,将海面划分成若干个微面元,结合每个面元的位置信息和斜率信息,仿真了海面的自发辐射偏振图像和红外辐射偏振图像。与传统的光线跟踪方法不同,本文中将每个微面元均当作一个具有粗糙度的粗糙面,而不是近似地当作镜面,从而减少了面元数的划分,在不影响仿真精度的情况下有效提高了仿真效率。仿真结果表明:本文仿真模型与微面元镜面反射模型计算得出的DOLPe相接近,但是本文的模型能计算出自发辐射在45°和135°方向的线偏振辐射,更加真实地反映海面自发辐射的实际情况。另外,本文对仿真生成的海面红外偏振图像进行了质量评价,并将海面与典型人造目标舰船的红外偏振图像进行了对比,证明了红外偏振成像在海上目标探测方面的有效性。