本發(fā)明屬于遙感影像處理技術(shù)領(lǐng)域，特別涉及了卡方變換結(jié)合MRF模型的多時相遙感影像變化檢測方法。

背景技術(shù)：

隨著多時相高分辨率遙感數(shù)據(jù)的不斷積累以及空間數(shù)據(jù)庫的相繼建立，如何從這些遙感數(shù)據(jù)中提取和檢測變化信息已成為遙感科學和地理信息科學的重要研究課題。根據(jù)同一區(qū)域不同時相的遙感影像，可以提取城市、環(huán)境等動態(tài)變化的信息，為資源管理與規(guī)劃、環(huán)境保護等部門提供科學決策的依據(jù)。我國“十二五”將加大拓展實施“十一五”已啟動實施的高分辨率對地觀測工程，關(guān)注包括高分辨率遙感目標與空間環(huán)境特征分析及高可靠性自動解譯等基礎(chǔ)理論與關(guān)鍵技術(shù)研究，正在成為解決國家安全和社會經(jīng)濟發(fā)展重大需求的研究焦點。

遙感影像的變化檢測就是從不同時期的遙感數(shù)據(jù)中，定量地分析和確定地表變化的特征與過程。各國學者從不同的角度和應(yīng)用研究提出了許多有效的檢測算法，如變化矢量分析法(Change Vector Analysis，CVA)、基于Fuzzy C-means(FCM)的聚類方法等。其中，傳統(tǒng)的基于卡方變換(Chi-Squared Transform,CST)的多時相光學遙感變化檢測，先計算差異影像的均值和方差矩陣，然后再基于置信水平，確定變化檢測的閾值，進而得到變化檢測結(jié)果。該類技術(shù)中，使用CST的不足是僅使用多時相高分辨率差異影像的光譜信息，沒有利用空間信息。另外，在計算差異影像的均值和方差矩陣時，估計的精度不高。

針對上述問題，有必要研究新的高分辨率可見光遙感圖像變化檢測技術(shù)來有效克服上述難點。

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

為了解決上述背景技術(shù)提出的技術(shù)問題，本發(fā)明旨在提供卡方變換結(jié)合MRF模型的多時相遙感影像變化檢測方法，克服了現(xiàn)有技術(shù)難以解決高空間分辨率遙感影像背景信息復雜、噪聲干擾嚴重的問題。

為了實現(xiàn)上述技術(shù)目的，本發(fā)明的技術(shù)方案為：

卡方變換結(jié)合MRF模型的多時相遙感影像變化檢測方法，包括以下步驟：

(1)輸入兩時相的高分辨率光學遙感影像，分別記為X₁和X₂；

(2)對X₁和X₂進行影像配準；

(3)利用多元變化檢測方法對X₁和X₂進行輻射歸一化校正；

(4)計算多時相差異影像D_X＝X₁-X₂；

(5)初始化多時相差異影像D_X中的非變化區(qū)域，計算非變化區(qū)域的均值矢量和方差矩陣，并計算多時相差異影像上每個點的卡方值；

(6)在給定的置信水平的基礎(chǔ)上，計算檢測閾值，并根據(jù)檢測閾值進行檢測，確定多時相差異影像中的變化區(qū)域和非變化區(qū)域；

(7)將步驟(6)確定的非變化區(qū)域與步驟(5)確定的非變化區(qū)域進行比較，如果兩者一致，則將步驟(6)得到的檢測結(jié)果作為變化檢測結(jié)果，如果兩者不同，則將步驟(6)確定的非變化區(qū)域作為新的非變化區(qū)域，返回步驟(5)，循環(huán)迭代；

(8)將步驟(7)得到的變化檢測結(jié)果作為MRF模型的輸入，并基于多時相差異影像D_X的模值得到最終的變化檢測結(jié)果。

進一步地，步驟(8)的具體過程如下：

(a)計算多時相差異影像D_X的模值：

上式中，X_1b和X_2b分別表示第X₁和X₂第b波段的影像，B表示每一個時相遙感影像的波段數(shù)目，(i,j)是影像的坐標；

(b)構(gòu)建MRF模型的能量函數(shù)：

上式中，U_data表示數(shù)據(jù)項，U_context表示空間局部能量項，M₁和M₂分別表示影像的高和寬，Y(i,j)表示變化檢測結(jié)果坐標(i,j)的值，Y^S(i,j)是坐標(i,j)的鄰域系統(tǒng)；

(c)采用ICM優(yōu)化算法求解U最小化，得到最終的變化檢測結(jié)果。

進一步地，U_data進一步表示為，

上式中，μ_Y(i,j)∈{μ_n,μ_c}，和μ_n分別表示非變化區(qū)域的方差和均值，和μ_c分別表示變化區(qū)域的方差和均值。

進一步地，其特征在于，U_context進一步表示為，

上式中，(p,q)是(i,j)的鄰域坐標，β是控制空間局部能量項的參數(shù)，I是指示函數(shù)，定義如下：

進一步地，在步驟(2)中，對X₁和X₂進行影像配準包括幾何粗校正和幾何精校正，所述幾何粗校正的過程：

(A)選擇X₁和X₂分別作為基準影像和待校正影像；

(B)在基準影像和待校正影像上分別采集地面控制點，地面控制點的數(shù)量大于等于9，且地面控制點均勻分布在影像上；

(C)計算基準影像和待校正影像各地面控制點處的均方誤差；

(D)采用多項式糾正法對待校正影像進行糾正；

(E)采用雙線性插值法對待校正影像進行重采樣；

所述幾何精校正是將經(jīng)過幾何粗校正的遙感影像，利用自動匹配與三角剖分法進行校正。

進一步地，在步驟(5)中，采用下式計算檢測閾值：

上式中，1-α是置信水平，C_ij表示D_X在坐標(i，j)處的卡方值，為檢測閾值。

采用上述技術(shù)方案帶來的有益效果：

本發(fā)明將MRF模型引入到CST檢測結(jié)果之后，使得基于CST變換的變化檢測具有空間約束能力。在變化檢測中，采用迭代的方法，估計非變化區(qū)域的均值和標準差，克服了僅僅利用整個多波段差異影像計算均值和方差矩陣的不足，使得變化檢測的結(jié)果更加可靠，也更加具有穩(wěn)健性。

附圖說明

圖1是本發(fā)明的方法流程圖；

圖2(a)、2(b)分別是2007年1月的沙特阿拉伯Mina地區(qū)高分辨率IKONOS圖像第3波段示意圖、2007年12月的沙特阿拉伯的Mina地區(qū)高分辨率IKONOS圖像第3波段示意圖；

圖3是變化檢測的參考圖像；

圖4(a)、4(b)、4(c)、4(d)分別是EM-CVA算法、ICST算法、RCST算法、本發(fā)明算法的檢測結(jié)果示意圖。

具體實施方式

以下將結(jié)合附圖，對本發(fā)明的技術(shù)方案進行詳細說明。

如圖1所示，卡方變換結(jié)合MRF模型的多時相遙感影像變化檢測方法，步驟如下：

步驟1：輸入同一區(qū)域、不同時相的兩幅高分辨率光學遙感影像，分別記為：X₁和X₂。

步驟2：對X₁和X₂進行影像配準，分為粗校正和精校正兩個步驟：

對于幾何粗校正，利用ENVI4.8軟件中的相關(guān)功能實現(xiàn)，具體操作步驟為：

(1)顯示基準影像和待校正影像；

(2)采集地面控制點GCPs，GCPs應(yīng)均勻分布在整幅圖像內(nèi)，GCPs的數(shù)目至少大于等于9；

(3)計算均方誤差；

(4)采用多項式糾正法對待校正影像進行糾正；

(5)采用雙線性插值進行重采樣輸出，若求未知函數(shù)f在點P＝(x,y)的值，假設(shè)已知函數(shù)f在Q₁₁＝(x₁,y₁),Q₁₂＝(x₁,y₂),Q₂₁＝(x₂,y₁),及Q₂₂＝(x₂,y₂)四個點的值，如果選擇一個坐標系統(tǒng)使得這四個點的坐標分別為(0,0)、(0,1)、(1,0)和(1,1)，那么雙線性插值公式就可以表示為：

f(x,y)≈f(0,0)(1-x)(1-y)+f(1,0)x(1-y)+f(0,1)(1-x)y+f(1,1)xy對于幾何精校正，將經(jīng)過幾何粗校正的多光譜遙感影像數(shù)據(jù)，利用自動匹配與三角剖分法進行幾何精校正。采用逐點插入法構(gòu)建Delaunay三角網(wǎng)，對每一個三角形，利用其三個頂點的行列號與其對應(yīng)的基準影像同名點的地理坐標來確定該三角形內(nèi)部的仿射變換模型參數(shù)，對待校正影像進行糾正，得到校正后的遙感影。

步驟3：利用多元變化檢測(Multivariate Alteration Detection，MAD)方法對X₁和X₂進行輻射歸一化校正。首先找到X₁和X₂各波段亮度值的一個線性組合，得到變化信息增強的差異影像，通過閾值確定變化和未變化區(qū)域，然后通過未變化區(qū)域?qū)?yīng)的兩時相像元對的映射方程，完成相對輻射校正。

步驟4：對輸入的多時相高分辨率影像X₁和X₂，計算多時相差異影像D_X：

D_X＝X₁-X₂

步驟5：將整個D_X視為非變化區(qū)域，并計算其均值矢量m和方差矩陣Σ，并計算差異影像每一個點的卡方值：

C_ij＝(x_ij-m)^TΣ^-1(x_ij-m)～χ²(b)

上式中，C_ij表示差異影像(i，j)坐標點的卡方值，其服從自由度為b的卡方分布；x_ij表示差異影像在(i，j)坐標點的矢量值；Σ^-1表示方差矩陣的逆矩陣；b表示差異影像的波段數(shù)目。

步驟6：給定置信水平1-α，利用下式計算檢測閾值：

當置信水平為1-α時，C_ij的值大于的概率為α。如果α取值較小，則大于的C_ij可以視為出界點(outlier)或者變化點，由此確定閾值為并根據(jù)該閾值獲得初步的檢測結(jié)果。

步驟7：將步驟6檢測結(jié)果中確定的非變化區(qū)域與步驟5確定的非變化區(qū)域(整個D_X)進行比較，如果兩者一致，則將步驟6得到的檢測結(jié)果作為變化檢測結(jié)果，如果兩者不同，則將步驟6確定的非變化區(qū)域作為新的非變化區(qū)域，返回步驟5，循環(huán)迭代.

步驟8：將步驟7得到的變化檢測結(jié)果作為MRF模型的輸入，并基于差異影像的模值獲取最終的變化檢測結(jié)果。具體實現(xiàn)步驟為：

(1)首先計算D_X的模值X_M：

上式中，X_1b和X_2b分別表示第一時相和第二時相多光譜影像X₁和X₂第b波段影像，B表示每一個時相遙感影像的波段數(shù)目，(i,j)是影像的坐標；

(2)構(gòu)建MRF模型的能量函數(shù)如下：

上式中，U_data表示數(shù)據(jù)項，U_context表示空間局部能量項，M₁和M₂分別表示影像的高和寬，Y(i,j)表示變化檢測結(jié)果坐標(i,j)的值，Y^S(i,j)是坐標(i,j)的鄰域系統(tǒng)。

其中，數(shù)據(jù)項U_data可以進一步表示為：

上式中，μ_Y(i,j)∈{μ_n,μ_c}，和μ_n分別表示非變化區(qū)域的方差和均值，和μ_c分別表示變化區(qū)域的方差和均值。

其中，局部局部空間能量項U_context可以進一步表示為：

上式中，(p,q)是鄰域坐標。本發(fā)明中采用8鄰域系統(tǒng)，β是控制空間局部能量項的參數(shù)，I為指示函數(shù)，定義如下：

(3)采用ICM優(yōu)化算法求解U最小化，得到最終的變化檢測結(jié)果。

本發(fā)明的效果可通過以下實驗結(jié)果與分析進一步說明：

本發(fā)明的實驗數(shù)據(jù)為沙特阿拉伯的Mina地區(qū)的多時相IKNOS高分辨影像數(shù)據(jù)，圖像大小為700×950，使用B1、B2和B3三個波段，圖2(a)和圖2(b)為兩時像B3波段的遙感影像。

為了驗證本發(fā)明的有效性，將本發(fā)明變化檢測方法與下述變化檢測方法進行比對：

(1)基于CVA的EM方法(CVA-EM)[意大利的Bruzzone L.等在文章“Automatic analysis of difference image for unsupervised change detection”(IEEE Transactions on Geoscience and Remote Sensing,2000,38(3):1171-1182.)中所提的檢測方法]。

(2)基于迭代的CST檢測(ICST)方法[B.Desclée,P.Bogaert,and P.Defourny在文章“Forest change detection by statistical object-based method”(Remote Sensing of Environment,2006,102(1-2):1-12.)中所提的方法]。

(3)基于魯棒估計的CST檢測(RCST)方法[Aiye Shi等在文章“Unsupervised change detection based on robust chi-squared transform for bitemporal remotely sensed images”(International of Remote Sensing,2006,102(1-2):1-12.)中所提的方法]。

(4)本發(fā)明方法。

圖3為變化檢測的參考圖像。在置信水平0.99的條件下，上述四種方法的檢測結(jié)果如圖4(a)、4(b)、4(c)、4(d)所示，檢測性能用錯檢數(shù)FP、漏檢數(shù)FN、總錯誤數(shù)OE和Kappa系數(shù)四個指標來衡量。FP、FN和OE越接近于0、Kappa系數(shù)越接近于1，表明變化檢測方法的性能越好。檢測結(jié)果如表1所示。由表1可見，本發(fā)明所提的檢測方法性能優(yōu)于其他三種檢測方法，這表明本發(fā)明所提的變化檢測方法是有效的。

表1

以上實施例僅為說明本發(fā)明的技術(shù)思想，不能以此限定本發(fā)明的保護范圍，凡是按照本發(fā)明提出的技術(shù)思想，在技術(shù)方案基礎(chǔ)上所做的任何改動，均落入本發(fā)明保護范圍之內(nèi)。