方凌江 粘永健 雷樹濤 倪志揚(yáng)

（濟(jì)南軍區(qū)聯(lián)勤部后勤信息中心 濟(jì)南 250022）

1 引言

高光譜遙感是指光譜分辨率高達(dá)10－2λ數(shù)量級(jí)的光譜成像技術(shù)，在地質(zhì)勘探、軍事偵察等領(lǐng)域得到了廣泛的應(yīng)用［1］。由于傳感器空間分辨率的限制，以及地物在地表空間分布的復(fù)雜性，使得所拍攝的圖像中存在大量的混合像元。端元提取就是從存在混合像元的高光譜圖像中提取出場(chǎng)景的基本組成成分，它是高光譜解混的重要步驟。當(dāng)物質(zhì)之間的混合方式屬于非致密混合時(shí)，光子在物質(zhì)內(nèi)部之間的相互作用可以忽略，可用線性混合模型來描述混合像元的形成機(jī)理［2］。近十年來，學(xué)者們?cè)诰€性混合模型的基礎(chǔ)上提出了多種端元提取算法，例如純像元指數(shù)（Pure Pixel Index，PPI）［3］、頂點(diǎn)成分分析算法（Vertex Component Analysis，VCA）［4］、N－FINDR［5～6］，迭 代 誤 差 分 析 （Iterative Error Analysis，IEA）［7］，自動(dòng)形態(tài)學(xué)端元提取方法（Automated Morphological Endmember Extrac－tion，AMEE）［8］、最小體積單形體分析（Minimum Volumn Simplex Analysis，MVSA）［9］，其中，VCA算法是一種端元逐次提取算法，前面提取的端元不依賴后面提取的端元，有效避免了在端元數(shù)目估計(jì)不正確的情況下，所提取的端元之間相互干擾的問題，而且該算法提取速度快，計(jì)算復(fù)雜度低。但是VCA算法是將高光譜圖像中的極值像元作為端元，易受異常像元的干擾，而且它以高光譜圖像中的某個(gè)像元作為端元，受噪聲影響也較大。

本文在VCA算法的基礎(chǔ)上，利用地物在地表空間成片分布這一特點(diǎn)，采取空間鄰域迭代計(jì)算來提取空間純區(qū)域，并將這些區(qū)域內(nèi)的所有像元的均值作為待提取的端元，該算法可以有效地消除異常像元的干擾，并降低噪聲對(duì)端元的影響，快速有效地提取出端元。

2 算法描述

2.1 線性混合模型

在線性光譜混合模型中，每個(gè)像元的光譜響應(yīng)由各個(gè)端元的光譜響應(yīng)按照其對(duì)應(yīng)的豐度線性混合而成。用x∈RL×1表示波段數(shù)為L(zhǎng)的高光譜圖像中某像元的光譜響應(yīng)，x可以表示為

其中M表示端元數(shù)量，N 表示像元數(shù)量，E＝［e1，…，eM］∈RL×M表示構(gòu)成整幅高光譜圖像的端元矩陣，ei表示第i個(gè)端元，各個(gè)端元仿射獨(dú)立。s［n］＝［s1［n］，…，sM［n］］T∈RM×1表示第n個(gè)像元的豐度向量。w［n］∈RL×1為誤差向量，表示該像元被噪聲污染的程度。假定噪聲為零均值加性高斯噪聲，則w［n］～N（0，diag（σ2j）），其中diag（σ2j）為L(zhǎng)×L的對(duì)角矩陣，σ2j（j＝1，2，…，L）為各波段的噪聲方差。根據(jù)豐度的物理含義，豐度還要滿足非負(fù)約束（Abundance Non－negative Constraint，ANC）及和為 一 約 束 （Abundance Sum－to－one Constraint，ASC），可以表示為

在沒有噪聲的干擾下，所有滿足式（1）～（3）的像元構(gòu)成一個(gè)以各端元為頂點(diǎn)的凸面單形體，VCA正是利用正交投影的方法尋找單形體的頂點(diǎn)來提取端元，詳細(xì)步驟可以參考文獻(xiàn)［4］。此時(shí)，所有的純像元也都位于單形體的頂點(diǎn)，純像元就是端元，但是在有噪聲的干擾下，單形體體積膨脹，純像元偏離端元的位置，位于端元的周圍，此時(shí)端元不能由單個(gè)純像元表示。

用單形體來近似擬合有噪聲干擾的高光譜數(shù)據(jù)，靠近單形體的頂點(diǎn)的像元純度較高，PPI算法正是利用像元純度指數(shù)來提取端元的，可以認(rèn)為，由VCA算法所得到的像元是純像元，但是由于VCA所提取的純像元只有一個(gè)，受噪聲影響較大。對(duì)某種典型地物而言，通常它在高光譜圖像中分布有一定數(shù)量的純像元，如果可以找到更多該地物的純像元，并以這些純像元的均值作為該地物端元提取結(jié)果，根據(jù)統(tǒng)計(jì)學(xué)知識(shí)，此時(shí)的提取精度要比以單個(gè)像元作為提取結(jié)果的精度要高。

2.2 空譜聯(lián)合約束下的端元優(yōu)化

典型地物在空間上分布連續(xù)，純像元的鄰近像元仍然為該類地物的純像元的可能性較大?？兆V聯(lián)合約束要求VCA提取的純像元位于光譜相似、且空間連續(xù)的一致區(qū)域。下面詳細(xì)描述空譜聯(lián)合約束的實(shí)現(xiàn)過程：以VCA提取的單個(gè)純像元為種子，采用鄰近像元搜索，逐步擴(kuò)大純像元的數(shù)量，其過程如圖1所示。

圖1 鄰近像元搜索過程

記v為VCA在第i步所提取的純像元，搜索鄰近像元記為FindNV＿Pixels，實(shí)現(xiàn)的流程如表1所示。

表1 鄰像元搜索的流程圖

在表2中，V表示純像元集合，Nb表示取鄰域，本文中采用8鄰域搜索，mean為取均值，dist為衡量?jī)蓚€(gè)像元相似度的準(zhǔn)則，τ為門限。由FindNV＿Pixels得到的V為純像元區(qū)域，記V＝｛v1，…，vK｝，K 為V 中純像元的數(shù)量，V 具有兩個(gè)性質(zhì)：

1）?vk1，vk2∈V（k1≠k2），都有dist（vk1，vk2）≤τ成立，體現(xiàn)了對(duì)純像元區(qū)域具有光譜相似的約束。

2）?vk1，vk2∈V（k1≠k2），存在一條路徑U＝｛vk，vk∈V｝，使得vk1、vk2連通，體現(xiàn)了對(duì)地物在空間上具有連續(xù)性的約束。

然而，還有兩個(gè)尚待解決的問題：一是dist的選擇；二是τ的選擇。下面根據(jù)統(tǒng)計(jì)檢測(cè)理論中的Neyman－Pearson（NP）準(zhǔn)則來解決這兩個(gè)問題。vi、vj這兩個(gè)像元的差ε為

其中μ是由vi、vj兩個(gè)像元成分不同所引起的非噪聲差值。定義隨機(jī)變量r

由式（4）與式（5）可知，如果xj為該類地物的純像元，則ε～N（0，2diag（）），r為自由度為L(zhǎng) 的chi平方分布；否則，ε～N（μ，2diag（）），r為自由度為L(zhǎng)，非中心參量為‖μ‖的非中心chi平方分布。定義兩個(gè)假設(shè)：

為提高純像元區(qū)域的純度，則應(yīng)使vj誤判為純像元的概率最小，也就是使P（H0｜H1）最小。預(yù)設(shè)虛警概率PFA（PFA＝P（H1｜H0）），根據(jù)NP準(zhǔn)則得到：

其中τ為判決門限，可通過χ2L分布及PFA得到。與SPA、AMEE僅利用光譜信息來選取相似像元不同，本文算法利用了空間上下文信息。

在高光譜圖像中，通常會(huì)有大量的異常點(diǎn)存在，而降維及去噪等預(yù)處理措施并不能完全去除異常點(diǎn)，所以在降維及去噪后的高光譜圖像中仍然存在一定數(shù)量的異常點(diǎn)。由于異常點(diǎn)大多都屬于數(shù)據(jù)點(diǎn)云的極值點(diǎn)，VCA會(huì)受異常點(diǎn)的干擾，甚至在提取結(jié)果中直接將異常像元作為端元提取的結(jié)果，從而導(dǎo)致算法錯(cuò)誤。而異常點(diǎn)通常都是孤立點(diǎn)，具有數(shù)量少、在空間上不連續(xù)以及與周圍點(diǎn)的差異比較明顯等特征，可以利用這些特征來消除端元提取過程中異常點(diǎn)的影響。設(shè)定純像元區(qū)域大小門限為η，當(dāng)K≤η時(shí)就判定為異常像元。綜上，算法的整體流程為

表2 算法的整體流程圖

3 實(shí)驗(yàn)結(jié)果

將本文提出的算法記為SS－VCA，利用 Hysime［10］來估計(jì)噪聲，端元數(shù)目已知。與 AMEE、SSEE、SPEE等算法一樣，本文算法也利用了空間上下文信息，但是上述算法的性能依賴于多個(gè)參數(shù)的設(shè)置，而且參數(shù)的設(shè)置與場(chǎng)景的空間復(fù)雜程度有關(guān)，很難用統(tǒng)一的標(biāo)準(zhǔn)來比較本文算法與上述三種算法的優(yōu)劣，而且本文的重點(diǎn)在于強(qiáng)調(diào)對(duì)VCA的改進(jìn)作用，所以在實(shí)驗(yàn)中并不將上述三種算法作為比較對(duì)象。比較準(zhǔn)則為各個(gè)端元與真實(shí)端元之間的角距離，用ei、＾ei分別表示參考端元與算法提取的端元，其定義為

3.1 仿真數(shù)據(jù)實(shí)驗(yàn)

從美國(guó)地質(zhì)調(diào)查局（USGS）光譜庫(kù)中任意選擇五種地物的光譜作為仿真數(shù)據(jù)的端元，該光譜庫(kù)中地物光譜的波長(zhǎng)范圍是0.38μm～2.5μm，光譜分辨率為10nm，共有224個(gè)波段。豐度的生成方法是，先利用 MRF（Markov Random Field）生成空間分辨率為400×400的含有三種類別地物的分類圖，再由線性混合模型模仿由于傳感器分辨率不足而引起的地物混合現(xiàn)象，生成50×50的高光譜圖像。下面用兩組實(shí)驗(yàn)來測(cè)試算法的性能。在所有的試驗(yàn)中，虛警概率PFA設(shè)置為10－2，異常像元大小η為1。

圖2 不同信噪比條件下的算法性能比較

實(shí)驗(yàn)1（測(cè)試算法的抗噪能力）。在仿真數(shù)據(jù)中加入加性零均值高斯白噪聲，并逐漸改變信噪比（SNR）的大小，由10dB逐漸增加到30dB，步長(zhǎng)為5dB。圖2給出了不同算法在不同信噪比條件下的平均角距離，記為A＿SAD。從圖2可以得出，與VCA相比，SS－VCA具有更強(qiáng)的抗噪能力。

實(shí)驗(yàn)2（測(cè)試異常像元對(duì)算法的影響）。在USGS光譜庫(kù)中另外選擇一種地物，并在空間中隨機(jī)選擇五個(gè)不相鄰的位置，將原來數(shù)據(jù)中的像元的光譜替換為當(dāng)前地物的光譜，以模仿高光譜圖像中的異常像元。此時(shí)SNR設(shè)為10dB。由圖3可以直觀看出，VCA算法提取的端元中包含了異常像元，而SS－VCA則正常提取端元。

圖3 異常像元存在條件下的算法性能比較

3.2 真實(shí)數(shù)據(jù)實(shí)驗(yàn)

實(shí)驗(yàn)中采用了美國(guó)加州印度松樹測(cè)試地的高光譜數(shù)據(jù)，該數(shù)據(jù)由機(jī)載可見光及紅外成像光譜儀（Airborne Visible／Infrared Imaging Spectrometer，AVIRIS）拍攝而成，該數(shù)據(jù)波長(zhǎng)范圍是0.4μm～2.5μm，光譜分辨率為10nm，共有224個(gè)波段，空間分辨率為3.7m，截取大小為86×83子塊作為實(shí)驗(yàn)對(duì)象，該塊主要分布椰菜、玉米、雜草以及不同類別的生菜等自然地物，雖然該高光譜圖像空間分辨率較高，但由于地物分布復(fù)雜，仍存在相當(dāng)數(shù)量的混合像元。由于信噪比太低或者水吸收，36個(gè)波段被移除（波段號(hào)為1～4，104～112，150～167，220～224），僅保留188個(gè)波段用于數(shù)據(jù)處理。

為評(píng)價(jià)算法的性能，在實(shí)驗(yàn)中根據(jù)地表真實(shí)標(biāo)記圖手工從高光譜數(shù)據(jù)中選擇各類地物的多個(gè)純像元，并取各類純像元的均值作為參考光譜?？紤]到部分地物內(nèi)部存在較大的光譜變異性，難以用一個(gè)端元來描述，在實(shí)驗(yàn)中取端元數(shù)量為10。在提取端元之后，由全約束解混FCLS［11］估計(jì)的豐度，由SS－VCA提取的端元所得到的豐度估計(jì)如圖4所示。

比較真實(shí)標(biāo)記圖與估計(jì)豐度得到端元對(duì)應(yīng)的地物，其中生菜（5wk）對(duì)應(yīng)（d）、（f）、（h）三個(gè)端元，說明該地物具有較大的內(nèi)部變異性，而雜草對(duì)應(yīng)（g）、（j）兩個(gè)端元，分別表示開始衰老的雜草及綠色雜草。SS－VCA、VCA的端元提取結(jié)果分別與參考光譜的比較如表3所示。

圖4 Salinas－A數(shù)據(jù)及SS－VCA豐度所提取端元對(duì)應(yīng)的豐度

表3 Salinas－A數(shù)據(jù)的角距離比較（粗體數(shù)字表示最好的結(jié)果）

在克服異常像元的影響方面，在VCA提取的端元中，有一個(gè)端元為異常像元，而SS－VCA沒有。在提取結(jié)果的精度方面，由表3可以得出SS－VCA所提取的端元具有更小的平均角距離，精度更高。

4 結(jié)語(yǔ)

本文算法幾乎不受異常像元的影響，并具有較強(qiáng)的抗噪能力。但它只能找出其中一塊屬于某類地物的純區(qū)域，不能將空間不相鄰但屬于該類地物的所有純區(qū)域同時(shí)找出來，抗噪性能的提高還有一定的空間，這是下一步要研究的方向。

［1］Keshava N，Mustard J F.Spectral Unmixing［J］.IEEE Signal Processing Magazine，2002，19（1）：44－57.

［2］Adams J B，Gillespie A R.Spectral－Mixture Analysis［M］.Remote Sensing of Landscapes with Spectral Images：A physical Modeling Approach，Cambridge University，2004：126－167.

［3］Boardman J W，Kruse F A，Green R O.Mapping Target Signatures via Partial Unmixing of AVIRIS Data［C］.1995.

［4］Nascimento J M P，Bioucas－Dias J M.Vertex Component Analysis A Fast Algorithm to Unmix Hyperspectral Data［J］.IEEE Transaction on Geoscience and Remote Sensing，2004，43（8）：898－910.

［5］Winter M E.N－FINDR：An Algorithm for Fast Autonomous Spectral Endmember Determination in Hyperspectral Data［C］／／Proceedings of SPIE，Denver，1999，3753：266－275.

［6］Plaza A，Chang C I.An Improved N－FINDR Algorithm in Implementation.Algorithms and Technologies for Multispectral，Hyperspectral and Ultraspectral Imagery［C］／／Proceedings of SPIE，Bellingham，WA，2005，5806：298－306.

［7］Neville R A，Staenz K，Szeredi T，et al.Automatic Endmember Extraction from Hyperspectral Data for Mineral Exploration［C］／／Proceedings of International Conference on Airborne Remote Sensing，Ottawa，Canada，1999：21－24.

［8］Plaza A，Martinez P，Perez R，et al.Spatial／Spectral Endmember Extraction by Multidimensional Morphological Operations［J］.IEEE Transaction on Geoscience and Remote Sensing，2002，40（9）：2025－2041.

［9］Li J，Bioucas－Dias J E M.Minimum Volume Simplex Analysis：a Fast Algorithm to Unmix Hyperspectral Data［C］／／Proceedings of IEEE International Conference on Geoscience and Remote Sensing，Boston，USA，2008：7－11.

［10］Bioucas－Dias J M，Nascimento J M P.Hyperspectral Subspace Identification［J］.IEEE Transaction on Geoscience and Remote Sensing，2008，46（8）：2435－2445.

［11］Heinz D C，Chang C.Fully Constrained Least Squares Linear Spectral Mixture Analysis Method for Material Quantification in Hyperspectral Imagery［J］.IEEE Transaction on Geoscience Remote Sensing，2001，39（3）：529－545.