混合像元分解在氧气吸收被动测距中的应用

周 斌，张 瑜，王海伟，王继红

（1. 郑州科技学院，河南 郑州 450064；2. 河南省智能信息处理与控制工程技术研究中心，河南 郑州 450064；3. 军事科学院系统工程研究院，北京 100010；4. 中国航天系统科学与工程研究院，北京 100045；5. 陆军装备部驻沈阳地区军事代表局驻长春地区第一军事代表室，吉林 长春 130000）

1 引 言

氧气吸收被动测距系统是利用大气中氧气对目标辐射光谱的吸收特性来解算目标距离的一种被动测距系统［1-3］。其中，氧气A 吸收带（758～778 nm）较其它两个主要吸收带即γ 吸收带（629～679 nm）和B 吸收带（685～699 nm）有更高的吸收强度［4］，利用该吸收带及其两侧带肩内的光谱独立性，通过左右带肩可拟合出吸收带内无氧气吸收时的吸收基线，而后采用多项式拟合等方法获得氧气A 吸收带的平均透过率，通过洛伦兹线型的Elasser 模型［5］、相关K 分布法（Correlated-K Distribution，CKD）［6］等 建 立 氧 气 平 均 透过率与路径长度的数学模型。该方法与其他气体成分被动测距相比，光谱带宽仅有50 nm 左右，可对未知辐射目标进行被动测距，具有隐蔽性好、稳定度高和测程远等特点。

实际应用中，无论是被动测距还是主动测距，都可能受到雨雪、雾霾、沙尘等极端天气的影响［7］。复杂天气条件下，大气传输路径上的悬浮粒子，如雨滴、雾滴、霾粒子、雪粒子及沙尘粒子等能够通过散射和吸收作用［8］，使位于近红外波段的氧气A 吸收带随传输距离的增加而衰减，从而对氧气吸收被动测距系统的测程和精度产生严重的不良影响。特别是，氧气平均吸收率的精确计算主要依赖于目标辐射光谱的准确提取，这对后续测距模型的解算尤为重要，但因空气中各类微小粒子对辐射的散射和吸收等叠加作用，氧气吸收被动测距系统采集的目标光谱必然夹杂着各类复杂的背景光谱，给目标的光谱提取带来更多的不确定性，降低了系统的测距能力。虽然采用多次循环采集平均法［9］、背景消除法［7］可在较大程度上消除被动测距系统噪声，以及大气各类气溶胶粒子散射、吸收作用的影响，但这些方法需要采集多次数据且精度并不理想，无法满足实际飞行目标的实时性测距要求。

混合像元分解技术是随着高光谱技术而发展起来的一种图像处理技术。该技术将大气中的各类微小粒子视为一种特殊的端元参与到解混过程，利用端元提取算法将解混后的粒子端元去掉，获取所需的目标光谱信息［10-12］。本文将被动测距系统获得的高光谱图像中的像元光谱视为目标与背景的混合光谱，以雨滴粒子为例建立复杂天气条件下目标与背景的像元混合模型，然后通过混合像元分解，抑制雨滴粒子端元丰度，以期提高测距精度，并选取降雨、雾霾和降雪等3种极端天气为例对固定点目标进行被动测距实验，研究该背景抑制方法在氧气吸收被动测距系统中的有效性。

2 复杂背景条件下光谱混合模型

2.1 光谱图像像元混合机理

本文所述的被动测距系统是一种高光谱成像光谱仪。它采集的目标辐射光谱信息以像元为单位进行记录，通过获取光谱信息可得到目标在一连续波段的图像。它在全色二维图像空间分布基础上，增加了第三维的光谱特征信息，具有分辨率高、图谱合一、相邻波段相关性高的特点。图1 描述了本文被动测距系统采集的不同波段图像叠加示意图，其中单一波段图像的空间尺寸为1 002 pixel×1 002 pixel，空间维X和空间维Y表示图像的空间位置，光谱维Z表示图像的光谱特征。通过采集光谱数据可得到任一像元在不同波段的光谱分布，或视场内单一波段的灰度图像。

图1 高光谱数据示意图Fig.1 Schematic diagram of hyperspectral data

利用光谱仪对目标辐射进行采集时，因光谱仪的空间分辨率和瞬时视场角以及目标与光谱仪之间远距离的影响，单个像元内除包括目标辐射光谱外，还包含目标附近其他物体反射或散射的太阳光谱，从而产生混合像元，导致光谱仪获取的目标光谱特征模糊。若简单地将这种混合像元归为单一物体辐射或反射光谱，将导致物体辐射或反射光谱基本特征无法真实反映。利用混合像元分解方法，计算光谱仪图像中各种亚像元（端元）的物体类别及其所占面积比例，可获得物体光谱端元。

2.2 目标与复杂天气背景的混合像元模型

雾霾、烟尘及微小颗粒雨雪等气溶胶粒子可视为一种特殊端元，利用丰度值可以刻画出它们对像元中目标光谱的影响［13］。在被动测距的测程内，光谱图像中的像元光谱可视为目标辐射光谱、雾霾、雨雪等粒子散射光谱的线性混合，据此建立目标与复杂天气背景光谱的混合像元模型。

首先不考虑雾霾、雨雪等天气的影响，设光谱仪单个像元所对应的物体在波长λ上的辐亮度为L(λ)，单个像元的光谱信号X(λ)为：

由式（2）可以看出，混合模型中并无显性表达出雾霾、雨雪天气所产生的附加反射、散射光谱以及它们对目标光谱所产生的衰减作用。因此，需要对光谱混合模型进行修改，引入雾霾、雨雪对目标光谱观测和精确采集的影响。由于大气传输中各类气溶胶粒子的复杂性和时变性，光子传输具有随机性，如考虑粒子的多次散射，则必须首先计算光子终止位置和随机迁移路径，模型复杂，计算量庞大。因此为简化分析过程，本文参照［13-14］等文献给出的假设条件，做如下假设：

第一，雾霾、雨雪等复杂天气条件下，在目标辐射传输路径和光谱采集时间内，雾霾区域变化缓慢，属于均匀漫反射体。在雨雪天气条件下，粒子谱分布在较短时间内保持不变，同时雾霾、雨雪天气在短时间内对氧气A 带附近不同波段的透过率恒定；

第二，目标辐射光在雾霾、雨雪粒子中不存在多次散射现象，雾霾、雨雪等气溶胶粒子对该波段辐射无吸收。

基于上述假设，针对我国大部分地区常见的小雨天气，建立细小雨滴条件下光谱仪接收的像元线性光谱混合模型。将细小雨滴端元引入单像元光谱信号，则式（1）可表示为：

式中：Lrain(λ)为雨滴端元在波长λ处的辐亮度；tj为第j种物质在雨滴区域中的透过率，按假设条件视为常数；rrain为雨滴本身的漫反射系数，满足t＋rrain=1。

此时光谱仪所形成的改进的像元线性混合模型描述为：

式中：erain(λ)表示“纯”雨滴端元的光谱，βrain表示雨滴自身的反射率。由式（4）可以看出，在考虑雨滴光谱对目标光谱的影响后，改进的混合像元模型在形式上没有发生变化，故重新定义。将雨滴端元作为物质端元后的新端元集合e'j(λ)（j=1，2，…，m，m＋1），将雨滴端元加入端元集合后形成新的丰度表示为α'j，则式（4）表示的线性混合模型改进为：

改进的线性混合模型仍然满足端元光谱丰度值非负、丰度值之和为1 的约束条件。可见改进的混合模型实际上是在光谱图像像元中引入一个特殊的等效端元，因此，现有的混合像元分解算法能够直接适用于式（5）的线性混合像元模型。

满足式（5）的混合光谱信号可表示为：

依此方法建模，将雨端元以雾霾、雪、沙尘等端元代替可得到各类复杂天气条件下目标与背景的混合像元模型。考虑到气溶胶粒子在不同谱分布下的散射和吸收效率并不相同，模型中需要对不同气溶胶粒子的透过率、漫反射系数等参量进行重新赋值。

需要说明的是，改进的线性混合模型是在忽略二次散射和多次散射，并假设光传输路径上降水的相对稳定基础上获得的。从严格意义上讲，各类气溶胶粒子对光谱观测值的影响极为复杂，这种复杂性源于各类粒子物理结构的复杂性。例如当光传输距离越来越远，降水强度越来越大，光路中雨滴粒子越来越多，雨滴空间重叠概率越来越高，光谱穿过雨滴粒子后的散射和吸收衰减越来越强，对光传输的影响不容忽视［14-16］。这是一个复杂的问题，这里不做讨论。线性混合像元模型旨在从光谱图像处理角度出发，利用简化的线性模型对复杂天气条件下各类气溶胶粒子的影响进行较大程度地近似表达，以便于光谱图像数据处理。

3 基于混合像元分解的目标光谱提取

3.1 混合像元分解技术原理

混合像元分解中线性光谱的解混主要分两个步骤完成：第一步，提取光谱图像矩阵中的端元，即确定端元矩阵E；第二步，对端元进行丰度反演，即确定丰度矩阵A。端元提取和丰度反演的基础是确定端元数目、对高光谱数据进行降维并解决端元光谱变异性的问题。线性光谱混合像元的分解流程如图2 所示，其中的光谱数据来自ENVI 软件数据库。

图2 混合像元分解技术流程Fig.2 Flow chart of pixel unmixing technology

混合像元分解技术包括3 个阶段：首先，确定端元数目；然后对高光谱数据进行降维处理，提高线性光谱解混的计算效率；最后对光谱进行解混。

3.1.1 端元数目确定

光谱端元数目m的确定是混合像元分解技术的基础，直接影响到端元提取精度。HFC（Harsanyi，Farrand and Chang）算 法 基 于Ney-Pearson 理论，通过分析光谱图像相关矩阵与协方差矩阵产生的特征值的差异性，对端元的数目进行判定，属于特征值方法，该算法可有效提高自动化程度。但HFC 算法未考虑噪声白化对弱信号的影响，造成数据丢失，因此进行端元数目确定前，需对图像进行噪声白化，且样本数量不是很大时，假设条件不成立［17］。本文采用的噪声子空间投影法（Noise Subspace Projection，NSP）是针对噪声白化的优化算法［18］，在噪声估计正确的前提下，尤其是样本个数不多或缺少先验知识的情况下，仅需计算样本协方差矩阵就可实现对端元估算，具有较好的稳定性和实用性。

3.1.2 数据降维

数据降维的目的是利用较低维数数据表示原始光谱数据的高维特征，除去图像中大量冗余信息，保留主要数据信息，但同时不影响数据完整性。对于本文讨论的被动测距光谱数据，其维数并不大，直接应用于被动测距计算的波段数在50 个波段以内；而且数据降维是根据光谱图像相邻波段之间存在的相关性进行的，降维极有可能将氧气吸收带内或两侧带肩用于计算吸收率的数据去除，不利于后续对氧气吸收率的计算，所以本文对光谱数据不降维，直接进行后续端元提取。

3.2 目标光谱提取方法

端元提取算法主要有基于单形体几何学的像元纯度指数算法（Pure Pixel Index，PPI）、内部最大体积（N-FINDR）算法、顶点成分分析（Vertex component analysis，VCA）算法等，多数端元提取算法都需要通过PCA 或非负矩阵分解（Non-negative Matrix Factorization，NMF）将 高光谱数据的维数降至N-1 维后进行处理，可较大程度上影响氧气吸收率的计算［19-21］。顺序最大角凸锥（Sequential Maximum Angle Convex Cone，SMACC）算法［22］为计算得到端元光谱提供了自动并且简便的方法，其端元是通过迭代依次提取的，每迭代一次获得一个端元，然后计算端元在像元中所占的比例系数，利用投影变换从各个像元中去除该端元的影响，直到提取所有的端元后停止。相对于VCA 算法而言，该算法的关键步骤是判断光谱图像中某个像元是否包含此端元以及是否进行斜交投影。目前，SMACC算法作为全自动端元提取方法，已集成到ENVI软件中，可通过IDL、Python 等软件进行直接调用。

SMACC 算法能够从光谱数据中直接获取端元光谱，可有效避免因降维造成的不同波段目标辐射光谱数据丢失，最大限度保持目标光谱的完整性和准确性，确保后续氧气吸收率的计算精度。

3.3 复杂天气氧气吸收被动测距流程

基于混合像元分解技术的背景抑制方法的基本流程如图3 所示。依据探测系统所处位置的海拔、大气温度和压强信息，结合探测的天顶角，利用相关K 分布法建立测距模型；利用NSP 法确定获得的光谱数据的端元数目，然后调用ENVI软件SMACC 模块提取目标端元光谱，根据氧气吸收率的计算方法计算氧气吸收率，最后结合测距模型解算目标距离。

图3 基于混合像元分解技术的被动测距基本流程Fig.3 Flow chart of passive ranging based on pixel unmixing technology

4 实 验

4.1 测距原理

设实验测得的大气中氧气A 吸收带的中心波长为λ1，其左右两侧无气体吸收带肩对应的中心波长分别为λ2和λ3，光谱仪测得的氧气吸收带内平均光谱辐射强度为Pb1，氧气吸收带左右带肩的平均光谱强度分别为Pb2和Pb3，通过插值可得到氧气吸收带内非吸收基线的辐射强度Pb为：

利用文献［6］给出的基于相关K 分布法建立的氧气平均吸收率与路径长度的关系模型，可解算出被测红外目标距离。

4.2 实验设备及条件

光谱测量设备采用高光谱成像光谱仪，其光谱响应范围为450～800 nm，光谱分辨率为1 nm，实验中设置采集波段为740～790 nm，步长为1 nm，带宽为3.6 nm，积分时间为0.5 s；距离标定设备为某型军用激光测距机（精度为±1 m）；大气强度及温度测量设备为BY-2003P 型数字大气压力表。

实验条件分别选取小雨、重度霾、中雪3 种典型极端天气情况。其中，小雨天气的地表能见度小于2 000 m，降雨率≤5 mm/h，实测温度为13 ℃，大气压强为1 007.0 hPa，采用300 W 卤钨灯作为目标光源，目标距离被动测距系统2 360 m，二者视在天顶角为89.34°；重度霾和中雪的地表能见度小于600 m，重度霾当天大气温度为3 ℃、大气压强为1 019.0 hPa，中雪当天大气温度为1 ℃、大气压强为1 018.0 hPa，因目标辐射衰减严重，故采用1 000 W 卤钨灯，目标距离被动测距系统550 m，二者视在天顶角约为90°。

4.3 实验结果

根据视在天顶角、大气温度和压强、被动测距系统海拔信息，依据CKD 法可建立不同实验条件下的测距模型，其中降雨天气测距模型如图4 所示。图5 分别给出了小雨、重度霾和中雪3 种典型极端天气条件下目标750 nm波长的场景图。

图4 降雨条件下的被动测距模型Fig.4 Distance measurement model under rainfall condition

图5 极端天气条件目标在750 nm 波长下的高光谱图像Fig.5 Hyperspectral images of target at 750 nm in extreme weather conditions

对于降雨天气，选择单次采集获得的光谱数据，首先通过NSP 算法确定雨天采集的光谱数据的端元数，而后利用SMACC 端元提取方法对单次采集的光谱数据进行端元提取，获得卤钨灯目标端元和雨粒子的端元，如图6 所示。

图6 降雨天气条件下端元的光谱分布曲线Fig.6 Curves of endmember spectra in rainy weather

将获得的目标端元光谱代入式（7）和式（8）计算目标辐射的氧气吸收率，并根据所建立的CKD 测距模型解算目标距离。表1 给出了测量时间间隔为5 min 的10 次测距解算平均值，单次测量与平均值的相对误差均在3%以内，显示了近1 h 内的距离解算值较为稳定。表1 还给出了单次采集与8 次循环采集法获得的氧气吸收率解 算出来的目标距离。

表1 降雨天气目标的解算距离及误差Tab.1 Resolving ranges and errors of target in rainy weather

由表1 可以看出，端元提取方法和多次循环采集法的测距精度明显高于单次采集下的测距精度。这是因为对于单次采集光谱而言，探测器的背景起伏噪声、散粒噪声及假信号噪声等给目标光谱提取带来较大不确定性；而端元提取后的距离解算精度与多次循环采集法相比又有较大提升，多次循环采集法在2 360 m 距离处的解算值为2 190 m，端元提取后的距离解算值为2 280 m，绝对误差由-170 m 降到-80 m，测距相对误差由7.20%提高到3.39%。

重度霾和中雪天气较雨天环境更为复杂，使得目标辐射传输路径上的气溶胶光学厚度增大。虽然算法上引入非线性基线可等效消除大气散射和大气湍流等因素对平均吸收率的影响，但实际测得的光谱信息中叠加了大量背景辐射光谱，仍需利用混合像元分解进行背景抑制，扣除背景影响。图7～图8 分别为雾霾天气和降雪天气条件下获得的不同端元的光谱分布情况。

图7 雾霾天气条件下端元的光谱分布曲线Fig.7 Curve of endmember spectra in haze weather

图8 降雪天气条件下端元的光谱分布曲线Fig.8 Curve of endmember spectra in snowfall weather

采用与降雨天气同样的测距方法，可解算出重度霾和中雪天气条件下的目标距离，所得结果分别如表2 和表3 所示。表中均给出了原始测距数据和利用背景消除法的测距数据。

表2 雾霾天气目标的解算距离及误差Tab.2 Resolving ranges and errors of target in haze weather

表3 降雪天气目标的解算距离及误差Tab.3 Resolving ranges and errors of target in snowfall weather

由表中数据可知，对于雾霾天气，在550 m利用背景消除法能够较大程度提高目标光谱氧气吸收率的拟合精度及距离解算精度，使测距误差下降到11.81%，而经过端元提取后，光谱信息中的混合背景光谱得到进一步去除，使距离解算误差下降到5.81%，测距精度得到提升；对于降雪天气，在550 m 经像元分解法得到的测距误差与背景消除法相比，从9.09%下降至4.36%，测距精度得到大幅改善。

从实验结果可以看出，相同距离上同一天气条件下的测距精度存在一定误差，且不同距离不同天气条件下的测距精度也不相同，雾霾天气的测距精度最小，小雨天气的测距精度最大。分析其原因主要有以下几个方面：首先，受大气中各类粒子复杂多样性和成像条件的影响，高光谱图像存在“异物同谱”或“同物异谱”的现象，因此同一目标的纯像元光谱并不完全相同；其次，实验环境设置上存在波动，如卤钨灯目标由市政供电，钨丝发热辐射光强可能在某一时刻发生细微变化，或光谱仪积分时间较短使大气辐射光和探测器内部噪声一定程度影响成像质量；最后，实际气象条件具有复杂性和不确定性，例如雨雾混合、雨雪混合等，不同粒子对光子的散射吸收差异较大，使得光在不同粒子中的传输特性十分复杂且多样。

即便如此，在上述3 种复杂天气条件下的实验结果表明，基于混合像元分解技术的背景抑制方法可有效满足复杂天气条件的目标提取和背景抑制，保证目标光谱的准确性。通过去除复杂天气叠加的背景光谱端元，提高氧气吸收率和距离解算值的精度，且相对于循环采集法、背景消除法耗时较长的不足，目标光谱提取方法更适合采集飞行中目标的辐射光谱。

5 结 论

本文根据目标辐射光谱采集的实时性要求，建立了复杂天气条件下目标与背景的像元混合模型，给出了基于混合像元分解的目标光谱提取方法，分别在降雨、雾霾、降雪3 种复杂天气，利用氧气吸收被动测距系统对固定点卤钨灯目标进行了被动测距实验，根据CKD 法建立的测距模型解算出目标距离。实验结果表明：复杂天气条件下，混合像元分解法可有效保证目标光谱的准确性，抑制复杂天气背景干扰，提高氧气吸收率的计算精度，进而降低测距误差，提高距离解算精度；降雨天气条件下，在2 360 m 与循环采集平均法的7.20%测距精度相比，像元分解后测距精度提高到3.39%；雾霾天气和降雪天气条件下，在550 m 与背景消除法的11.81%和9.09%相比，像元分解后的测距精度分别提高到5.81%、4.36%；且只需处理单次采集数据，可满足飞行运动目标光谱实时性采集的要求。

需要明确指出的是，基于混合像元分解的目标光谱提取方法，其效果很大程度上依赖于混合像元模型以及混合像元分解算法，因此如何在混合像元模型中引入各类大气粒子对光谱观测值的非线性影响，如何通过特征提取降低光谱异变对解混结果的影响，以及如何准确选定样本区域获得更为准确的端元光谱以提高给定混合像元模型下被动测距系统的距离解算精度是需要进一步研究的问题。