[發明專利]一種基于隨機投影的高光譜圖像稀疏解混方法

有效

• 申請號: 201110207433.0
• 申請日: 2011-07-23
• 公開（公告）號: CN102314685A
• 公開（公告）日: 2012-01-11
• 發明（設計）人: 史振威;翟新雅;都仁扎那
• 申請（專利權）人: 北京航空航天大學
• 主分類號: G06T7/00
• 分類號: G06T7/00
• 代理公司: 北京慧泉知識產權代理有限公司 11232
• 代理人: 王順榮;唐愛華

1.一種基于隨機投影的高光譜圖像稀疏解混方法，其特征在于：該方法具體步驟如下：

步驟一：用計算機讀取數據：計算機在MATLAB?R2008b環境下讀取高光譜圖像數據，數據來源于成像光譜儀采集到的遙感圖像，得到數據立方體，高光譜圖像數據去除被水汽吸收的波段和信噪比較低的波段，將高光譜圖像數據逐像素點排列，得到原始的高光譜圖像矩陣V；設高光譜圖像有m個波段，共有n個像素點，V＝[v₁，v₂，...，v_n]，v_i是高光譜圖像第i個像素點的光譜列向量，是一個m維向量；讀取現有的高光譜數據庫數據，選擇光譜庫中純物質光譜數據構建光譜庫矩陣S；設純物質光譜的數量為q，S＝[s₁，s₂，...，s_q]，s_i是光譜庫中第i個純物質的光譜列向量，同樣是一個m維向量；

步驟二：高光譜數據隨機投影：用計算機讀取數據后，對高光譜圖像數據和高光譜庫數據進行隨機投影，就是用一個所包含的元素符合正態分布的矩陣R_d與矩陣V，S分別相乘；設矩陣R_d的大小為d×m，d＜m，得到降維后的數據R_dV和R_dS，其矩陣大小分別為d×n，d×q，代替原始的高光譜圖像數據V和高光譜數據庫S參與下述步驟三的計算，其中d是隨機投影后光譜波段的個數；

步驟三：高光譜圖像稀疏解混：數據進行隨機投影后，就對降維后的高光譜圖像數據進行光譜解混；設高光譜數據庫中含有足夠豐富的純物質光譜特征，那么高光譜圖像中含有的端元只占極少一部分，也就是說高光譜圖像中每個像素點的光譜曲線是由高光譜數據中少數的純物質的光譜曲線線性組合構成，這體現了端元的稀疏性表達；端元在高光譜數據庫中是稀疏的，他們所對應的豐度也具有稀疏性，即豐度矩陣是稀疏的，以豐度的稀疏性作為正則化項，符合實際的物理意義；在這里，采用線性混合模型，即將一條檢測得到的光譜曲線，分解為光譜庫中純物質光譜線性組合的形式，其系數即為相應豐度；設豐度矩陣為H，其大小為q×n，其所有元素滿足非負性，線性混合模型如下表示：

V＝SH

對原始高光譜圖像稀疏解混的目標函數為：

minH||H||1]]>????s.t.V＝SH，H≥0

||H||₁表示H包含的所有元素絕對值的和，因為豐度具有非負性，||H||₁越小，豐度矩陣H含有的非零元素值將越少，其表現出的稀疏性就較強；將上式轉變為無約束優化問題：

minH12||V-SH||22+λ||H||1]]>

在上式目標函數中，l₂項是擬合項，體現的是解混的結果與原高光譜圖像的擬合程度；l₁即稀疏性正則化項；λ是正則化參數，調節著擬合項與正則化項的權重，其大小根據實際高光譜數據設置；因為V，S數據量很大，在迭代求解的過程中占用內存較大，并且計算時間較長；采用隨機投影算法，能夠有效保持上式中的l₂距離項，實現數據量壓縮，隨機投影后的目標函數為：

minH12||RdV-RdSH||22+λ||H||1]]>

針對上述目標函數，采用分裂Bregman算法求解上述函數極值，從而得到該優化模型下的最優解——最優豐度矩陣H；

引入輔助變量D，并且D＝H，上式目標函數轉化為：

minH,d12||RdV-RdSH||22+λ||D||1s.t.D=H]]>

轉化無約束形式為：

minH,dλ||D||1+μ2||RdV-RdSH||22+12||D-H||22]]>

其中μ是正則化參數，其大小根據實際高光譜數據設置。，引入誤差反饋變量B，上述優化問題轉變為迭代求解過程，如下：

Hk+1=argmin{μ2||RdV-RdSH||22+12||Dk-H-Bk||22}Dk+1=argmin{λ||D||1+12||D-Hk+1-Bk||22}]]>

具體求解步驟如下：

(1)初始化豐度矩陣H以及輔助變量矩陣B，D，V_d，H⁰＝B⁰＝D⁰＝0；V_d⁰＝R_dV

(2)迭代變量H，D，B，V_d，交替迭代過程如下：

Hk+1=[μ(RdS)T(RdS)+I]-1[μ(RdS)TVdk-Bk+Dk]Dk+1=max(Hk+1+Bk-λ,0)Bk+1=Bk+Hk+1-Dk+1Vdk+1=Vdk+RdV-RdSHk+1]]>

這里k為迭代次數，V_d^k+1的迭代過程可以加速收斂；

(3)停止迭代條件：當豐度矩陣H相鄰兩次迭代的相對誤差不大時，迭代停止；迭代停止條件為：

||Hk+1-Hk||22||Hk||22<tol]]>

這里，取tol＝10^-6，若迭代步數少于10次就能收斂，使用迭代次數控制迭代過程，設迭代次數為1000步；

步驟四：獲取豐度圖及真正的端元：在獲得豐度矩陣后，設置合適的閾值，豐度矩陣中小于閾值的元素設為零，不小于閾值的元素不做處理；找出豐度矩陣中含有非零元素的行標號，豐度矩陣的行標號與光譜庫的列標號對應；取出光譜矩陣中列標號所對應的列，即為真正的端元光譜；取出豐度矩陣中行標號所對應的行，即為端元所對應的豐度圖；在豐度圖中顏色較亮的區域表示端元所占的比例較大，較暗的區域對應的比例小，顏色為黑色的區域表示不含有該端元。