本發(fā)明屬于圖像處理技術領域，特別涉及極化sar圖像的變化檢測方法，可應用于城市規(guī)劃、生態(tài)環(huán)境考察以及自然災害的評測。

背景技術：

極化合成孔徑雷達polsar圖像變化檢測是一種從不同時間同一地點的兩幅極化sar圖像中提取變化信息，生成差異圖，確定地物變化信息的遙感圖像處理技術。極化sar能接收四個通道的回波信號，它能夠更加全面地表示目標的散射機理，所以其圖像所包含的信息量遠大于單通道sar圖像。近年來，極化sar圖像變化檢測已成為圖像處理研究的一個新研究方向，被廣泛的應用于各個領域，比如土地覆蓋和利用的變化監(jiān)測、城市規(guī)劃、環(huán)境監(jiān)測分析、自然災害估計等。

目前，極化sar圖像變化檢測的研究還處于初步階段，大致分為三類。

第一類是以強度信息為特征，將傳統(tǒng)變化檢測技術提取極化sar影像的差異信息，如將主成分分析、閾值分割及矩陣分解的等方法。但這類方法的不足之處是不能充分利用極化sar的散射信息。

第二類是以極化sar數(shù)據(jù)為基礎，利用極化sar統(tǒng)計分布來提取極化sar的差異信息，如基于極化協(xié)方差矩陣的似然比檢驗的變化檢測方法。不過，該方法的應用前提是地物目標的極化協(xié)方差矩陣滿足wishart分布，但實際地物散射特性比較復雜，有時很難滿足此條件，因此方法的通用受到限制。

第三類是利用一些極化目標分解模型對極化sar數(shù)據(jù)進行特征提取，然后尋找合適方法提取差異信息。極化目標分解是通過若干已知的基本散射解譯目標散射。典型的方法有pauli分解和cloude分解，其中pauli分解將將目標散射機制分解為三個基本散射矩陣的線性組合；cloude分解將極化目標分解為三個基本散射分量的組合。然后再進行差異信息提取，如差值、比值、小波融合等，最后對差異信息進行優(yōu)化，得到變化檢測結果。由于極化sar包含信息的多樣性、和地物散射特性的復雜性。

目前，不少的研究是沿著第三類方法的方向進行的，而這類方法由于在利用極化目標分解模型時，需要在分解前確定極化目標是否相干，以便得到好的極化特征，然而，在分解前是不易判斷極化目標是否相干；同時由于該類方法在提取差異圖時，沒有利用圖像的空間信息，也未考慮噪聲的影響，從而使檢測結果不準確。

技術實現(xiàn)要素：

本發(fā)明的目的在于針對上述已有技術的不足，提出一種基于散射功率特征與低秩稀疏模型的極化sar圖像變化檢測方法，以在分解前，判斷出極化目標的相干性，降低變化檢測的漏檢率，提高檢測的準確率。

本發(fā)明的技術方案是：對讀入的多時相極化sar數(shù)據(jù)分別進行freeman分解得到較好的散射功率特征，并用其構造極化圖像序列作為稀疏低秩模型的輸入，利用低秩稀疏模型對輸入進行分解獲得稀疏圖像序列，并對其進行融合，得到可分性較好的差異圖，最后使用模糊c均值聚類算法進行分割，得到變化檢測的結果圖，其實現(xiàn)步驟包括如下：

(1)分別從第一時相極化數(shù)據(jù)中提取第一相干矩陣t1，第二時相極化數(shù)據(jù)中提取第二相干矩陣t2，作為freeman分解的兩個輸入；

(2)根據(jù)第一相干矩陣t1和第二相干矩陣t2，得大小均為c＝m×n的第一時相的輸入圖像i1和第二時相的輸入圖像i2；

(3)利用第一時相圖像i1，第二時相圖像i2構造(k-2)幅變化圖像ii；并用i1、i2和ii組成k幅變化圖像序列i＝[i1,…ii,…ik]，其中，i＝2,3,…,k-1,k為圖像序列中圖像的個數(shù)，k≥30，i1對應第一時相圖像i1經(jīng)變換后的列向量，ii對應構造圖像ii經(jīng)變換后的列向量，ik對應第二時相圖像i2經(jīng)變換后的列向量，且i1，ii，ik∈r^c×1，i∈r^c×k，r^c×k表示大小為c×k的實數(shù)空間；

(4)用低秩稀疏分解法將圖像序列i分解為三個子圖像序列之和：i＝l+s+g，其中，s為稀疏圖像序列，l為低秩圖像序列，g為噪聲圖像序列，l,s,g∈r^c×k，且s＝{s1,…sl,…,sk}，sl是稀疏圖像序列s中第l個列向量，l＝1,2,…,k；

(5)用加權均值融合的方法，對稀疏圖像序列s進行融合，得到差異圖，差異圖的大小為c×1；

(6)用模糊c均值方法，對差異圖進行聚類，得到最后的變化檢測結果圖。

本發(fā)明與現(xiàn)有的技術相比具有以下優(yōu)點：

1.本發(fā)明針對極化sar的散射特性，利用freeman分解將多時相極化sar數(shù)據(jù)分解為三個散射功率特征，并針對不同地物，選擇不同的散射功率作為變化檢測的輸入，提高了變化檢測的準確率；

2.本發(fā)明中使用的低秩稀疏方法，既考慮了多時相圖像內(nèi)的鄰域信息，又考慮了多時相圖像之間的差異信息，降低了變化檢測的漏檢率；

3.本發(fā)明中使用的加權平均的方法，由于稀疏圖像序列中的每一幅圖像僅能檢測中部分的變化信息，因而，提高了差異圖的可分性，進而提高了變化檢測的準確率。

附圖說明

圖1是本發(fā)明的實現(xiàn)流程圖；

圖2是本發(fā)明中第一時相和第二時相的pauli圖；

圖3是本發(fā)明仿真所使用的日本東京機場跑道極化sar圖像數(shù)據(jù)集；

圖4是通過人工標記得到圖2的標準變化檢測結果圖；

圖5是用傳統(tǒng)和本發(fā)明變化檢測方法對圖2的變化檢測結果圖；

圖6是本發(fā)明仿真所使用的日本東京部分植被區(qū)遙感圖像數(shù)據(jù)集；

圖7是通過人工標記得到圖5的標準變化檢測結果圖；

圖8是用傳統(tǒng)和本方法變化檢測方法對圖5的變化檢測結果圖。

具體實施方式

以下結合附圖，對本發(fā)明的實施例和效果做進一步說明：

參照圖1，本發(fā)明基于散射特征與低秩稀疏模型的極化sar變化檢測方法，包括如下步驟：

步驟1:從兩時相的極化sar數(shù)據(jù)中提取相應的兩時相相干矩陣。

本實例使用從空載衛(wèi)星alos獲取的兩時相極化sar數(shù)據(jù)，即第一時相極化sar數(shù)據(jù)和第二時相極化sar數(shù)據(jù)；

1a)從第一時相極化sar數(shù)據(jù)中提取第一相干矩陣t1；

1b)從第二時相極化sar數(shù)據(jù)中提取第二相干矩陣t2。

步驟2:根據(jù)第一相干矩陣t1和第二相干矩陣t2，得大小均為c＝m×n的第一時相的輸入圖像i1和第二時相的輸入圖像i2。

2a)利用如下公式分別對第一相干矩陣t1和第二相干矩陣t2分別進行freeman分解：

pv＝8fv/3，其中fv為體散射分量的分解系數(shù)

pd＝fd(1+α²)，其中fv為二面角散射分量的分解系數(shù)，α為常數(shù)

ps＝fs(1+β²)，其中fs為平面散射分量的分解系數(shù)，β為常數(shù)

得第一時相散射功率圖像pd1、pv1、ps1和第二時相散射功率圖像pd2、pv2、ps2；

2b)對第一時相散射功率圖像和第二時相散射功率圖像分別進行配準，得到配準后第一時相散射功率圖像pd1'、pv1'、ps1'和第二時相散射功率圖像pd2'、pv2'、ps2'；

2c)根據(jù)不同區(qū)域的不同地物對應不同的散射功率圖像特征進行特征圖像選擇：

從第一時相的散射功率圖像pd1'、pv1'、ps1'中選出能反應變化信息的散射功率特征圖像，記為第一時相的輸入圖像i1；

從第二時相的散射功率圖像pd2'、pv2'、ps2'中選出能反應變化信息的散射功率特征圖像，記為第二時相的輸入圖像i2。

步驟3:構造k-2幅變化圖像。

3a)用對數(shù)比值法，從第一時相圖像i1和第二時相圖像i2中提取初始變化區(qū)域i0；

3b)利用第一時相圖像i1，第二時相圖像i2以及初始變化區(qū)域i0，得到k-2幅變化圖像ii：

其中，fi表示變化圖像ii中發(fā)生變化的區(qū)域，ui表示變化圖像ii中沒有發(fā)生變化的區(qū)域,表示將初始變化區(qū)域i0劃分為k等份，并取其中的i份，在式中，表示將第二時相圖像i2中的部分變化區(qū)域的灰度值，遷移到變化圖像ii中的對應位置上；表示第一時相圖像i1中除與部分變化區(qū)域位置相同的區(qū)域之外的所有區(qū)域，在式中，表示將第一時相圖像i1中除與部分變化區(qū)域位置相同的區(qū)域之外的所有區(qū)域的灰度值，遷移到變化圖像ii中的對應位置上。

步驟4:利用第一時相圖像i1，第二時相圖像i2和(k-2)幅變化圖像ii，組成k幅變化圖像序列i＝[i1,…ii,…ik]，其中，i＝2,3,…,k-1,k為圖像序列中圖像的個數(shù)，k≥30，i1對應第一時相圖像i1經(jīng)變換后的列向量，ii對應構造圖像ii經(jīng)變換后的列向量，ik對應第二時相圖像i2經(jīng)變換后的列向量，且i1，ii，ik∈r^c×1，i∈r^c×k，r^c×k表示大小為c×k的實數(shù)空間。

步驟5:用低秩稀疏分解法將圖像序列i分解為三個子圖像序列之和。

5a)用雙邊隨機投影的方法，對圖像序列i進行低秩逼近，得到迭代低秩矩陣lt為：

lt＝q1[r1(a2t^ty1t)^-1r2^t]^1/(2q+1)q2^t，

其中，q為收斂系數(shù)，t為迭代次數(shù)，a2t＝ia1(t-1)，當t＝1時，a1(t-1)是初始化秩為r的高斯矩陣，q1,r1分別為對左投影矩陣y1t進行奇異值分解得到的中間參數(shù)：y1t＝ia1t＝q1r1，當t≥2時，a1t＝i^ta2(t-1)，t為轉置符號，q2,r2分別是對右投影矩陣y2t進行奇異值分解得到的中間參數(shù)：y2t＝i^ta2t＝q2r2；

5b)用圖像序列i和迭代低秩矩陣lt去逼近迭代稀疏矩陣st：

st＝pω(i-lt)，

其中，ω表示矩陣(i-lt)的前k個最大非零值子集空間，pω表示將非零數(shù)據(jù)從ω空間投影到迭代稀疏矩陣st上；

5c)設定閾值ε為3.2e-5，重復步驟5a)到5b)，每重復一次，迭代次數(shù)t加1，直至滿足如下終止條件為止：

||i-lt-st||²/||i||²<ε，其中，||||²表示2范數(shù)，

將滿足終止條件的迭代稀疏矩陣st記為稀疏圖像序列s，將迭代低秩矩陣lt記為低秩圖像序列l(wèi)；

5d)由圖像序列i，稀疏圖像序列s和低秩圖像序列l(wèi)，得到噪聲圖像序列g：

g＝i-s-l,

其中，l,s,g∈r^c×k，s＝{s1,…sl,…,sk}，sl是稀疏圖像序列s中第l個列向量，l＝1,2,…,k；

步驟6:根據(jù)稀疏圖像序列s得到變化檢測結果圖。

6a)用加權均值融合的方法，對稀疏圖像序列s進行融合，得到差異圖，差異圖的大小為c×1；

6b)用模糊c均值方法，對差異圖進行聚類，得到每個像素對應的類別標簽，得到最終的變化檢測結果圖。

本發(fā)明的效果可以通過以下仿真實驗進一步說明：

1.實驗條件

實驗環(huán)境為：window7,cpuinteli3,基本頻率為3.2ghz,軟件平臺為matlabr2012b.

仿真使用的數(shù)據(jù)集為兩幅四視的l波段極化sar數(shù)據(jù)，此數(shù)據(jù)集分別拍攝于2006年7月和2009年4月的日本東京地區(qū)，該數(shù)據(jù)集中包含豐富的地物，如海面、城市和植被等。為了便于觀察，從數(shù)據(jù)集中合成出了兩時相極化sar數(shù)據(jù)的rgb圖像:圖2(a)和圖2(b)，它們的尺寸為2300×1048，其中r代表|hh-vv|，g代表|hv|+|vh|，b代表|hh+vv|,h表示衛(wèi)星水平發(fā)送接收方式，v表示垂直發(fā)送接收方式，hv表示水平發(fā)送，垂直接收。

選擇圖2的區(qū)域1和區(qū)域2進行分析。區(qū)域1主要是修建飛機場跑道引起的變化，所選變化區(qū)域的大小為600×500；區(qū)域2主要是植被隨季節(jié)變化引起的變化，所選變化區(qū)域大小為318×230。

將區(qū)域1對應的極化sar數(shù)據(jù)作為第一個數(shù)據(jù)集，如圖3所示，其中圖3(a)為新機場跑道的第一時相的體散射功率圖像，圖3(b)為新機場跑道第二時相的體散射功率圖像。圖像大小均為600×500像素，灰度級為256，其包括290420個非變化像素和9580個變化像素。

對上述東京新機場跑道進行人工標記，得到此跑道標準圖如圖4所示。

將區(qū)域2對應的極化sar數(shù)據(jù)作為第二個數(shù)據(jù)集，如圖6所示，其中圖6(a)為東京部分植被區(qū)域的第一時相二面角散射功率圖像，圖6(b)為東京部分植被地區(qū)的第二時相二面角散射功率圖像，這些圖像大小均為318×230，灰度級為256，其包括69627個非變化像素和3513個變化像素。

對上述東京部分植被區(qū)域進行人工標記，得到此區(qū)域的標準圖如圖7所示。

2.實驗內(nèi)容和實驗結果

實驗1：用本發(fā)明方法和現(xiàn)有差值法sm及現(xiàn)有對數(shù)比值法lr，對圖3進行變化檢測。實驗結果如圖5所示，其中5(a)是sm差值法對圖3進行變化檢測的結果圖，5(b)是lr對數(shù)比值法對圖3進行變化檢測的結果圖，5(c)是本發(fā)明方法對圖3進行變化檢測的結果圖。

實驗2：用本發(fā)明方法與兩種傳統(tǒng)的變化檢測方法：差值法sm和對數(shù)比值法lr，對圖6進行變化檢測，實驗結果如圖8所示。其中8(a)是差值法sm對圖6進行變化檢測的結果圖，8(b)是對數(shù)比值法lr對圖6進行變化檢測的結果圖，8(c)是本發(fā)明方法對圖6進行變化檢測的結果圖。

從圖5(a)和圖8(a)中可以看出，sd差值法對圖像進行變化檢測的實驗結果中有很多的雜點；

從圖5(b)和圖8(b)中可以看出，lr對數(shù)比值法的實驗結果產(chǎn)生了很多的細節(jié)丟失，邊緣很模糊，且圖5(b)丟失了大部分的變化區(qū)域；

從圖5(c)和圖8(c)中可以看出，在較好保持邊緣細節(jié)的情況下，本發(fā)明能減少雜點的個數(shù)，進而降低漏檢率。

對本發(fā)明方法和上述兩種傳統(tǒng)方法在兩個變化區(qū)域上的變化檢測結果進行統(tǒng)計，如表1所示。表1中有四種評價指標：分別為虛警數(shù)fa，漏檢數(shù)ma，錯誤數(shù)oe和正確率pcc，其中，虛警數(shù)fa為實際沒有發(fā)生變化但被當作變化檢測出來的像素，漏檢數(shù)ma為沒有檢測出來的實際發(fā)生了變化的像素，錯誤數(shù)oe＝fa+ma，正確率pcc＝錯誤數(shù)/圖像像素個數(shù)，其中錯誤數(shù)為所用方法的結果圖與標準圖的差值。

表1實驗結果數(shù)據(jù)

從表1可以看出，本發(fā)明與兩種傳統(tǒng)變化檢測方法相比：

首先，從表1的虛檢數(shù)指標中可以看出，本發(fā)明對區(qū)域2的檢測有相對較高的虛檢數(shù)，這與本發(fā)明中的散射功率特征的選取和低秩稀疏分解好壞有一定的關系；

其次，從表1的漏檢指標中可以看出，本發(fā)明方法在區(qū)域1數(shù)據(jù)集以及區(qū)域2數(shù)據(jù)集上與兩種傳統(tǒng)變化檢測方法相比，均獲得較高的準確率，并且均獲得較少的漏檢數(shù)。而且本發(fā)明方法在區(qū)域1數(shù)據(jù)集上與差值法相比，漏檢數(shù)降低128個像素點，在區(qū)域2數(shù)據(jù)集上的；漏檢數(shù)降低457個像素點。

綜上，本發(fā)明方法與兩種傳統(tǒng)變化檢測方法相比，均有較少的漏檢數(shù)，較少的錯誤數(shù)，且無論是在區(qū)域1數(shù)據(jù)集以及區(qū)域2數(shù)據(jù)集中，均獲得較高的準確率，提高了變化檢測精度。