針對高光譜端元提取的空譜聯合預處理方法

吳銀花,王鵬沖,吳慎將,張發強

(1.西安工業大學 光電工程學院，陜西 西安 710021;2.中國科學院 光譜成像技術重點實驗室，陜西 西安 710029)

1 引 言

高光譜遙感實現了遙感圖像光譜分辨率的突破性提高，它利用很多較窄的波段成像，將觀測到的各種地物與電磁波相互作用后的輻射信號以完整的光譜曲線記錄下來，使得本來在常規全色、多光譜遙感中不能識別的地物，在高光譜遙感中能夠得到有效的識別[1-9]。因此，高光譜遙感在軍事、醫療、農業、海洋、考古等多個領域得到了廣泛的應用，主要用于實現目標分類、目標識別、目標探測等。然而，成像光譜儀往往具有較大的瞬時視場角，使得獲取的高光譜圖像空間分辨率較低，導致存在大量的混合像元，即一個像元所對應的地面區域內往往包含兩種或更多種類的特征地物。這嚴重阻礙了高光譜遙感應用精度的提高。而大部分高光譜目標分類、識別、探測算法的前提是假設每個像元所對應的地面區域內只包含一種特征地物。即每個像元是純像元，又稱為端元。因此，為了提高高光譜遙感應用的精度，就必須進行高光譜解混合，使高光譜遙感應用由像元級達到亞像元級[10-12]。

高光譜解混合通常分為兩個步驟：端元提取和豐度反演。即先從整個高光譜圖像中提取含某一種地物的比例較高的像元作為端元，再針對每個像元求解對應各端元的豐度系數。其中如何從高光譜數據中準確提取端元是高光譜解混合的關鍵，也一直是國內外學者的熱點研究問題。目前已經有不少成熟的端元提取方法被提出和引用，其中利用高光譜數據凸面幾何特性的方法由于符合光譜混合過程的物理原理，且形式簡單、技術成熟，目前得到了非常廣泛的應用，比較經典的算法有頂點成分分析(VCA)[13]、順序最大角凸錐(SMACC)[14]等。然而，這些方法主要是從光譜特征角度出發，進行端元提取，而忽視了像元在空間上的相關性，易受噪聲和異常點的干擾，端元提取精度不高，且每個端元的提取過程都需遍歷圖像中的所有像元，導致運行時間較長。后來也有一些結合光譜特性和空間特性的端元提取方法[15-16]被提出，然而這些算法大多未充分利用真實地物空間分布特性，要么仍然易受噪聲和異常點的干擾，要么還需要為不同圖像人為設定不同閾值來輔助完成端元提取，導致算法的魯棒性較低。

為了提高高光譜端元提取效率，提出了一種空譜聯合預處理方法。該預處理方法先采用文中定義的光譜純度指數SPI，從光譜特征和空間特征兩方面綜合預估高光譜圖像中每個像元成為純像元的可能性；之后基于SPI值在高光譜圖像中的分布，利用空間相關性，有效去除光譜純度較低的冗余像元，從而形成精簡的候選端元集。該候選端元集與常用端元提取算法結合，不僅可提高抗干擾性和端元提取精度，同時可大幅降低時間復雜度。

2 算法介紹

2.1 光譜純度指數SPI

由于真實地物在圖像空間分布的連續性，相鄰像元間總存在一定的相關性，尤其是空間分辨率較高的高光譜圖像。從像元的角度來說，某一地物對應的純像元周圍一般是同一地物對應的純像元或者是含有該種地物成分的混合像元，而不同地物對應的純像元一般不會是相鄰像元。即純像元分布在同質區域中，而同質區域顯然是局部區域，局部區域內可認為同一地物的光譜特征一致。因此同質區域中相鄰像元之間光譜特征往往比較相似，越靠近純像元，越相似。

在給定的局部區域內，某一像元與其相鄰像元的光譜特征越相似，說明該局部區域的同質性越高，從而該像元成為純像元的可能性越高。而光譜特征主要體現在其光譜形狀和幅值：(1)同一地物對應的兩個純像元，二者光譜形狀和幅值均很相近，光譜特征相似度很高；(2) 不同地物對應的兩個純像元，二者的光譜形狀和幅值均差異較大，光譜特征相似度很低；(3) 不同地物邊界中的兩個混合像元，由于兩個像元中含有各種地物的比例不同，二者的光譜形狀和幅值一般很難都相近，光譜特征相似度較低。

光譜角度填圖法(SAM)是通過計算兩個像元光譜之間的夾角θ來衡量二者的光譜特征相似度，如式(1)所示。θ值越小，說明二者越相似。θ值與光譜向量的模無關，因此SAM法主要測量光譜形狀的差異。歐式距離法(ED)是通過計算兩個像元之間的歐式距離d來衡量二者的光譜特征相似度，如式(2)所示。d值越小，說明二者越相似。與SAM法不同，ED法主要測量的是各波段幅值差異的累積。

(1)

(2)

式中：X和Y分別表示兩個像元的光譜向量，xi和yi分別表示X和Y的第i個波段幅值，p表示波段數，θ的取值范圍是[0,π/2]。

綜合考慮光譜形狀和幅值，這里定義了一種新的參數s來測量像元間光譜特征相似度，進而給出了光譜純度指數SPI的概念及其計算方法，用于表示每個空間像元稱為純像元的可能性。具體方法如下：

(1)新的參數s計算方法如式(3)所示。其中，θ′和d′分別是利用式(1)、(2)求出的θ和d經歸一化后的值，以避免角度θ和距離d的數值差異影響。由于參數s既考慮了光譜形狀差異，又考慮了光譜幅值差異，且將這兩種差異視為同樣重要，從而更能準確衡量光譜特征相似度。與θ和d類似，s值越小，說明兩個像元光譜形狀和幅值的綜合差異越小，光譜特征越相似。

(3)

(2)在新的參數s基礎上，定義了光譜純度指數SPI的概念：當前像元Xi,j(i和j分別表示該像元在圖像空間中的橫坐標和縱坐標)與其為中心的N×N大小窗口T內相鄰像元之間參數s的最大值定義為當前像元的光譜純度指數SPIi,j，計算方法如式(4)所示。其中，i和k表示圖像中像元的橫坐標，j和l表示圖像中像元的縱坐標，k的取值范圍是[i-N/2+1,i+N/2-1]，l的取值范圍是[j-N/2+1,j+N/2-1]。

SPIi,j=max{s(Xi,j，Xk,l)}

.

(4)

可以看出，SPI值越大，說明對應像元的局部區域窗口同質性越低，也就是光譜純度越低，從而成為端元光譜的可能性也越低；反之，SPI值越小，說明對應像元的局部區域窗口同質性越高，也就是光譜純度越高，從而成為端元光譜的可能性也越高。

2.2 基于SPI的空間去冗余方法

進行光譜解混合時，對于每種地物，只需要此種地物中一個純度最高的像元作為端元即可。如前所述，由于真實地物在圖像空間分布的連續性，純像元往往分布在圖像的同質區域中，且越靠近純像元，相鄰像元間光譜特征越相似，像元的光譜純度指數SPI值越小。即每個局部區域中，光譜純度指數SPI最小值對應的像元成為純像元的可能性最高，而其他像元由于光譜純度指數SPI值較大，成為純像元的可能性較低，可以認為是冗余像元。顯然，這些冗余像元在端元提取過程中，不僅會提高時間復雜度，還會增加端元提取的錯誤率。因此，需要一種有效的去冗余方法，在保持光譜多樣性不變的基礎上，準確判斷冗余像元，并將其去除，使得僅將成為純像元可能性較高的像元參與端元提取，從而節省運算時間，也降低端元提取的錯誤率。

通過利用真實地物在圖像空間分布的連續性，這里給出了一種基于SPI的空間去冗余方法。根據SPI值的空間分布，尋找每個局部區域中SPI最小值對應的像元，以此來判斷每個像元是否為冗余像元，并將冗余像元去除，形成最終參與端元提取的候選端元集。具體方法如下：

(1)利用式(3)和(4)，計算高光譜圖像中所有像元的SPI值，并生成與圖像同等大小的二維SPI矩陣，如式(5)所示，其中m和n分別表示圖像中像元行數和列數。

(5)

(2)針對每個像元，以該像元為中心的Q×Q大小窗口L范圍內，尋找SPI最小值。這一步可采用Q×Q大小窗口的最小值濾波器，對(1)中獲得的二維SPI矩陣進行濾波的方式實現，既簡單又快捷。

(3)如果某一像元的SPI值與其為中心的窗口L內SPI最小值相等，則認為該像元在窗口L內成為純像元的可能性最高，判定為非冗余像元；反之，判定為冗余像元。由所有非冗余像元組成候選端元集。

由于候選端元集是由每個局部區域中SPI值最小的像元，也就是成為純像元可能性最高的像元組成，從而有助于避免噪聲和異常點的干擾，以提高端元提取精度。同時由于參與端元提取運算的像元數量有效減少，可大幅縮短運算時間。

2.3 基于本文預處理方法的高光譜端元提取

根據2.1～2.2節的描述，結合提出的空譜聯合預處理方法，高光譜端元提取整體框架如圖1所示，具體步驟如下：

圖1 基于預處理方法的高光譜端元提取框架Fig.1 Framework of the hyperspectral endmember extraction with the proposed preprocessing method

步驟1：根據2.1節中的方法，利用式(1)～(4)，以5×5窗口大小，計算高光譜圖像中每個像元的SPI值。

步驟2：基于步驟1中計算得到的SPI值空間分布，利用2.2節中的方法，以3×3窗口大小，進行空間去除冗余，并生成候選端元集。

步驟3：采用常用端元提取方法，在步驟2中獲得的候選端元集范圍內進行端元提取處理。

3 實驗數據

為了驗證該算法的可行性和有效性，本文將分別利用模擬高光譜圖像和真實高光譜圖像進行相關實驗。

3.1 模擬高光譜圖像

模擬高光譜圖像是由USGS光譜數據庫中任意選取的5條光譜曲線作為端元光譜，以一定比例線性混合生成的高光譜數據立方體，空間大小為60像元×60像元。作為端元選取的每條光譜曲線波長范圍為0.395 1～2.56 μm，包含420個波段，光譜分辨率為5 nm。為了模擬真實地物在圖像空間的分布特征，本文按以下步驟生成了模擬高光譜圖像：

(1)將60像元×60像元大小的空間區域分割為25塊12像元×12像元大小的子塊，每個子塊隨機設定為某一種端元，此時對于每個像元，只有設定端元對應的豐度系數為1，其余端元對應的豐度系數為0。

(2)采用窗口大小為15×15、均值為0的高斯低通濾波器，對每一種端元對應的豐度系數圖進行濾波，從而模擬混合像元。

(3)進行歸一化，使得每個像元對應豐度系數之和為1。

(4)從USGS光譜數據庫中另選取2條光譜曲線，同樣波長范圍為0.395 1～2.56 μm，包含420個波段，光譜分辨率為5 nm。

(5)在60像元×60像元大小的空間區域中，隨機選擇兩個像元位置，并分別用步驟(4) 中選取的2條光譜曲線替換這兩個像元位置所對應的原始混合光譜，從而模擬2個異常點。

(6)在步驟(5) 得到的模擬數據中，添加信噪比SNR為30 dB的零均值高斯噪聲，以模擬實際數據采集時生成的各種噪聲和干擾。

圖2是按照上述步驟生成的模擬高光譜圖像第30波段對應的二維灰度圖，其中藍圈里的紅點表示添加的兩個異常點，空間位置坐標分別是(8,28)和(23,35)。圖3是模擬高光譜數據中端元和異常點對應的光譜曲線，其中5條實線對應端元光譜，2條虛線對應異常點光譜。

圖2 模擬高光譜圖像(第30波段)Fig.2 Simulated hyperspectral image (Band 30)

圖3 模擬高光譜圖像中端元和異常點的光譜曲線Fig.3 Spectral curves of the endmembers and outliers in the simulated hyperspectral image

3.2 真實高光譜圖像

選用的真實高光譜圖像是1994年采集的加州Jasper Ridge自然保護區，具有224個波段，波段范圍為380～2 500 nm，光譜分辨率高達9.46 nm，移除受水蒸氣和大氣影響的26個波段(1～3 nm、108～112 nm、154～166 nm和220～224 nm)后，剩余198個有效波段數據。由于該高光譜圖像比較復雜，截取以原始圖像中第(105,269)像元為第一像元的100像元×100像元的子圖像進行端元提取性能評估實驗，該子圖像對應空間區域主要由樹、水、土壤和道路4種端元成分組成。圖4是由該子圖像的第50、75、125有效波段組成的偽彩色圖。

圖4 Jasper Ridge偽彩色圖Fig.4 False color map of Jasper Ridge

4 實驗結果與討論

本文共進行了3組仿真實驗，第1、2組是利用模擬高光譜圖像，第3組是利用真實高光譜圖像。第2、3組實驗中采用經典方法VCA和SMACC進行了預處理之后端元提取，并分別用SPI+VCA和SPI+SMACC來表示。

4.1 第1組實驗

為了驗證候選端元集的有效性，針對3.1中模擬高光譜圖像，先利用2.1的方法計算了其SPI值分布圖，接著利用2.2的方法去除了空間冗余像元，實驗結果如圖5所示。

圖5 候選端元分布圖Fig.5 Distribution of candidate pixels

圖5中紅點表示經過2.1～2.2方法得到的候選端元；而背景灰度圖中各像元的灰度對應其最大豐度系數值，用于表示實際光譜純度，灰度值越大(越亮)，表示該像元實際光譜純度越高。

從圖5可以看出：(1) 大部分候選端元分布在各端元對應的空間區域中光譜純度較高的中心區域，且兩個異常像素點均未被選為候選像元；(2) 僅小部分候選端元分布在不同端元對應空間區域的相鄰邊界，說明這些候選端元雖然實際光譜純度低，但與相鄰像元光譜相似度較高，因此一般都是多個端元均勻混合后的混合像元，進而后續在候選像元集范圍內進行端元提取處理時，較容易去除；(3) 候選端元數量是116，遠少于原始像元數量3 600，大幅減少了后續端元提取過程中時間復雜度。實驗結果表明，提出的候選端元獲取方法合理有效。

4.2 第2組實驗

為了進一步驗證候選端元集的有效性，比較了未經任何預處理時端元提取效果與采用提出的預處理方法時端元提取效果。分別利用VCA和SPI+VCA、SMACC和SPI+SMACC算法，對模擬高光譜圖像進行了端元提取處理。不同算法提取的端元對應空間位置坐標如表1所示，對應光譜曲線如圖6所示。

表1 提取的端元對應空間位置坐標Tab.1 Spatial coordinates of the extracted endmembers

從表1 可以看出，VCA和SMACC算法提取的端元中，端元1和端元5對應的空間位置坐標與模擬高光譜圖像中添加的異常點坐標相同，即VCA和SMACC算法中，兩個異常點被提取為端元光譜。而SPI+VCA和SPI+SMACC提取的端元中，沒有包含異常點或者其空間相鄰像元。這表明提出的預處理方法可有效避免異常點干擾。

圖6表明，SPI+VCA和SPI+SMACC算法提取的5個端元光譜曲線均與對應原始端元光譜曲線一致。而VCA和SMACC算法提取的端元中，只有端元2、端元3、端元4的光譜曲線與對應原始端元光譜曲線一致，端元1和端元5的光譜曲線與圖3中給出的異常點光譜曲線一致，這與表1的結果相符合。另外，由表1可知，SPI+VCA算法提取的端元2、端元4、端元5的空間坐標與SPI+SMACC算法提取的相應端元空間坐標相同，而圖6中兩種算法提取的端元光譜曲線之間卻有所差異。相比SPI+SMACC算法，SPI+VCA算法提取的端元光譜曲線與原始端元光譜曲線更加一致，且更加平滑。這是因為VCA算法在端元提取過程中，采用數據降維算法PCA去除了圖像中的噪聲，而SMACC算法沒有進行任何去噪處理。

為了定量評價不同算法的端元提取性能，采用SAD法計算了提取的各端元光譜與其對應的原始端元光譜之間的夾角，以衡量不同算法的端元提取精度，計算結果如表2所示。

表2 提取的端元與對應原始端元之間光譜夾角Tab.2 Spectral angle between the extracted endmembers and original endmembers

由表2可以看出，SPI+VCA和SPI+SMACC算法提取的5個端元與其對應的原始端元之間光譜夾角均較小，且相比SPI+SMACC，SPI+VCA對應光譜夾角平均小1.23°，這與圖6相符合；而VCA和SMACC算法提取的5個端元中，端元2、端元3、端元4與其對應的原始端元之間光譜夾角較小，端元1和端元5由于受異常點干擾，導致對應光譜夾角較大，這與表1結果相符合。VCA算法對應的光譜夾角比SPI+VCA算法對應的光譜夾角平均大10.1231°，SMACC算法對應的光譜夾角比SPI+SMACC算法對應的光譜夾角平均大7.921 6°。這表明經過預處理后端元提取精度得到了顯著提高。

4.3 第3組實驗

為了充分驗證預處理方法的性能，針對3.2中真實高光譜圖像Jasper Ridge進行了仿真實驗。圖7是SPI+VCA算法提取的Jasper Ridge中4個端元光譜曲線，與相應典型地物光譜曲線基本一致。在SPI+VCA算法中，參與端元提取的候選端元數量為726，遠少于原始像元數量10 000。

圖7 Jasper Ridge的4個端元光譜曲線Fig.7 Spectral curves of Jasper Ridge’s four endmembers

然而，由于較難確定真實高光譜圖像中精確的原始端元光譜，第3組實驗對真實高光譜圖像進行了光譜解混合處理，并比較了不同算法求解的豐度系數圖，以驗證預處理方法的性能。與第2組實驗類似，第3組實驗中仍分別比較了SPI+VCA和VCA算法，SPI+SMACC和SMACC算法。其中，對于SPI+VCA、SPI+SMACC、VCA算法，采用全約束最小二乘法FCLS[17]求解各端元對應的豐度系數；而對于SMACC算法，算法本身在提取端元的同時求解各端元對應的豐度系數。實驗結果如圖8所示。

通過對比分析圖8(a)和圖8(b)可知，SPI+VCA算法對應的各端元豐度圖與真實地物分布是一致的，且相比VCA算法對應的各端元豐度圖，顯然更好地區別出了樹、水和土壤，更清楚地識別出了道路線條，即每個像元的端元屬性更加顯著，各端元對應地物的區域性更加明顯，空間連續性更好，這說明相比VCA算法，SPI+VCA算法提取的端元更加精確，且抗噪能力更強。

(a) SPI+VCA、VCA與原始端元(a) SPI+VCA vs.VCA vs.original endmember

通過對比分析圖8(c)和圖8(d)可知，SPI+SMACC算法對應的各端元豐度圖與真實地物分布也是一致的，且相比SMACC算法對應的各端元豐度圖，顯然更準確、更清楚地識別出了各端元對應地物，尤其是水和土壤。SMACC算法中未能識別出水，且把土壤識別為兩種地物。這說明相比SMACC算法，SPI+SMACC算法提取的端元更加精確，且抗噪能力更強。

圖8 豐度圖Fig.8 Abundance maps

5 結 論

本文針對高光譜端元提取，提出了一種空譜聯合的預處理方法，先利用文中定義的光譜純度指數SPI測量高光譜圖像中每個像元的光譜純度，并在此基礎上，利用SPI的空間分布去除光譜純度較低的冗余像元，從而形成精簡的候選端元集。經過預處理后，采用常用端元提取方法，在候選端元集范圍內提取端元。

實驗結果表明，采用提出的預處理方法后，對于模擬高光譜圖像，提取的端元與原始端元之間夾角平均減少了9.022 3°，候選端元數量少于原始像元數量的10%。即提出的預處理方法顯著提高了抗干擾性，有助于有效提高端元提取精度；且大幅減少了參與端元提取的像元數量，明顯降低了時間復雜度；另外，該預處理算法不需要設定任何與圖像相關的閾值，因而具有較高的算法魯棒性。