Infrared remote sensing imaging simulation method for earth’s limb scene

近年来，利用红外卫星临边观测获取的地球临边场景遥感图像对空中运动目标进行探测和识别是遥感技术的重要应用方向之一。临边观测模式得到的图像具有高空间覆盖范围、高垂直分辨率的特点^[1]。然而，公开可获取的图像较少，且难以覆盖不同观测角度、气候条件、云区分布、观测区域的情况。地球临边场景红外仿真可以为空中运动目标检测识别算法的研究提供不同条件下的数据样本。

目前，地球临边场景仿真的相关研究包括海洋场景仿真，大气场景红外辐射特性研究两个方面。针对海洋场景仿真，王彦等^[2]根据海面微面法向量计算海洋辐射特性，渲染海面纹理，计算出红外图像；宋波等^[3]提出了一种高分辨率海洋背景成像仿真方法，通过亚米级海面零视距辐射场，模拟不同海况下的卫星遥感图像。这些研究集中在如何增强海洋红外图像的波纹细节信息，没有考虑到大气环境对海洋辐射模型的复杂影响，不适用于远距离、大范围的卫星遥感领域。针对大气场景红外辐射特性，史波等^[4]建立了临边大气红外辐射计算模型，研究了影响大气红外辐射特性的主要因素；王飞翔等^[5]研究了不同海拔下的大气中波红外透过率。针对大气场景红外辐射的研究主要集中于地面观测角度，不同于卫星遥感领域需要考虑地表出射的红外辐射对大气、云层的影响。

地表包括陆地和海洋，不同于在陆地上可以使用陆地-卫星联合探测，当飞行器目标飞过海洋上空时，对空中目标的观测主要依赖于卫星。深空背景是地球大气层外附近的宇宙空间，其中的红外辐射主要来源于大气层反射的太阳辐射和大气层逸散的辐射。宇宙中除太阳以外的恒星等红外辐射源距离较远，其影响可以忽略不计。因此，文中的研究对象集中在海洋、大气和云。文中研究了海洋和云的红外辐射模型、地球-空间坐标系与相机坐标系的转换关系、大气传输模型，提出了地球临边场景红外遥感成像仿真方法，可以为空中高速目标探测提供有力的数据支撑。

文中主要贡献有：

（1）提出了一种根据云覆盖域的变化计算海洋红外辐射和大气辐射的图像仿真模型，能够反映云区对海洋接收的红外辐射造成的削弱效果；

（2）提出了一种根据空中目标实时位置进行场景定位的图像仿真模型，能够仿真卫星对不同地理位置、不同飞行高度的空中目标进行遥感观测的情况；

（3）根据太阳、地球、卫星探测器的位置变化，实施图像仿真实验，能够快速提供卫星在不同探测条件下的地球临边场景图像。

地球临边场景的红外辐射包含海洋和云的红外辐射。海洋出射的红外辐射包括对太阳辐射的反射，对大气辐射的反射和海洋自发辐射。大气层对红外图像的贡献体现在云的成像、大气传输效应和大气自身热辐射。

地球临边场景图像中每个像元对应的地表面积较大，海洋的纹理细节不明显。和传统的海洋图像仿真算法重视三维形态的仿真不同，地球临边场景中的海洋红外辐射模型主要考察海面在区域范围的辐射和反射特性，其中海洋自发辐射在海洋总体红外辐射中的比重最大。

为了计算海面自身热辐射，可把海面近似作为灰体处理。首先计算同温度下的黑体辐射亮度，再将黑体辐射亮度与海面自发辐射率相乘，得到海面自发辐射亮度：

式中：${\varepsilon }_{sea}=0.98\left[1-{\left(1-\mathrm{cos}\;\theta \right)}^{5}\right]$为海面自发辐射率，$ \theta $表示海平面法向量与探测方向的夹角；$ {L}_{BB}\left(\lambda ,T\right) $为黑体红外波段的自发辐射亮度。

海面同温黑体的光谱辐射出射度为：

式中：${c}_{1}=\left(3.741\;5\pm 0.000\;3\right)\times {10}^{8}\;\mathrm{W}\cdot {\text{μ}\mathrm{m}}^{4}\cdot {\mathrm{m}}^{-2}$为第一辐射常数；${c}_{2}=\left(1.438\;79\pm 0.000\;19\right)\times {10}^{4}\;\text{μ}\mathrm{m}\cdot \mathrm{K}$为第二辐射常数；$\lambda$是辐射波长，单位：μm；T 代表海面的绝对温度，单位：K。在红外波段，黑体自发辐射亮度为：

式中：${\lambda }_{1}$、${\lambda }_{2}$分别为红外波段范围的下限和上限，单位：μm。海洋在$8 \sim 14\;\text{μ}\mathrm{m}$的自发热辐射强度比在$3 \sim 5\;\text{μ}\mathrm{m}$波段大一个数量级左右。

计算海面的反射模型，首先计算海洋接收的红外辐照度，再将接收的红外辐照度与海面在探测方向的反射率相乘，得到海面反射辐射亮度。

利用辐射传输软件MODTRAN计算在不同入射角度、不同地理条件下的太阳光谱辐照度$ {E}_{sun,\lambda } $，从而计算红外波段海面接收的太阳辐照度：

平静的海面对光源的反射包括漫反射和镜面反射。使用双向反射分布函数（Bidirectional Reflectance Distribution Function, BRDF）可以用来描述海面的反射特性。BRDF描述了出射方向上的反射光强相对于入射辐照度的比例$ f $：

式中：$ {\theta }_{i} $为入射天顶角；$ {\theta }_{r} $为出射天顶角；${\varphi }_{i}$为入射方位角；${\varphi }_{r}$为出射方位角；$ \lambda $为入射光波长；${L}_{r}$为出射方向的辐照度，单位：$\mathrm{W}/\left({\mathrm{m}}^{2}\cdot {\rm{sr}}\cdot {\text{μ}}\mathrm{m}\right)$；$ {E}_{i} $为入射方向的辐照度，单位：$\mathrm{W}/\left({\mathrm{m}}^{2}\cdot {\text{μ}}\mathrm{m}\right)$。BRDF函数表示为：

式中：D为坡度分布系数；G为几何衰减因子。D和G的数值与海面的粗糙度、光线的入射方向与反射方向有关^[6-7]。

根据海面单位面积接收的太阳辐照度与海面对太阳光的反射率，求得海面对太阳辐射的反射光强为：

同理可计算海面对大气背景辐射的反射模型。大气背景光来自各个方向，为了简化计算，使用4个方位上的45°天顶角入射的大气辐射强度代替整个大气背景光。根据双向反射分布原理，计算海面对大气背景光的反射光强，结合海面对太阳辐射的反射，计算海洋对红外辐射的反射光强。大气辐射在$8 \sim 14\;{\text{μ}}\mathrm{m}$的自发热辐射强度比在$3 \sim 5\;{\text{μ}}\mathrm{m}$波段大一个数量级。此外，还需要考虑云对太阳辐射的遮挡与对大气辐射的衰减^[8]。

云的红外辐射包括云表层反射的太阳辐射、穿过云层出射的背景辐射、以及云的自发辐射。根据云的分布规律生成多层云图，结合云的辐射基本原理模型，计算红外辐射强度。

云的自发辐射集中在$8 \sim 14\;{\text{μ}}\mathrm{m}$波段，采用灰体近似方法来计算。云的辐亮度公式为：

式中：$ {t}_{c} $为云的透过率；$ {L}_{bb}\left({T}_{cloud}\right) $是与云温度相同的黑体辐射亮度；$ {T}_{cloud} $为云的温度。黑体辐射计算公式同公式（3）。灰体的温度稳定时，吸收率等于发射率，吸收率与透过率之和为1,因此$ 1-{t}_{c} $可以看作是云的发射率。

假设云的密度均匀且各向同性，视线上云的透过率可以通过下式来计算：

式中：$ {\sigma }_{e} $为云的容积消光系数；d为视线方向上的云层厚度。

云对太阳辐射的反射遵循朗伯体公式：

式中：$ {\rho }_{c} $为云的反射系数；E为辐射源照射到云上的辐照度；$ \theta $为太阳辐射传输方向与云表面法向量的夹角。与自发辐射相比，云的反射在零视距总体红外辐射中所占比重较小。

穿过云层出射的背景辐射包括海洋背景辐射和近地大气背景辐射。二者在云层中传输时经过吸收、散射等作用后大幅衰减，能够穿出云层出射的辐射强度非常小，可以忽略不计。

采用中点偏移算法（Midpoint Displacement, MPD）来生成云的外观模型数据场。数据场中，数据分布范围代表云的分布范围，各点的数值代表云层厚度。算法基本过程如下：

（1）给四边形的4个角点赋初始高度值，如图1 (a)所示。

Figure 1.  Schematic diagram of midpoint displacement

（2）在正方形的中点生成一个随机值，中点值为4个角点高度平均值叠加随机量，如图1 (b)所示。

（3）取正方形中点和任意一条边的两个端点，计算高度平均值，叠加随机偏移量，分别得到四条边上的中点值，如图1 (c)所示。

上述过程为0级递归，能将2×2区域扩充为3×3区域。重复上述过程，通过多级递归求得二维数据场d。

将数据场d转换为云的厚度场H，点$ \left(x,y\right) $处的云厚度值H为：

式中：$ {d}_{max} $和$ {d}_{avg} $分别为$ d $的最大值和平均值；$ d\text{'} $是由$ d $计算得到的数据场，$ d\text{'} $中的0值像素是$ d $中小于平均值$ {d}_{avg} $的像素，代表无云点；${H}_{min} \sim {H}_{max}$代表云厚度的分布范围，$ \Delta H={H}_{max}-{H}_{min} $。

由此得到单层云的三维分布信息。地球上的云分布在不同的海拔高度。在高云、中云、低云的高度范围内分别选取代表值作为云层高度，仿真单层云的外观模型，结合云的红外辐射模型计算多个云层叠加后的仿真图像。

地球临边场景红外图像仿真的实质是计算探测器各像素的灰度值。首先根据探测器的地理位置、探测方向，判断探测器各像元的探测对象是海洋或云层，或无对象（深空场景）；其次，根据第1节所述的海洋场景和云场景红外辐射模型，结合红外辐射的大气传输理论计算探测对象的红外辐射在传输路径上的衰减，从而得到探测器各像元接收的红外辐射强度；最后将红外辐射强度信号转换为图像灰度值，得到地球临边场景红外仿真图像。仿真过程如图2所示。

Figure 2.  Schematic diagram of infrared image simulation of earth limb scene

地球-空间坐标系以地心为原点，可以直观地描述海洋与云层的分布情况，卫星的经纬度、高度和太阳直射点的经纬度。卫星坐标系以卫星为原点，光轴方向指向地球中心。相机坐标系以卫星搭载的红外传感器为原点，光轴方向指向空中虚拟目标。通过卫星姿态控制系统^[9]，可以调整传感器方向，将其对准虚拟目标，从而将卫星坐标系转换为相机坐标系。虚拟目标不参与红外成像的过程。

通过计算地球-空间坐标系、卫星坐标系和相机坐标系之间的位置和方向，得到地球-空间坐标系与相机坐标系的转换关系，从而确定地球-空间坐标系中的任意点在相机坐标系中的位置。结合成像投影公式，得到该点在像面的投影位置，即可确定红外相机单个像元的探测对象。在计算中，将面元细分为更小面元，再计算面元成像位置，能够提高仿真图像的准确性^[10]。坐标系之间的转换关系如图3所示，以地球中心（Center of earth and O-earth）为原点的地球-空间坐标系；以卫星（O-satellite）为原点，Z轴指向地球中心的卫星坐标系；以相机（O-camera）为原点，Z轴指向虚拟目标（Virtual target）的相机坐标系。卫星坐标系的中心与相机坐标系的中心可近似视为同一个点。

Figure 3.  Schematic diagram of transformation between coordinate systems

相机坐标系中的点$\left(x,y,{\textit{z}}\right)$在像素坐标系中的位置$ \left(u,v\right) $遵循相机投影公式：

式中：$ H=\left[\begin{array}{ccc}{f}_{x}& 0& {u}_{0}\\ 0& {f}_{y}& {v}_{0}\\ 0& 0& 1\end{array}\right] $表示相机内部参数；f_x、f_y分别表示x轴、y轴上的焦距；（u₀，v₀）表示像面的中心坐标；RT是一个3×4矩阵，描述了地球-空间坐标系与相机坐标系的转换关系；s为归一化常数。

根据地球-空间坐标系、卫星坐标系和相机坐标系之间的转换关系，结合相机投影公式，判断卫星红外传感器单个像元的探测对象。根据探测对象的不同，计算探测器接收的红外辐射强度。当红外卫星传感器的探测对象是海洋时，传感器接收的辐射能量包括海洋表面零视距辐射出射度穿过整个大气层到达大气层外的部分，以及视线路径上各点的大气背景辐射到达大气层外的部分^[11]：

式中：$ {\tau }_{sea} $为海洋辐射在大气传输路径上的衰减率；$ {L}_{h} $为海面自发辐射；$ {L}_{r,sun} $为海面对太阳辐射的反射；$ {L}_{r,sky} $为海面对大气背景辐射的反射；$ {\tau }_{sky} $表示大气中某点在传输路径上的辐射衰减率；$ {E}_{sky} $为该点处的大气背景辐照度。

当探测对象是云时，传感器接收的辐射能量为云层表面零视距辐射出射度和视线路径上各点的大气背景辐射到达大气层外的部分：

式中：$ {\tau }_{cloud} $为云辐射在大气传输路径上的衰减率；$ {L}_{cloud} $表示云的自发辐射；$ {L}_{reflection} $为云对太阳辐射的反射；$ {\tau }_{sky} $和$ {E}_{sky} $与公式（13）中的同名参数意义相同，而数值不同，这是因为不同层流的大气辐射强度有差异 ^[12]；$\int {\tau }_{sky}\cdot {E}_{sky}$中包含单次和多次散射的太阳光谱辐射^[13]。

传感器单像元接收的辐射强度光学信号经过光电转换和灰度映射转化为灰度值：

式中：$ {V}_{s} $表示单像元的电压值；$ {L}_{d\left(i,j\right)} $表示像面上$ \left(i,j\right) $处的像元接收的辐射强度，根据探测对象的不同分别为$ {L}_{s} $（海洋）、$ {L}_{c} $（云）或0 （深空）；$ {R}_{d} $为传感器的光谱响应率；$ {A}_{d} $为像元接收面积。电压经过灰度映射转化为灰度值：

式中：$ F\left(i,j\right) $表示$ \left(i,j\right) $处的像元灰度值；$ {V}_{n} $为随机噪声；$ {V}_{max} $和$ {V}_{min} $分别为附加随机噪声后的传感器最大电压和最小电压；M表示图像的灰度等级。随机噪声包括散粒噪声、热噪声，以及传感器的随机振动造成的像点模糊等。其中散粒噪声和热噪声都具有白噪声的特性，采用均值为0的高斯分布来模拟。

像点模糊效果采用点扩散函数（Point Spread Function, PSF）来计算：

式中：$ p $为光学响应指数；$\left({i}{\mathrm{\text{'}}},{j}{\mathrm{\text{'}}}\right)$表示像面上的点$ \left(i,j\right) $经过点扩散映射的位置。利用点扩散函数对图像进行卷积，可以模拟传感器随机振动的影响。

根据海洋场景和云场景的红外辐射模型计算探测对象的零视距辐射出射度，经过坐标系转换、大气传输效应计算、传感器效应仿真，得到地球临边场景仿真图像。

为了验证仿真模型的真实性和准确性，选取几组不同拍摄时间、探测方向、虚拟目标位置的仿真参数，模拟不同日照条件下的仿真图像。算法运行环境为Win10 64位操作系统，Intel i7-9700CPU，主存32 GB，软件平台为MATLAB 2012b。模型中的探测器像素尺寸为320×256，像元间距为20 μm，探测波段为8~14 μm。卫星处于东经60°上空。仿真参数设置见表1，仿真效果如图4、图5所示。

Figure 4.  Schematic diagram of infrared simulation effect of earth limb scene (Without noise)

Figure 5.  Schematic diagram of infrared simulation effect of earth limb scene (Including noise)

图4所示为不考虑噪声的仿真图像，图5所示为考虑传感器传输噪声与镜头振动的仿真图像。图中可以明显区分海洋、云和深空因红外辐射特性的不同而导致的图像灰度差异。云的图像存在单方向的条纹，离地球边缘越近，条纹现象越明显。在倾斜角度拍摄真实的云图时，存在单方向的条纹特性，而在地球临边观测中，同样采用倾斜角度来观测云，同样出现了条纹，符合观测经验。

使用拉普拉斯算子和（Laplacian Sum, LS）、灰度梯度(Grayscale Mean Gradient, GMG)评价仿真图像的清晰度^[14]。拉普拉斯算子和比较了3×3 pixel邻域内的像素偏差程度。首先使用拉普拉斯算子对图像求卷积，该算子定义如下：

基于拉普拉斯算子和的图像清晰度定义为：

式中：$ \left|G\left(i,j\right)\right| $为像素点$ \left(i,j\right) $处的拉普拉斯算子卷积；M、N为图像的行数和列数；LS值越大，图像信息越丰富，边缘越锐利，图像质量越好。

灰度梯度函数反映了图像的细节和纹理变化特征，定义为：

式中：$ f\left(i,j\right) $为像素点$ \left(i,j\right) $处的图像灰度；GMG值越大，高频信息越丰富，图像越清晰。

对仿真实验结果进行图像清晰度评价，分别计算图像的拉普拉斯算子和、灰度梯度，其结果见表2。不含噪声的仿真结果图像的平均拉普拉斯算子和（LS值）为0.15，平均灰度梯度(GMG值)为0.70。含噪声图像的LS值为0.14，GMG值为0.68。仿真图像边缘易检测，高频信息丰富，画质清晰。传感器噪声对图像质量的影响较小。图像仿真时间见表2第6列，平均单帧图像仿真时间为123.45 s，运行速度较快。

文中研究了地球临边场景的红外卫星遥感成像仿真方法。首先，根据场景组分的差异，分别建立海洋分布模型、多层云分布模型，并根据海洋和云层的红外辐射和反射特性，建立地球临边场景红外辐射模型。其次，建立了地球-空间坐标系与红外相机坐标系的转换关系，结合大气传输理论和传感器效应，计算各观测角度的地球临边场景红外遥感仿真图像。文中提出了根据云覆盖域的变化计算海洋红外辐射和大气辐射的图像仿真模型，能够反映云区对海洋接收红外辐射造成的削弱效果；根据虚拟目标进行场景定位，能够仿真卫星对不同地理位置、不同飞行高度的空中目标进行遥感观测的情况。仿真实验结果表明，该方法能够仿真不同时间下卫星在各轨道高度探测的地球临边场景，图像平均拉普拉斯算子和可达0.15，平均灰度梯度可达0.70，图像画质清晰；单帧图像平均仿真时间为123 s，可以快速模拟地球临边场景。