对于飞机等运动类遥感物体，由于其运动速度高、飞行成本大，获得多种参数下的高光谱图像用于检测研究的难度远远大于地面静物，因此采用仿真的方法获得多种条件下的图像是一种经济、高效可行的方式。高光谱图像仿真技术可以在不用实际飞行的情况下获得遥感器在不同场景、不同实验条件下的近似图像，不仅节约成本，而且数据生成效率高。通用的高光谱图形仿真方法一般从地物反射特性出发，通过数值模拟或者图像融合生成高光谱图像^[1]，仿真方法灵活，输入参数多，国、内外学者在此方面进行了较多的工作。Allen等^[2]提出了AGR、SAIL、Prospect、FCR等模型研究了森林叶冠的图像仿真问题；文献[3-6]从场景地物反射率出发研究背景的图像仿真问题，获得了基于数学模型的地物场景仿真方法。

但传统的高光谱场景图像仿真方法由于涉及的元素较多，仿真链路较长，造成高光谱图像仿真的周期长，同时由于变量的误差传递链条较长仿真结果与实际应用的场景往往难以较好的符合。因此采用实际测量校正得到的背景数据和目标数据，通过数据合成的方式生成高光谱仿真图像，能够产生与实际应用场景符合性更好的高光谱仿真图像。本文基于上述考虑，通过基于实际测量的背景场景，结合仿真获得目标特性数据通过数据重采样、光谱耦合、随机噪声叠加及遥感器相对定标误差校正等计算，生成高光谱仿真图像，可以更加快速、高效的获得高光谱的图像数据，服务于特定场景下的异常检测算法性能评价等应用。

1 图像仿真的数学模型 1.1 飞机辐射特性计算方程

飞行飞机的目标特性可由流动与传热计算控制方程、光谱辐射传输方程^[7]进行计算。飞机流场计算时假设其流动为稳态流动，考虑高温燃气中燃烧产物浓度分布等条件^[8-9]，飞机流场与传热可由下述方程表示：

$ \frac{\partial \rho U_{i}}{\partial x_{i}}=0 $

(1)

$ \frac{\partial \rho U_{i} U_{j}}{\partial x_{i}}=-\frac{\partial P}{\partial x_{j}}+\frac{\partial}{\partial x_{i}}\left(\mu\left(\frac{\partial U_{i}}{\partial x_{j}}+\frac{\partial U_{j}}{\partial x_{i}}\right)-\rho\;\overline{u_{i} u_{j}}\right), i, j=1,2,3 $

(2)

$ \frac{\partial \rho U_{i} T}{\partial x_{i}}=\frac{\partial}{\partial x_{i}}\left(\Gamma \frac{\partial T}{\partial x_{i}}-\rho\;\overline{u_{i} \boldsymbol{\tau}}\right)-\frac{1}{C_{p}} \nabla \cdot q_{\mathrm{R}} $

(3)

$ \frac{\partial \rho U_{i} Y_{i}}{\partial x_{i}}=\frac{\partial}{\partial x_{i}}\left(D_{i} \frac{\partial Y_{i}}{\partial x_{i}}-\rho\;\overline{u_{i} y_{i}}\right) $

(4)

式中：ρ为流体密度，U_i、U_j为流体速度，P为流体压力，μ为流体的动力黏性系数，T为流体温度，Г为流体热扩散系数，C_p为流体定压比热容，u_i、u_j为湍流的脉动速度，τ为湍流的脉动温度，▽·q_R为辐射热源，Y_i为i组分的质量浓度，D_i为该组分的质量扩散系数，$ - \rho\; \overline{{u_i}{y_i}} $为流场中的湍流扩散项。

飞机红外辐射计算通过辐射亮度在参与性介质内的传输方程计算获得，即

$ \begin{gathered} \underbrace{\frac{\mathrm{d} L_{\lambda}(S)}{\mathrm{d} S}}_{(c)}=-\underbrace{\alpha_{\lambda} L_{\lambda}(S)}_{(E x)}+\underbrace{\alpha_{\lambda} L_{\lambda b}(S)}_{(E m)}-\underbrace{\sigma_{s \lambda} L_{\lambda}(S)}_{(E s)}+ \\ \underbrace{\frac{\sigma_{\lambda s}}{4 \pi} \int_{\omega_{i}=4 \pi} L_{\lambda}\left(S, \omega_{i}\right) {\mathit{\boldsymbol{ \boldsymbol{\varPhi}}}}\left(\lambda, \omega, \omega_{i}\right) \mathrm{d} \omega_{i}}_{(S)} \end{gathered} $

(5)

式中：(C)为单位时间内、单位面积、光谱能量在ω方向的单位立体角中传输单位距离后的变化率，(Ex)为被介质吸收引起的光谱能量衰减，(Em)为由于介质发射引起的光谱能量增加，(Es)为被介质散射引起的光谱能量衰减，(S)为各个方向投射在S处的能量的散射引起的ω方向上的光谱能量的增加，Φ(λ, ω, ω_i)为相函数，ω_i方向上的入射辐射亮度引起的ω方向上的光谱能量的增加。

求解上述方程的辐射边界条件与界面的辐射性质有关，对于灰体界面，辐射由壁面自身辐射的光谱辐射亮度与壁面对入射辐射的反射形成光谱辐射亮度共同组成，如下式所示。

$ {L_{\mathit{\boldsymbol{\lambda }},m}}(0) = {\mathit{\boldsymbol{\varepsilon }}_m}{L_{b\mathit{\boldsymbol{\lambda }}}}\left( {{T_m}} \right) + \frac{{1 - {\mathit{\boldsymbol{\varepsilon }}_m}}}{{\rm{ \mathit{ π} }}}H_{\mathit{\boldsymbol{\lambda }},m}^0 $

(6)

式中：ε_m为壁面的发射率，T_m为壁面的温度，L_bλ(T_m)为黑体温度为T_m时的光谱辐射亮度，H_{λ, m}⁰为壁面的入射光谱辐射照度。

1.2 高光谱亚像元图像合成仿真模型

按照上述公式可获得目标特性，假设背景特性已经由观测获得起伏数据，则目标与背景按分辨率混合为亚像元时，可表示为

$ L_{\text {混合}}(x, y, \lambda)=\left\{\begin{array}{l} \boldsymbol{L}_{i} \cdot \boldsymbol{A}_{i}+e_{\text {noise }} \text {, 线性条件 } \\ f\left(\boldsymbol{L}_{i}, \boldsymbol{A}_{i, b}\right)+e_{\text {noise }}, \text { 非线性条件 } \end{array}\right. $

(7)

式中：L_i为高光谱影像中任意一个L维的光谱向量，A=[a₁, a₂, a₃，…，a_i]为对应的组分，e_noise为噪声项。由于实际应用中线性条件下计算应用较为方便^[10]，因此可采用线性混合模型描述目标与背景的耦合，考虑目标在图像中的空间位置分布，混合像元辐射表示为

$ \begin{aligned} &L_{\text {混合}}(x, y, \lambda)=A_{\mathrm{T}}(x, y) \times L_{\mathrm{T}}(\lambda)+ \\ &{\left[1-A_{\mathrm{T}}(x, y)\right] \times L_{B}(x, y, \lambda)+e} \end{aligned} $

(8)

式中：L_混合(x, y, λ)为图像上混合位置的像元值；A_T(x, y)为目标在像元(x, y)中的丰度，对于面积为S_plan的飞机，在遥感器分辨率为R_esol观测条件下，A_T(x, y)=S_plan/Re_sol²；L_T(λ)为目标的光谱；L_B(x, y, λ)为背景的光谱。一般通过数值仿真会获得目标的辐射强度L_T及背景的辐射亮度值H_b^[11]，当考虑不同的丰度混合时，由于目标的尺寸不变，应考虑遥感器分辨率变化带来的背景与目标像元强度的仿真变化，因此仿真图像矩阵表示为

$ L_{\text {合成 }}(x, y, \lambda)=\left[\begin{array}{ccccc} \boldsymbol{\delta}_{\text {new }} H_{\mathrm{b}, 1}+e & \boldsymbol{\delta}_{\text {new }} H_{\mathrm{b} 1,2}+e & \cdots & \boldsymbol{\delta}_{\text {new }} H_{\mathrm{bl}, n-1}+e & \boldsymbol{\delta}_{\text {new }} H_{\mathrm{b} 1, n}+e \\ & \cdots & & & \\ \cdots & & & L_{\text {混和}}+e & \cdots \\ \boldsymbol{\delta}_{\text {new }} H_{\mathrm{bm}, 1}+e & \boldsymbol{\delta}_{\text {new }} H_{\mathrm{bm}, 2}+e & \cdots & \boldsymbol{\delta}_{\text {new }} H_{\mathrm{bm}, n-1}+e & \boldsymbol{\delta}_{\text {new }} H_{b n, n}+e \end{array}\right] $

(9)

式中，δ_new为新的混合像元丰度下A_new下其他像元的校正系数，δ_new=R_esol²×A_T/A_new。

在特定遥感器光谱响应特性下，目标和背景的高光谱数据需要进行校正，统一到遥感器的响应特性下，对目标光谱曲线和背景光谱曲线按照遥感器的光谱响应函数f(λ_i)进行光谱重采样校正采用如下：

$ L_{k}\left(\lambda_{j}\right)=\frac{\int_{\lambda_{j, L}}^{\lambda_{j, U}} \rho_{k}(\lambda) \times f_{j}(\lambda) \mathrm{d} \lambda}{\int_{\lambda_{j, L}}^{\lambda_{j, U}} f_{j}(\lambda) \mathrm{d} \lambda} $

(10)

式中：L_k(λ_j)为目标或者背景光谱k按照波段λ_j的遥感器光谱响应等效后在波段λ_j处的光谱值，f_j(λ)为波段λ_j的遥感器光谱响应函数，λ_j,U、λ_j,L分别为波段λ_j的上、下波长范围，L_k(λ)为k的光谱曲线。

高光谱遥感器由于谱段较多，一般没有完整提供每个波段的光谱响应函数，一般提供中心波段和半带宽(Full width half maximum，FWHM)^[12]，其光谱响应函数通常可以根据中心波长与半带宽通过高斯函数模型进行模拟^[13]。

$ f\left(\lambda_{i}\right)=\frac{1}{\sigma \sqrt{2 \pi}} \mathrm{e}^{\frac{\left(\lambda_{i}-\mu\right)^{2}}{2 \sigma^{2}}} $

(11)

$ \sigma=\mathrm{FWHM} /(2 \times \sqrt{2 \ln 2}) $

(12)

式中：μ为中心波长，σ为标准差，FWHM为光谱分辨率。

实际上，随机噪声、遥感器相对辐射定标误差也将影响成像质量^[14-15]，在上述目标和背景项上叠加噪声随机误差和相对辐射定标误差，得到DN值图像描述为

$ \begin{aligned} \operatorname{DN}(x, y, \lambda)=&\left[\operatorname{gain}_{\mathrm{sig}} \times L_{\text {合成 }}(x, y, \lambda)+\right.\\ &\text { gain } \left._{\text {noise }} \times \mathrm{e}_{\text {nosise }}(x, y, \lambda)\right] \times \\ & n_{\text {cam }}(x, y, \lambda) \end{aligned} $

(13)

式中：DN(x, y, λ)为仿真图像像元DN值，gain_sig为信号的增益，gain_noise为噪声增益，L_合成(x, y, λ)为入瞳辐射值，e_nosise(x, y, λ)为噪声辐射值，n_cam(x, y, λ)为考虑遥感器相对定标的误差系数。

通过目标光谱和背景光谱生成高光谱亚像元仿真图像的基本流程如图 1所示，假设已经通过仿真的方式获得目标和背景经过大气校正的数据，则可采用上述公式开展图像模拟。

2 仿真条件 2.1 目标特性的仿真

仿真前试验测定飞机的平均反射率曲线并与发动机边界条件一并输入仿真模型, 见表 1。

10 km和0.1 km的大气吸收系数如图 2所示，输入经仿真计算获得目标辐射特性如图 3、4所示。

2.2 图像合成仿真

背景起伏数据采用Landsat卫星在轨获取的海背景图像如图 5所示，将准备的图 4所示的背景光谱仿真数据赋予每个背景像元即可得到背景起伏高光谱三维图像矩阵^[16-17]。为模拟在不同像元丰度下的异常检测情况，以飞机目标0.42亚像元丰度为起点，依次序列形成0.12~0.42的丰度图像。其他仿真参数见表 2及图 6所示。

3 仿真结果与异常检测分析

根据亚像元图像仿真输入，以0.004的丰度步进值进行仿真，在每个谱段的高光谱图像序列上一共获得了75个亚像元的目标，按照5×15排列，30 dB信噪比下检测结果如图 7所示。从图 7(a)中可以发现，亚像元的辐射强度较弱，肉眼已经较难分辨，图 7(b)为局部放大图，可以看出亚像元较大丰度像元与周围像元的灰度差异。

经典的检测方法如RX算法、CEM算法分别代表了未知飞机光谱下的异常检测算法和已知光谱下的光谱匹配算法，目前两种算法的仍是高光谱检测中的基础算法。其中RX异常检测算法和CEM异常检测算法分别表示为:

$ \operatorname{RX}\left(\boldsymbol{r}_{x}\right)=\left(\boldsymbol{r}_{x}-\mu_{b}\right)^{\mathrm{T}} C_{b}^{-1}\left(\boldsymbol{r}_{x}-\mu_{b}\right) $

(14)

$ \operatorname{CEM}\left(\boldsymbol{r}_{x}\right)=\boldsymbol{r}_{x}^{\mathrm{T}} \frac{R^{-1} d}{d^{T} R^{-1} d} $

(15)

式中: RX(r_x)为RX检测算子，CEM(r_x)为CEM异常检测算子，r_x为待检测像素光谱向量，μ_b为背景窗口均值，C_b为背景窗口协方差，d为待检测飞机的光谱。

采用RX异常检测方法对结果进行异常检测，获得RX算法的ROC曲线如图 7(c)所示，可以发现部分亚像元目标已经难以在低虚警率下检出。

针对丰度为0.2和0.3的图像分别计算不同信噪比下检出异常目标最小的检测阈值，结果如图 8(a)所示，从图中可以看出，在两种不同的丰度下，异常目标检测所需的最小检测阈值随信噪比的升高而提升，在信噪比较小时，目标检测阈值明显减小。设置信噪比为10 dB，采用RX算法和CEM检测算法^[19-20]开展异常检测，检测结果如图 8(b)、图 8(c)所示，从图中可以看出，RX算法已经无法从仿真图中检出预设的异常像元位置，CEM仍然可以检出异常像元结果。

4 结论

1) 采用仿真目标特性和背景起伏数据通过数值计算可以高效获得多种目标丰度、多信噪比的亚像元仿真图像。

2) 高光谱亚像元图像仿真中噪声会对检测结果产生较大影像，当噪声较大时，可以采用已知目标光谱匹配异常检测提高检测概率。

3) 本文的仿真方法可以用于亚像元异常目标检测中不同类型检测算法性能的比较。