混合像元分解在氧氣吸收被動測距中的應用

周 斌，張 瑜，王海偉，王繼紅

（1. 鄭州科技學院，河南 鄭州 450064；2. 河南省智能信息處理與控制工程技術研究中心，河南 鄭州 450064；3. 軍事科學院系統(tǒng)工程研究院，北京 100010；4. 中國航天系統(tǒng)科學與工程研究院，北京 100045；5. 陸軍裝備部駐沈陽地區(qū)軍事代表局駐長春地區(qū)第一軍事代表室，吉林 長春 130000）

1 引 言

氧氣吸收被動測距系統(tǒng)是利用大氣中氧氣對目標輻射光譜的吸收特性來解算目標距離的一種被動測距系統(tǒng)［1-3］。其中，氧氣A 吸收帶（758～778 nm）較其它兩個主要吸收帶即γ 吸收帶（629～679 nm）和B 吸收帶（685～699 nm）有更高的吸收強度［4］，利用該吸收帶及其兩側帶肩內的光譜獨立性，通過左右?guī)Ъ缈蓴M合出吸收帶內無氧氣吸收時的吸收基線，而后采用多項式擬合等方法獲得氧氣A 吸收帶的平均透過率，通過洛倫茲線型的Elasser 模型［5］、相關K 分布法（Correlated-K Distribution，CKD）［6］等 建 立 氧 氣 平 均 透過率與路徑長度的數(shù)學模型。該方法與其他氣體成分被動測距相比，光譜帶寬僅有50 nm 左右，可對未知輻射目標進行被動測距，具有隱蔽性好、穩(wěn)定度高和測程遠等特點。

實際應用中，無論是被動測距還是主動測距，都可能受到雨雪、霧霾、沙塵等極端天氣的影響［7］。復雜天氣條件下，大氣傳輸路徑上的懸浮粒子，如雨滴、霧滴、霾粒子、雪粒子及沙塵粒子等能夠通過散射和吸收作用［8］，使位于近紅外波段的氧氣A 吸收帶隨傳輸距離的增加而衰減，從而對氧氣吸收被動測距系統(tǒng)的測程和精度產(chǎn)生嚴重的不良影響。特別是，氧氣平均吸收率的精確計算主要依賴于目標輻射光譜的準確提取，這對后續(xù)測距模型的解算尤為重要，但因空氣中各類微小粒子對輻射的散射和吸收等疊加作用，氧氣吸收被動測距系統(tǒng)采集的目標光譜必然夾雜著各類復雜的背景光譜，給目標的光譜提取帶來更多的不確定性，降低了系統(tǒng)的測距能力。雖然采用多次循環(huán)采集平均法［9］、背景消除法［7］可在較大程度上消除被動測距系統(tǒng)噪聲，以及大氣各類氣溶膠粒子散射、吸收作用的影響，但這些方法需要采集多次數(shù)據(jù)且精度并不理想，無法滿足實際飛行目標的實時性測距要求。

混合像元分解技術是隨著高光譜技術而發(fā)展起來的一種圖像處理技術。該技術將大氣中的各類微小粒子視為一種特殊的端元參與到解混過程，利用端元提取算法將解混后的粒子端元去掉，獲取所需的目標光譜信息［10-12］。本文將被動測距系統(tǒng)獲得的高光譜圖像中的像元光譜視為目標與背景的混合光譜，以雨滴粒子為例建立復雜天氣條件下目標與背景的像元混合模型，然后通過混合像元分解，抑制雨滴粒子端元豐度，以期提高測距精度，并選取降雨、霧霾和降雪等3種極端天氣為例對固定點目標進行被動測距實驗，研究該背景抑制方法在氧氣吸收被動測距系統(tǒng)中的有效性。

2 復雜背景條件下光譜混合模型

2.1 光譜圖像像元混合機理

本文所述的被動測距系統(tǒng)是一種高光譜成像光譜儀。它采集的目標輻射光譜信息以像元為單位進行記錄，通過獲取光譜信息可得到目標在一連續(xù)波段的圖像。它在全色二維圖像空間分布基礎上，增加了第三維的光譜特征信息，具有分辨率高、圖譜合一、相鄰波段相關性高的特點。圖1 描述了本文被動測距系統(tǒng)采集的不同波段圖像疊加示意圖，其中單一波段圖像的空間尺寸為1 002 pixel×1 002 pixel，空間維X和空間維Y表示圖像的空間位置，光譜維Z表示圖像的光譜特征。通過采集光譜數(shù)據(jù)可得到任一像元在不同波段的光譜分布，或視場內單一波段的灰度圖像。

圖1 高光譜數(shù)據(jù)示意圖Fig.1 Schematic diagram of hyperspectral data

利用光譜儀對目標輻射進行采集時，因光譜儀的空間分辨率和瞬時視場角以及目標與光譜儀之間遠距離的影響，單個像元內除包括目標輻射光譜外，還包含目標附近其他物體反射或散射的太陽光譜，從而產(chǎn)生混合像元，導致光譜儀獲取的目標光譜特征模糊。若簡單地將這種混合像元歸為單一物體輻射或反射光譜，將導致物體輻射或反射光譜基本特征無法真實反映。利用混合像元分解方法，計算光譜儀圖像中各種亞像元（端元）的物體類別及其所占面積比例，可獲得物體光譜端元。

2.2 目標與復雜天氣背景的混合像元模型

霧霾、煙塵及微小顆粒雨雪等氣溶膠粒子可視為一種特殊端元，利用豐度值可以刻畫出它們對像元中目標光譜的影響［13］。在被動測距的測程內，光譜圖像中的像元光譜可視為目標輻射光譜、霧霾、雨雪等粒子散射光譜的線性混合，據(jù)此建立目標與復雜天氣背景光譜的混合像元模型。

首先不考慮霧霾、雨雪等天氣的影響，設光譜儀單個像元所對應的物體在波長λ上的輻亮度為L(λ)，單個像元的光譜信號X(λ)為：

由式（2）可以看出，混合模型中并無顯性表達出霧霾、雨雪天氣所產(chǎn)生的附加反射、散射光譜以及它們對目標光譜所產(chǎn)生的衰減作用。因此，需要對光譜混合模型進行修改，引入霧霾、雨雪對目標光譜觀測和精確采集的影響。由于大氣傳輸中各類氣溶膠粒子的復雜性和時變性，光子傳輸具有隨機性，如考慮粒子的多次散射，則必須首先計算光子終止位置和隨機遷移路徑，模型復雜，計算量龐大。因此為簡化分析過程，本文參照［13-14］等文獻給出的假設條件，做如下假設：

第一，霧霾、雨雪等復雜天氣條件下，在目標輻射傳輸路徑和光譜采集時間內，霧霾區(qū)域變化緩慢，屬于均勻漫反射體。在雨雪天氣條件下，粒子譜分布在較短時間內保持不變，同時霧霾、雨雪天氣在短時間內對氧氣A 帶附近不同波段的透過率恒定；

第二，目標輻射光在霧霾、雨雪粒子中不存在多次散射現(xiàn)象，霧霾、雨雪等氣溶膠粒子對該波段輻射無吸收。

基于上述假設，針對我國大部分地區(qū)常見的小雨天氣，建立細小雨滴條件下光譜儀接收的像元線性光譜混合模型。將細小雨滴端元引入單像元光譜信號，則式（1）可表示為：

式中：Lrain(λ)為雨滴端元在波長λ處的輻亮度；tj為第j種物質在雨滴區(qū)域中的透過率，按假設條件視為常數(shù)；rrain為雨滴本身的漫反射系數(shù)，滿足t＋rrain=1。

此時光譜儀所形成的改進的像元線性混合模型描述為：

式中：erain(λ)表示“純”雨滴端元的光譜，βrain表示雨滴自身的反射率。由式（4）可以看出，在考慮雨滴光譜對目標光譜的影響后，改進的混合像元模型在形式上沒有發(fā)生變化，故重新定義。將雨滴端元作為物質端元后的新端元集合e'j(λ)（j=1，2，…，m，m＋1），將雨滴端元加入端元集合后形成新的豐度表示為α'j，則式（4）表示的線性混合模型改進為：

改進的線性混合模型仍然滿足端元光譜豐度值非負、豐度值之和為1 的約束條件?？梢姼倪M的混合模型實際上是在光譜圖像像元中引入一個特殊的等效端元，因此，現(xiàn)有的混合像元分解算法能夠直接適用于式（5）的線性混合像元模型。

滿足式（5）的混合光譜信號可表示為：

依此方法建模，將雨端元以霧霾、雪、沙塵等端元代替可得到各類復雜天氣條件下目標與背景的混合像元模型?？紤]到氣溶膠粒子在不同譜分布下的散射和吸收效率并不相同，模型中需要對不同氣溶膠粒子的透過率、漫反射系數(shù)等參量進行重新賦值。

需要說明的是，改進的線性混合模型是在忽略二次散射和多次散射，并假設光傳輸路徑上降水的相對穩(wěn)定基礎上獲得的。從嚴格意義上講，各類氣溶膠粒子對光譜觀測值的影響極為復雜，這種復雜性源于各類粒子物理結構的復雜性。例如當光傳輸距離越來越遠，降水強度越來越大，光路中雨滴粒子越來越多，雨滴空間重疊概率越來越高，光譜穿過雨滴粒子后的散射和吸收衰減越來越強，對光傳輸?shù)挠绊懖蝗莺鲆暎?4-16］。這是一個復雜的問題，這里不做討論。線性混合像元模型旨在從光譜圖像處理角度出發(fā)，利用簡化的線性模型對復雜天氣條件下各類氣溶膠粒子的影響進行較大程度地近似表達，以便于光譜圖像數(shù)據(jù)處理。

3 基于混合像元分解的目標光譜提取

3.1 混合像元分解技術原理

混合像元分解中線性光譜的解混主要分兩個步驟完成：第一步，提取光譜圖像矩陣中的端元，即確定端元矩陣E；第二步，對端元進行豐度反演，即確定豐度矩陣A。端元提取和豐度反演的基礎是確定端元數(shù)目、對高光譜數(shù)據(jù)進行降維并解決端元光譜變異性的問題。線性光譜混合像元的分解流程如圖2 所示，其中的光譜數(shù)據(jù)來自ENVI 軟件數(shù)據(jù)庫。

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

混合像元分解技術包括3 個階段：首先，確定端元數(shù)目；然后對高光譜數(shù)據(jù)進行降維處理，提高線性光譜解混的計算效率；最后對光譜進行解混。

3.1.1 端元數(shù)目確定

光譜端元數(shù)目m的確定是混合像元分解技術的基礎，直接影響到端元提取精度。HFC（Harsanyi，F(xiàn)arrand and Chang）算 法 基 于Ney-Pearson 理論，通過分析光譜圖像相關矩陣與協(xié)方差矩陣產(chǎn)生的特征值的差異性，對端元的數(shù)目進行判定，屬于特征值方法，該算法可有效提高自動化程度。但HFC 算法未考慮噪聲白化對弱信號的影響，造成數(shù)據(jù)丟失，因此進行端元數(shù)目確定前，需對圖像進行噪聲白化，且樣本數(shù)量不是很大時，假設條件不成立［17］。本文采用的噪聲子空間投影法（Noise Subspace Projection，NSP）是針對噪聲白化的優(yōu)化算法［18］，在噪聲估計正確的前提下，尤其是樣本個數(shù)不多或缺少先驗知識的情況下，僅需計算樣本協(xié)方差矩陣就可實現(xiàn)對端元估算，具有較好的穩(wěn)定性和實用性。

3.1.2 數(shù)據(jù)降維

數(shù)據(jù)降維的目的是利用較低維數(shù)數(shù)據(jù)表示原始光譜數(shù)據(jù)的高維特征，除去圖像中大量冗余信息，保留主要數(shù)據(jù)信息，但同時不影響數(shù)據(jù)完整性。對于本文討論的被動測距光譜數(shù)據(jù)，其維數(shù)并不大，直接應用于被動測距計算的波段數(shù)在50 個波段以內；而且數(shù)據(jù)降維是根據(jù)光譜圖像相鄰波段之間存在的相關性進行的，降維極有可能將氧氣吸收帶內或兩側帶肩用于計算吸收率的數(shù)據(jù)去除，不利于后續(xù)對氧氣吸收率的計算，所以本文對光譜數(shù)據(jù)不降維，直接進行后續(xù)端元提取。

3.2 目標光譜提取方法

端元提取算法主要有基于單形體幾何學的像元純度指數(shù)算法（Pure Pixel Index，PPI）、內部最大體積（N-FINDR）算法、頂點成分分析（Vertex component analysis，VCA）算法等，多數(shù)端元提取算法都需要通過PCA 或非負矩陣分解（Non-negative Matrix Factorization，NMF）將 高光譜數(shù)據(jù)的維數(shù)降至N-1 維后進行處理，可較大程度上影響氧氣吸收率的計算［19-21］。順序最大角凸錐（Sequential Maximum Angle Convex Cone，SMACC）算法［22］為計算得到端元光譜提供了自動并且簡便的方法，其端元是通過迭代依次提取的，每迭代一次獲得一個端元，然后計算端元在像元中所占的比例系數(shù)，利用投影變換從各個像元中去除該端元的影響，直到提取所有的端元后停止。相對于VCA 算法而言，該算法的關鍵步驟是判斷光譜圖像中某個像元是否包含此端元以及是否進行斜交投影。目前，SMACC算法作為全自動端元提取方法，已集成到ENVI軟件中，可通過IDL、Python 等軟件進行直接調用。

SMACC 算法能夠從光譜數(shù)據(jù)中直接獲取端元光譜，可有效避免因降維造成的不同波段目標輻射光譜數(shù)據(jù)丟失，最大限度保持目標光譜的完整性和準確性，確保后續(xù)氧氣吸收率的計算精度。

3.3 復雜天氣氧氣吸收被動測距流程

基于混合像元分解技術的背景抑制方法的基本流程如圖3 所示。依據(jù)探測系統(tǒng)所處位置的海拔、大氣溫度和壓強信息，結合探測的天頂角，利用相關K 分布法建立測距模型；利用NSP 法確定獲得的光譜數(shù)據(jù)的端元數(shù)目，然后調用ENVI軟件SMACC 模塊提取目標端元光譜，根據(jù)氧氣吸收率的計算方法計算氧氣吸收率，最后結合測距模型解算目標距離。

圖3 基于混合像元分解技術的被動測距基本流程Fig.3 Flow chart of passive ranging based on pixel unmixing technology

4 實 驗

4.1 測距原理

設實驗測得的大氣中氧氣A 吸收帶的中心波長為λ1，其左右兩側無氣體吸收帶肩對應的中心波長分別為λ2和λ3，光譜儀測得的氧氣吸收帶內平均光譜輻射強度為Pb1，氧氣吸收帶左右?guī)Ъ绲钠骄庾V強度分別為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 型數(shù)字大氣壓力表。

實驗條件分別選取小雨、重度霾、中雪3 種典型極端天氣情況。其中，小雨天氣的地表能見度小于2 000 m，降雨率≤5 mm/h，實測溫度為13 ℃，大氣壓強為1 007.0 hPa，采用300 W 鹵鎢燈作為目標光源，目標距離被動測距系統(tǒng)2 360 m，二者視在天頂角為89.34°；重度霾和中雪的地表能見度小于600 m，重度霾當天大氣溫度為3 ℃、大氣壓強為1 019.0 hPa，中雪當天大氣溫度為1 ℃、大氣壓強為1 018.0 hPa，因目標輻射衰減嚴重，故采用1 000 W 鹵鎢燈，目標距離被動測距系統(tǒng)550 m，二者視在天頂角約為90°。

4.3 實驗結果

根據(jù)視在天頂角、大氣溫度和壓強、被動測距系統(tǒng)海拔信息，依據(jù)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

對于降雨天氣，選擇單次采集獲得的光譜數(shù)據(jù)，首先通過NSP 算法確定雨天采集的光譜數(shù)據(jù)的端元數(shù)，而后利用SMACC 端元提取方法對單次采集的光譜數(shù)據(jù)進行端元提取，獲得鹵鎢燈目標端元和雨粒子的端元，如圖6 所示。

圖6 降雨天氣條件下端元的光譜分布曲線Fig.6 Curves of endmember spectra in rainy weather

將獲得的目標端元光譜代入式（7）和式（8）計算目標輻射的氧氣吸收率，并根據(jù)所建立的CKD 測距模型解算目標距離。表1 給出了測量時間間隔為5 min 的10 次測距解算平均值，單次測量與平均值的相對誤差均在3%以內，顯示了近1 h 內的距離解算值較為穩(wěn)定。表1 還給出了單次采集與8 次循環(huán)采集法獲得的氧氣吸收率解 算出來的目標距離。

表1 降雨天氣目標的解算距離及誤差Tab.1 Resolving ranges and errors of target in rainy weather

由表1 可以看出，端元提取方法和多次循環(huán)采集法的測距精度明顯高于單次采集下的測距精度。這是因為對于單次采集光譜而言，探測器的背景起伏噪聲、散粒噪聲及假信號噪聲等給目標光譜提取帶來較大不確定性；而端元提取后的距離解算精度與多次循環(huán)采集法相比又有較大提升，多次循環(huán)采集法在2 360 m 距離處的解算值為2 190 m，端元提取后的距離解算值為2 280 m，絕對誤差由-170 m 降到-80 m，測距相對誤差由7.20%提高到3.39%。

重度霾和中雪天氣較雨天環(huán)境更為復雜，使得目標輻射傳輸路徑上的氣溶膠光學厚度增大。雖然算法上引入非線性基線可等效消除大氣散射和大氣湍流等因素對平均吸收率的影響，但實際測得的光譜信息中疊加了大量背景輻射光譜，仍需利用混合像元分解進行背景抑制，扣除背景影響。圖7～圖8 分別為霧霾天氣和降雪天氣條件下獲得的不同端元的光譜分布情況。

圖7 霧霾天氣條件下端元的光譜分布曲線Fig.7 Curve of endmember spectra in haze weather

圖8 降雪天氣條件下端元的光譜分布曲線Fig.8 Curve of endmember spectra in snowfall weather

采用與降雨天氣同樣的測距方法，可解算出重度霾和中雪天氣條件下的目標距離，所得結果分別如表2 和表3 所示。表中均給出了原始測距數(shù)據(jù)和利用背景消除法的測距數(shù)據(jù)。

表2 霧霾天氣目標的解算距離及誤差Tab.2 Resolving ranges and errors of target in haze weather

表3 降雪天氣目標的解算距離及誤差Tab.3 Resolving ranges and errors of target in snowfall weather

由表中數(shù)據(jù)可知，對于霧霾天氣，在550 m利用背景消除法能夠較大程度提高目標光譜氧氣吸收率的擬合精度及距離解算精度，使測距誤差下降到11.81%，而經(jīng)過端元提取后，光譜信息中的混合背景光譜得到進一步去除，使距離解算誤差下降到5.81%，測距精度得到提升；對于降雪天氣，在550 m 經(jīng)像元分解法得到的測距誤差與背景消除法相比，從9.09%下降至4.36%，測距精度得到大幅改善。

從實驗結果可以看出，相同距離上同一天氣條件下的測距精度存在一定誤差，且不同距離不同天氣條件下的測距精度也不相同，霧霾天氣的測距精度最小，小雨天氣的測距精度最大。分析其原因主要有以下幾個方面：首先，受大氣中各類粒子復雜多樣性和成像條件的影響，高光譜圖像存在“異物同譜”或“同物異譜”的現(xiàn)象，因此同一目標的純像元光譜并不完全相同；其次，實驗環(huán)境設置上存在波動，如鹵鎢燈目標由市政供電，鎢絲發(fā)熱輻射光強可能在某一時刻發(fā)生細微變化，或光譜儀積分時間較短使大氣輻射光和探測器內部噪聲一定程度影響成像質量；最后，實際氣象條件具有復雜性和不確定性，例如雨霧混合、雨雪混合等，不同粒子對光子的散射吸收差異較大，使得光在不同粒子中的傳輸特性十分復雜且多樣。

即便如此，在上述3 種復雜天氣條件下的實驗結果表明，基于混合像元分解技術的背景抑制方法可有效滿足復雜天氣條件的目標提取和背景抑制，保證目標光譜的準確性。通過去除復雜天氣疊加的背景光譜端元，提高氧氣吸收率和距離解算值的精度，且相對于循環(huán)采集法、背景消除法耗時較長的不足，目標光譜提取方法更適合采集飛行中目標的輻射光譜。

5 結 論

本文根據(jù)目標輻射光譜采集的實時性要求，建立了復雜天氣條件下目標與背景的像元混合模型，給出了基于混合像元分解的目標光譜提取方法，分別在降雨、霧霾、降雪3 種復雜天氣，利用氧氣吸收被動測距系統(tǒng)對固定點鹵鎢燈目標進行了被動測距實驗，根據(jù)CKD 法建立的測距模型解算出目標距離。實驗結果表明：復雜天氣條件下，混合像元分解法可有效保證目標光譜的準確性，抑制復雜天氣背景干擾，提高氧氣吸收率的計算精度，進而降低測距誤差，提高距離解算精度；降雨天氣條件下，在2 360 m 與循環(huán)采集平均法的7.20%測距精度相比，像元分解后測距精度提高到3.39%；霧霾天氣和降雪天氣條件下，在550 m 與背景消除法的11.81%和9.09%相比，像元分解后的測距精度分別提高到5.81%、4.36%；且只需處理單次采集數(shù)據(jù)，可滿足飛行運動目標光譜實時性采集的要求。

需要明確指出的是，基于混合像元分解的目標光譜提取方法，其效果很大程度上依賴于混合像元模型以及混合像元分解算法，因此如何在混合像元模型中引入各類大氣粒子對光譜觀測值的非線性影響，如何通過特征提取降低光譜異變對解混結果的影響，以及如何準確選定樣本區(qū)域獲得更為準確的端元光譜以提高給定混合像元模型下被動測距系統(tǒng)的距離解算精度是需要進一步研究的問題。