本發(fā)明屬于遙感影像處理技術(shù)領(lǐng)域，特別涉及了基于空間約束的卡方變換的多時相遙感影像變化檢測方法。

背景技術(shù)：

隨著多時相遙感數(shù)據(jù)的不斷積累以及空間數(shù)據(jù)庫的相繼建立，如何從這些遙感數(shù)據(jù)中提取和檢測變化信息已成為遙感科學(xué)和地理信息科學(xué)的重要研究課題。根據(jù)同一區(qū)域不同時相的遙感影像，可以提取城市、環(huán)境等動態(tài)變化的信息，為資源管理與規(guī)劃、環(huán)境保護等部門提供科學(xué)決策的依據(jù)。

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

針對上述問題，有必要在檢測過程中加入空間約束來提高檢測精度，另一方面，對關(guān)鍵參數(shù)的自適應(yīng)選擇也是檢測技術(shù)中的重點。

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

為了解決上述背景技術(shù)提出的技術(shù)問題，本發(fā)明旨在提供基于空間約束的卡方變換的多時相遙感影像變化檢測方法，采用空間約束的卡方變換結(jié)合自適應(yīng)置信水平的變化檢測方法，克服了多時相多光譜遙感影像背景信息復(fù)雜、噪聲干擾嚴(yán)重的問題。

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

基于空間約束的卡方變換的多時相遙感影像變化檢測方法，包括以下步驟：

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

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

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

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

(5)計算D_X的模值X_M，利用貝葉斯原理，并基于最大期望算法獲取最優(yōu)分割閾值T，將|X_M-T|≤δ的區(qū)域作為偽訓(xùn)練樣本集，其中δ取X_M動態(tài)范圍的百分比；

(6)設(shè)定置信水平的搜索范圍和搜索步長；

(7)確定D_X中的非變化區(qū)域，計算非變化區(qū)域的均值矢量和方差矩陣，并計算D_X上每個點的卡方值；初始化時，將整個D_X作為非變化區(qū)域；

(8)在給定的置信水平A的基礎(chǔ)上，計算檢測閾值，并根據(jù)該閾值進行卡方檢測，獲得初步的檢測結(jié)果M₀；

(9)給定濾波窗口大小，對M₀進行眾數(shù)濾波，即濾波窗口內(nèi)，如果變化的像素數(shù)目大于非變化的，則窗口中心點像素為變化的，反之，窗口中心點像素為非變化的，記眾數(shù)濾波的結(jié)果為M₁；

(10)判斷M₁中的非變化區(qū)域相較于M₀中的非變化區(qū)域是否有改變，如果兩者沒有改變，將給定的置信水平A加上步驟(6)設(shè)定的搜索步長作為新的置信水平，并返回步驟(7)；如果兩者改變了，則將M₁中的非變化區(qū)域作為新的非變化區(qū)域，返回步驟(7)；

(11)當(dāng)置信水平達到步驟(6)設(shè)定的搜索范圍的上界，終止循環(huán)迭代；

(12)針對每一個置信水平，計算偽訓(xùn)練樣本集的精度，選擇精度最高的偽訓(xùn)練樣本集對應(yīng)的置信水平，并在此基礎(chǔ)上，輸出最終的變化檢測結(jié)果。

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

(a)假設(shè)X_M影像上的非變化類ω_n和變化類ω_c服從如下的高斯分布：

ω_i∈{ω_n,ω_c}，m_i∈{m_n,m_c}，σ_i∈{σ_n,σ_c}上式中，m_n和σ_n分別為非變化類的均值和方差，m_c和σ_c分別為變化類的均值和方差；

(b)采用最大期望算法算法估計m_n、σ_n、m_c和σ_c這四個參數(shù)；

(c)依據(jù)貝葉斯最小誤差準(zhǔn)則，求解X_M的分割閾值T，在高斯分布的情況下，等同于求解下式：

上式中，為分割閾值T的估計值，p(ω_n)、p(ω_c)分別為非變化類先驗概率和變化類先驗概率；

(d)根據(jù)估計的分割閾值，構(gòu)建偽訓(xùn)練樣本集：

非變化區(qū)域偽訓(xùn)練樣本集：

變化區(qū)域偽訓(xùn)練樣本集：

其中，δ取X_M動態(tài)范圍的15％。

進一步地，在步驟(6)中，設(shè)定置信水平的搜索范圍為0.95-0.999，搜索步長為0.001。

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

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

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

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

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

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

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

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

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

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

(1)本發(fā)明將眾數(shù)濾波引入到CST檢測結(jié)果之后，使得使用CST變換具有空間約束能力；

(2)本發(fā)明采用迭代的方法，估計非變化區(qū)域的均值和標(biāo)準(zhǔn)差，克服了僅僅利用整個多波段影像Dx計算均值和方差矩陣的不足，可以使變化檢測的結(jié)果更加可靠，也更加具有穩(wěn)健性；

(3)在本發(fā)明中，最優(yōu)置信水平是基于偽訓(xùn)練樣本集獲取，使得該技術(shù)通用性較好。

附圖說明

圖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所示，基于空間約束的卡方變換的多時相遙感影像變化檢測方法，具體步驟如下：

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

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

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

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

(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₂)四個點的值，如果選擇一個坐標(biāo)系統(tǒng)使得這四個點的坐標(biāo)分別為(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)的基準(zhǔn)影像同名點的地理坐標(biāo)來確定該三角形內(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的模值X_M，利用貝葉斯原理，并基于最大期望算法獲取最優(yōu)分割閾值T，將|X_M-T|≤δ的區(qū)域作為偽訓(xùn)練樣本集，其中δ取X_M動態(tài)范圍的百分比。具體過程如下：

(a)假設(shè)X_M影像上的非變化類ω_n和變化類ω_c服從如下的高斯分布：

ω_i∈{ω_n,ω_c}，m_i∈{m_n,m_c}，σ_i∈{σ_n,σ_c}上式中，m_n和σ_n分別為非變化類的均值和方差，m_c和σ_c分別為變化類的均值和方差；

(b)采用最大期望算法算法估計m_n、σ_n、m_c和σ_c這四個參數(shù)，下面以非變化類的參數(shù)估計為例進行說明，變化類參數(shù)估計類似：

上式中，I和J分別表示X_M的行數(shù)和列數(shù)，上標(biāo)t表示迭代次數(shù)；

(c)依據(jù)貝葉斯最小誤差準(zhǔn)則，求解X_M的分割閾值T，在高斯分布的情況下，等同于求解下式：

上式中，為分割閾值T的估計值，p(ω_n)、p(ω_c)分別為非變化類先驗概率和變化類先驗概率；

(d)根據(jù)估計的分割閾值，構(gòu)建偽訓(xùn)練樣本集：

非變化區(qū)域偽訓(xùn)練樣本集：

變化區(qū)域偽訓(xùn)練樣本集：

其中，δ取X_M動態(tài)范圍的15％。

步驟6：設(shè)定置信水平的搜索范圍為0.95-0.999，以及置信水平的搜索步長0.001。

步驟7：確定D_X中的非變化區(qū)域(初始化時，將整個D_X作為非變化區(qū)域)，計算非變化區(qū)域的均值矢量m和方差矩陣Σ，并計算D_X上每個點的卡方值：

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

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

步驟8：在給定的置信水平A的基礎(chǔ)上，計算檢測閾值，并根據(jù)該閾值進行卡方檢測，獲得初步的檢測結(jié)果M₀。檢測閾值的計算方法如下：

當(dāng)置信水平為1-α?xí)r，C_ij的值大于的概率為α。如果α取值較小，則大于的C_ij可以視為出界點(outlier)或者變化點，由此確定檢測閾值為

步驟9：給定濾波窗口大小(一般設(shè)置為3×3或者5×5的窗口)，對M₀進行眾數(shù)濾波，即濾波窗口內(nèi)，如果變化的像素數(shù)目大于非變化的，則窗口中心點像素為變化的，反之，窗口中心點像素為非變化的，記眾數(shù)濾波的結(jié)果為M₁。

步驟10：判斷M₁中的非變化區(qū)域相較于M₀中的非變化區(qū)域是否有改變，如果兩者沒有改變，將給定的置信水平A加上步驟6設(shè)定的搜索步長作為新的置信水平，并返回步驟7；如果兩者改變了，則將M₁中的非變化區(qū)域作為新的非變化區(qū)域，返回步驟7。

步驟11：當(dāng)置信水平達到0.999，終止循環(huán)迭代。

步驟12：針對每一個置信水平，計算偽訓(xùn)練樣本集的精度，選擇精度最高的偽訓(xùn)練樣本集對應(yīng)的置信水平，并在此基礎(chǔ)上，輸出最終的變化檢測結(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為變化檢測的參考圖像。上述四種方法的檢測結(jié)果如圖4(a)、4(b)、4(c)、4(d)所示，檢測性能用錯檢數(shù)FP、漏檢數(shù)FN、總錯誤數(shù)OE和Kappa系數(shù)四個指標(biāo)來衡量。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)。