同質區共享端元變異性的高光譜混合像元分解

王寧，保文星*，屈克文*，馮偉

（1.北方民族大學 計算機科學與工程學院，寧夏 銀川 750021；2.西安電子科技大學 電子工程學院，陜西 西安 710071）

1 引言

高光譜遙感圖像包含著豐富的光譜信息和空間信息，已經廣泛應用于資源勘探、災害預警、國家安全、城市建設等諸多領域［1］。因為高光譜圖像的空間分辨率一般較低，而地物的分布比較復雜，圖像中一個像元體現出來多種地物的混合特征，這阻礙了以像元為單位的圖像處理算法的發展［2］。高光譜圖像混合像元分解也稱高光譜解混（Hyperspectral Unmixing，HU），它根據混合的光譜特征去反演真實的地物分布比例，對地物的空間分布進行更精細的解譯［3］?；旌舷裨纸鈱⒐庾V圖像分解為端元和與端元對應的豐度，端元為構成高光譜圖像的基本地物，豐度為地物在某一混合區域中所占的比例［4］。端元提取算法主要包括：純像元指數（Pixel Purity Index，PPI）［5］、內部最大體積法［6］和頂點成分分析（Vertex Component Analysis，VCA）［7］等。豐度估計的經典算法為同時考慮豐度非負并且和為一約束的全約束的最小二乘方法（Fully Constrained Least Squares，FCLS）［8］。豐度的稀疏性是豐度常見的特性，它說明高光譜圖像每個像元中只含有少量的地物［9］。

端元的光譜由于照明條件或大氣條件的差異及其他因素導致相同地物的光譜出現差異，出現“同物異譜”的現象，導致端元光譜出現變異性［10］。線性混合模型（Linear Mixing Model，LMM）中的一個隱含的假設是在一幅圖像中每一種材料通常由一個單獨的光譜特征表示，認為總的光譜曲線是各個地物的光譜曲線按照其占比線性疊加［11］。從物理意義上來說，不能用單一的端元光譜曲線來進行混合像元分解，這也是LMM 的一個缺陷［12］。由于大部分解混算法都是迭代求解端元和豐度，如果忽略端元的光譜變異性，會導致在解混的過程中傳播錯誤的建模誤差，從而影響解混結果［13］。

相同地物選取多個端元光譜參與到混合像元分解中，可以有效處理端元的變異性。使用多端元解混中經典的算法有多端元光譜混合分析法（Multiple Endmember Spectral Mixture Algorithm，MESMA）［14］、光譜束解混法［15］及局部解混法［16］，這些算法本質上在圖像的不同空間位置上允許選用不同的端元進行解混。而基于統計模型的方法利用先驗概率知識求解端元的分布，主要包括基于正態模型的（Normal Compositional Model，NCM）方法［17］、基于Bayesian 的方法［18］和基于β模型的（Bayesian Compositional Model，BCM）方法［19］。

在LMM 的基礎上擴展的模型可以處理解混過程中端元變異性帶來的不利影響，這些模型主要包括：擾動線性混合模型（Perturbed Linear Mixing Model，PLMM）［20］、擴展線性混合模型（Extended Linear Mixing Model，ELMM）［21］、廣義線性混合模型（Generalized Linear Mixing Model，GLMM）［22］、增廣線性混合模型（Augmented Linear Mixing Model，ALMM）［23］和縮放擾動線性混合模型（Scaled and Perturbed Linear Mixing Model，SPLMM）［24］。PLMM 改進傳統的LMM，添加空間變化的擾動端元來描述端元的光譜變異性。ELMM 在解混中考慮空間信息和端元光譜變異性的同時，引入的尺度因子可以有效處理由光照條件引起的端元的光譜變異性。GLMM 認為不同波長可能會受到不均勻的影響，對每一個像元都引入一個與端元矩陣相同維度的矩陣與端元矩陣進行笛卡爾積運算，對端元施加更加精細的放縮。ALMM 使用共享的尺度因子，在一個像元中使用相同的尺度因子對所有的端元進行放縮，使用字典學習的策略構建尺度因子之外的光譜變異性，SPLMM 在PLMM 中添加共享尺度因子，一個像素中的所有端元共享相同的尺度因子。

在深度學習領域，一些處理端元光譜變異性的算法也取得較好的效果，經典算法主要包括：深度生成端元模型（Deep Generative Unmixing，DeepGUn）［25］、概率生成模型（Probabilistic Generative Model for Spectral Unmixing，PGMSU）［26］、空間光譜解混深度生成模型（Deep Generative Model for Spatial-Spectral Unmixing，DGMSSU）［27］和三維卷積神經網絡（3D-Convolutional Neural Network，3D-CNN）［28］等。上述算法主要利用神經網絡的非線性建模能力，聯合空間信息與光譜信息，更好地擬合端元的分布情況。

照明條件的差異是引起端元變異性的一個重要因素，而尺度因子對光譜曲線的縮放能力已被證實可以有效處理由照明條件不同導致的端元變異性［13］。此外，高光譜圖像中存在著大量的局部同質區，一個像元與其周圍的像元有著相似的地物分布，端元的變異性也相對較小。針對PLMM 難以有效處理縮放效應造成的端元變異性這個問題，本文對PLMM 進行改進，提出了基于局部同質區共享端元變異性的解混算法（Shared Endmember Variability in Unmixing，SEVU），具體工作內容如下：

（1）在PLMM 中添加尺度因子以提高模型處理端元等比例變異性的能力，并提出一個針對尺度因子的正則項以約束端元縮放效應的強度；

（2）在解混的過程中引入超像素分割算法劃分局部同質區域，在局部同質區內各個像元共享相同的變異端元，將端元變異性的尺度由像元級轉變為超像素級；

（3）設計了基于交替方向乘子法（Alternating Direction Method of Multipliers，ADMM）的解混算法以求解每個局部同質區的變異端元，并在一組合成數據集和兩個真實數據集上進行實驗，驗證了算法的有效性。

2 高光譜解混方法

2.1 線性混合模型

2.1 線性混合模型

由于模型簡單，物理意義明確，LMM是應用最為廣泛的模型，它將高光譜圖像看作由標準端元光譜矩陣與豐度比例矩陣的乘積。LMM中總的光譜曲線是由各個地物的光譜按照其占比線性疊加得到的。對于LMM，通常采用兩個約束：豐度非負約束（Abundance Nonnegativity Constrained，ANC）和豐度和為一約束（Abundance Sum-to-one Constrained，ASC）。

LMM的向量形式為如公式（1）所示。在LMM中，Y=[y1|y2|…|yn]∈RL×N表示整幅圖像，yn∈RL×1表示圖像中的第n個像元，M∈RL×K表示端元，an∈RK×1表示yn對應豐度，rn∈RL×1表示第n個像元的重構誤差，R=[r1|r2|…|rn]∈RL×N表示整幅圖像的重構誤差。

其中：L表示光譜波段數目，K表示端元數目，N表示這幅圖像包含的像元個數。

2.2 擾動線性混合模型

擾動線性混合模型（PLMM）根據提取的端元可以被認為是參考端元的變化實例，改進LMM，添加逐像元的空間變化的擾動因子來描述端元的光譜變異性。PLMM如公式（2）所示，dMn∈RL×K表示yn對應的擾動端元矩陣。擾動端元的具體作用如圖1（a）所示：

圖1 不同因素對端元的影響Fig.1 Effect of different elements on endmembers

2.3 擴展線性混合模型

擴展線性混合模型（ELMM）在解混時用尺度因子來表示各種地物端元在所有像元中的變化，每種地物變異后的端元為對應參考端元矩陣乘一個尺度系數，這些尺度系數構成尺度因子矩陣。在ELMM中，Sn∈RK×K表示yn對應的尺度因子矩陣，是一個非負的對角矩陣。尺度因子可以對光譜曲線進行等比例放縮，它的具體作用如圖1（b）所示：

3 改進算法

3.1 添加尺度因子的擾動線性混合模型

改進后的模型如公式（4）所示，Sn和dMn同時構造端元的變異性，Sn構造端元的等比例變異，而dMn構造端元的非等比例變異，其中第n個像元處的變異端元為MSn+dMn，尺度因子S和擾動端元dM對光譜曲線的影響如圖1（c）所示。

端元矩陣和變異后的端元矩陣是非負的，需要對其進行非負約束，本算法也考慮豐度的物理意義，使用豐度非負約束（ANC）和豐度和為一約束（ASC）。綜上，SEVU 算法的約束條件為：

公式（6）中α，β，γ，η是四個正則化參數，式（7）～式（10）分別為關于豐度、端元、擾動端元和尺度因子的四個正則化函數。

大多數混合像元中往往含有少量的地物種類，故針對豐度施加稀疏約束，通過公式（7）可以控制豐度的稀疏性：

因為目標函數的求解較為復雜，需要對多個變量進行迭代求解，有必要對端元進行合適約束以減小算法的解空間，所以使用端元之間的相互距離進行約束，如公式（8）所示：

其中：ek表示RK中正交基中的第k個向量，M=表示所有端元第l波段上的行向量，IK為維度為K的單位矩陣。

端元光譜變異性是有限的，變異后的光譜與原始光譜差距不會太大，為了使得模型具有更好的性能，有必要對變異性能量進行約束，其中對擾動端元的懲罰為：

當變異后光譜與原始光譜差值越小，端元變異性能量應當越小，受公式（9）的啟發，對于尺度因子矩陣，本文提出針對等比例放縮的變異能量的懲罰項為：

3.2 超像素分割

PLMM 和ELMM 中端元變異性都是逐像元構建的，即一個像元對應著一個變異的端元，考慮到局部同質區中端元的變異性較小，所以在SEVU 算法中，一個局部同質區內多個相似的像元共享一個變異的端元。SEVU 算法使用超像素分割（Simple Linear Iterative Clustering，SLIC）［29］劃分局部同質區。具體的像元間的相似度如公式（11）所示：

其中：ix，cx，iy，cy為兩個像元的空間位置，Di，c是總的相似度，λ用于平衡光譜距離和空間距離的權重。其中SAM和SID計算如公式（12）和公式（13）所示，pil表示像元i在第l波段的輻射值：

對于劃分出的任意一個超像素yslic在迭代求解時看作最小的單位，yslic=其中num為此超像素中包含的像元個數。所以對于整幅圖像Y=[y1|y2|…|yn]∈RL×N，進行超像素變換，可得到變換后的圖像YSLIC=[y1|y2|…|yq]∈RL×q，其中q為超像素的個數，因為q?N，超像素變換可以看作高光譜圖像的空間域的降維操作。同理，對于豐度矩陣A=[a1|a2|…|an]∈RK×N，按照上述變換規則，得到Aslic∈RK×q。

3.3 SEVU 算法

SEVU 算法在適度的局部同質區內共享一致的端元變異性，一塊超像素代表一塊局部同質區。給定一幅圖像，超像素的塊數影響同質區的大小，適度則表明這幅圖像中超像素的大小并非固定的。共享端元變異性指在每一塊超像素（同質區）內，所有像元使用相同的由尺度因子和擾動構造的變異端元求解最終的豐度，即SEVU 算法只構建同質區間的端元變異性，忽略同質區內的端元變異性。圖2 為SEVU 算法示意圖，首先輸入高光譜圖像，然后先由VCA 算法和FCLS算法對端元和豐度進行初始化，然后利用空間上下文信息將高光譜圖像劃分為若干個超像素，然后將超像素作為基本單位使用ADMM 求解變異后的各個局部同質區的端元，最后使用FCLS 算法并利用像元與超像素的隸屬關系求解每個像元的豐度。

圖2 SEVU 算法示意圖Fig.2 Diagram of SEVU algorithm

3.4 算法求解

算法采用ADMM 進行求解。對于本模型中各個變量的求解，等價于最小化SEVU 算法的目標函數f(A，M，dM，S)，使用ADMM 將目標函數的求解劃分為四個子問題：豐度矩陣A，端元矩陣M，擾動端元矩陣dM，尺度因子矩陣S，然后采用拉格朗日乘子法（Lagrange Multiplier）求解各個子問題。

3.4.1 豐度

然后得到關于an的增廣拉格朗日形式的函數，如公式（15）所示：

3.4.2 端元

P=[P1|P2|…|Pn]，其中Pn表示按照如下規則將Sn中的對角元素轉化為一個列向量：

3.4.3 擾動端元

3.4.4 尺度因子

對Sn求偏導數，令偏導數等于0，得到以Sn為自變量的西爾維斯特方程，如公式（31）所示，可以使用Matlab 中Sylvester 函數直接求解：

3.4.5 算法流程

SEVU 算法的流程如算法1 所示，首先，局部同質區的劃分是依賴超像素分割算法SLIC 實現，超像素分割的塊數q控制超像素的大小。然后利用ADMM 迭代求解每塊超像素的變異端元，最后使用已求解的變異端元遍歷超像素內的所有像元，使用FCLS 算法求得最終的豐度，其中停止條件為達到最大循環次數或目標函數接近收斂

4 實驗結果與分析

4.1 實驗數據集及評價指標

實驗采用開源代碼生成高光譜合成數據集（空間維度：50×50，光譜維度：198）［13］，此數據集包含三個端元：植被（Vegetation），泥土（Soil）和水（Water）。端元光譜采用Simplified Illumination Model，Prospected-D Model，Hapke Model和Simple Atmospheric Compensation Model 來構建，豐度圖像使用高斯隨機場生成，豐度采用高斯隨機場生成，然后基于LMM 逐像素生成高光譜圖像，最后再將合成的圖像添加20～40 dB 的高斯噪聲。真實數據集使用解混中常用的數據集：Jasper Ridge 數據集（空間維度：100×100，光譜維度：198）和Cuprite 數據集（空間維度：150×270，光譜維度：188）。三個實驗數據集的偽彩色圖像及其端元光譜曲線如圖3 所示（彩圖見期刊電子版）。

圖3 數據集Fig.3 Datasets

評價指標包括平均端元光譜角距離（mean Spectral Angle Distance，mSAD）、豐度均方根誤差（abundance Root Mean Square Error，aRMSE）和重構均方根誤差（reconstruction overall Root Mean Square Error，rRMSE），以上三個評價指標均越低越好。需要說明，在端元變異場景下，需要估計每個像元的端元，但在真實數據集中端元的變異性是未知的，僅有一個可供參考的端元矩陣，所以計算真實數據集的端元提取精度時使用算法求解出的多個端元矩陣的平均值與真實數據集的端元矩陣計算mSAD［27］。

4.2 實驗

為了驗證SEVU 算法的有效性，實驗使用的數據包括一組合成數據集和兩個真實數據集，對比算法包括FCLS，SCLSU，PLMM，ELMM，GLMM，ALMM，所有實驗均使用VCA 進行端元初始化。為了保證實驗的公平性，在實驗時已經對PLMM，ELMM，GLMM，ALMM 中的參數進行調優處理，均使用FCLS 進行豐度初始化。在接下來的所有實驗中，高光譜圖像中端元的數目已知，最終實驗結果為10 次實驗的平均值，最好的結果以粗體顯示。

4.2.1 魯棒性實驗

為了驗證算法的魯棒性，實驗一為不同噪聲（20～40 dB）的場景下所有實驗算法的對比實驗。通過圖4 可以看出，SEVU 算法在mSAD 和aRMSE 取得較好的效果。在端元估計中，ELMM 表現的效果最差，其次是GLMM，而添加擾動端元的PLMM 取得較好的效果，而FCLSU，SCLSU 和ALMM 沒有求解端元；在豐度估計中，基于尺度因子的ELMM 和SCLSU 表現較差，而FCLSU 的效果較好；對于yRMSE，基于字典學習的ALMM 和使用更多尺度因子建模的GLMM 效果較好。隨著信噪比的增加，所有算法的性能逐漸變強，由于超像素有降噪的作用，SEVU 算法在信噪比較低的情形下具有優秀的解混性能。

圖4 信噪比對算法性能的影響Fig.4 Effect of signal-to-noise ratio on algorithm performance

圖5 參數對SEVU 算法性能的影響Fig.5 Effect of parameters on SEVU algorithm performance

4.2.2 參數實驗

為了確定模型中的參數選擇對混合像元分解精度的影響，實驗二在Jasper Ridge 數據集上驗證四個正則化參數的值對算法性能的影響，所有參數選擇的數量級從10-3到101。圖 5為α，β，γ，η的值對mSAD 和aRMSE的影響，mSAD 和aRMSE 對參數設定比較敏感，所以參數實驗時主要考慮mSAD 和aRMSE。γ的值對算法性能最敏感，針對此數據集，γ=1，其余的參數等于0.1 時SEVU 算法性能較好。

4.2.3 超像素實驗

在SEVU 算法中，超像素的塊數是一個變量，代表著局部同質區劃分的大小，同一幅圖像，分割的塊數越多，超像素包含的像元越少，表明同質區的同質性越強。超像素分割的塊數影響SEVU 算法的精度和運行時間，實驗數據集為：Jasper Ridge 數據集，實驗時控制其他變量為：[α=0.1，β=0.1，γ=1，η=0.1]。圖6 為此數據集分割塊數后的平均光譜圖像，結合表1 可以看出，yRMSE 始終與超像素塊數成反比；算法運行時間與超像素塊數呈正比。超像素分割的塊數太少，從圖像中使用的特征也少，求解的端元誤差較大，故SEVU 算法的準確度很差；分割的塊數越多，會忽視像元間的空間信息，難以取得較好的效果。當超像素分割塊數為200～600 時，合理利用圖像的空間信息，分割塊數適中，SEVU 算法的性能較好。

表1 超像素塊數對mSAD，aRMSE，yRMSE 的影響Tab.1 Effect of number of superpixel blocks on mSAD，aRMSE，yRMSE

圖6 超像素分割后圖像的塊數Fig.6 Number of blocks in the image after superpixel segmentation

4.2.4 合成數據集實驗

在合成數據集（30 dB）上超像素分割的塊數為624，SEVU 算法的各個參數設定為：[α=1，β=0.01，γ=1，η=0.1]。通過表2可以看出，SEVU 算法在端元估計時準確性最高，對三種端元在平均光譜角距離上優于其他算法，其次是PLMM，GLMM，ELMM?？偟膩碚f，各個算法對水（Water）端元的處理能力最弱，說明考慮光譜變異性的解混算法均無法有效處理水這種弱信號端元的變異性。圖7 顯示了SEVU 構造的三個端元的光譜束，通過圖8 和表3 可以看出，SEVU 求得的豐度更接近真實豐度，aRMSE更低。

表2 各個算法在mSAD 上的評價結果（合成數據集SNR_30 dB）Tab.2 Evaluation results of each algorithm for mSAD（Synthetic dataset SNR_30 dB）

表3 各個算法在aRMSE 和yRMSE 上的評價結果（合成數據集 SNR_30 dB）Tab.3 Evaluation results of each algorithm for aRMSE and yRMSE（Synthetic dataset SNR_30 dB）

圖7 SEVU 算法構造出的端元束（合成數據集 SNR_30 dB）Fig.7 Endmember bundles constructed by SEVU algorithm（Synthetic dataset SNR_30 dB）

圖8 各個算法估計的豐度圖像（合成數據集 SNR_30 dB）Fig.8 Abundance estimated by each algorithm（Synthetic dataset SNR_30 dB）

4.2.5 真實數據集實驗

對于Jasper Ridge 數據集，超像素分割塊數：601塊，參數設定為[α=0.1，β=0.1，γ=1，η=0.1]。圖9 顯示了SEVU 算法估計的四組相似的光譜束。通過表4 可以看出，除水之外其余端元的求解，SEVU 算法在mSAD 上均取得最好的效果；而對于水這種端元的求解，單獨使用擾動端元的PLMM 效果最好，也反映擾動端元處理弱信號的優越性。在此數據集上，考慮端元變異性的算法求解的端元均比VCA 效果好，也說明在真實場景下考慮端元變異性的必要性。

表4 各個算法在mSAD 上的評價結果（Jasper Ridge）Tab.4 Evaluation results of each algorithm for mSAD（Jasper Ridge）

圖9 SEVU 算法構造出的端元束（Jasper Ridge）Fig.9 Endmember bundles constructed by SEVU algorithm（Jasper Ridge）

結合圖10 和表5 可以看出，SEVU 算法總體的豐度估計也是最優的，尤其對道路（Road）這種端元的求解的誤差最低，進而對道路豐度的估計最接近真實值，說明在自然場景下基于局部同質區共享端元變異性兩階段分別求解端元塊和豐度的可行性和有效性。

表5 各個算法在aRMSE 和yRMSE 上的評價結果（Jasper Ridge）Tab.5 Evaluation results of each algorithm for aRMSE and yRMSE（Jasper Ridge）

圖10 各個算法估計的豐度（Jasper Ridge）Fig.10 Abundance estimated by each algorithm（Jasper Ridge）

對于Cuprite 數據集，超像素分割塊數：10 124 塊，參數設置：[α=0.01，β=1，γ=10，η=0.1]。圖11 可以看出不同端元變異的強度是不同的，變異性弱的端元束幾乎是單一端元的重復，表現出單一曲線的形態。本數據集極為復雜，因此超像素分割的塊數較多，平均5 個像元構成一個超像素，對于本數據集，超像素主要起到降維的作用。觀察表6 和表7 可以看出，SEVU 算法求解的mSAD 更加準確。所有算法在本數據集上的mSAD 較大，但各個算法之間的差異較小，大致在0.10～0.12 之間，也可以通過圖12 看出各個算法求得的豐度圖相似性較高。

表6 各個算法在yRMSE 上的評價結果（Cuprite）Tab.6 Evaluation results of each algorithm for yRMSE（Cuprite）

表7 各個算法在mSAD 上的評價結果（Cuprite）Tab.7 Evaluation results of each algorithm for mSAD（Cuprite）

圖11 SEVU 算法構造出的部分端元束（Cuprite）Fig.11 Partial endmember bundles constructed by SEVU algorithm（Cuprite）

圖12 各個算法估計的部分豐度（Cuprite）Fig.12 Partial abundance estimated by each algorithm（Cuprite）

4.2.6 算法收斂性分析

大多數基于LMM 的變體算法的求解都是非凸的，因此需要大量的先驗信息或施加較多的正則項來限制解空間的區域，使算法的優化盡可能收斂。目標函數值隨算法迭代次數的曲線可以從主觀上判斷算法的收斂情況，通過圖13 可以看出對于不同噪聲場景下的合成數據集和真實數據集，SEVU 算法的收斂性較好。

圖13 目標函數隨算法迭代次數的曲線Fig.13 Profiles of the objective function with the number of iterations of the algorithm

圖14 算法運行時間對比（s）Fig.14 Comparison of algorithm run times（s）

4.2.7 算法運行時間

盡管SEVU 算法同時使用尺度因子和擾動構建端元的變異性，由于使用超像素分割算法作為空間上的預處理來劃分局部同質區，而求解超像素內各像元的豐度向量使用豐度全約束（ASC，ANC）的FCLS 算法效率較高，這極大降低SEVU 算法運行的復雜度。通過圖 14 可以看出，SEVU 算法的運行時間在圖像較為復雜的時候有較大優勢，而真實的場景中圖像往往比較龐大且復雜，說明SEVU 算法在真實場景中應用的可行性。

5 結論

本文改進了PLMM，引入尺度因子的同時并利用超像素分割算法劃分局部同質區，最后設計出一個基于局部同質區內共享端元變異性的SEVU 算法。在合成數據集和真實數據集上進行一系列實驗，相較于PLMM，SEVU 算法在各個數據集上的端元誤差降低了1.78%～3.23%，豐度誤差降低了0.26%～2.42%，并且在大多數場景下SEVU 算法的解混效果優于對比算法，實驗結果驗證SEVU 算法的有效性和可行性。