本發(fā)明涉及圖像處理技術(shù)，特別涉及一種結(jié)合dt-cwt和mrf的遙感圖像變化檢測方法與裝置。

背景技術(shù)：

隨著遙感數(shù)據(jù)獲取手段不斷提高和更新周期不斷縮短，變化檢測技術(shù)作為遙感圖像處理與分析中的一項重要應(yīng)用，已廣泛應(yīng)用于土地利用和土地覆蓋變化監(jiān)測，森林或植被變化檢測，災(zāi)害救援與治理、城鎮(zhèn)和道路變化、地理數(shù)據(jù)庫更新等眾多領(lǐng)域。由于無監(jiān)督變化檢測方法除了利用兩時相的原始圖像外，不需要任何額外的先驗知識就能得到最終變化檢測結(jié)果，方法較為高效，故成為研究的熱點(diǎn)。通常無監(jiān)督變化檢測是基于像素差異的變化檢測，可分為閾值法、變換法、圖像建模法等。閾值法是通過選擇合適的閾值來區(qū)分差異圖像中的變化類與非變化類像素，從而得到最終的變化區(qū)域，具有簡單直觀的優(yōu)點(diǎn)，但是往往只進(jìn)行單一尺度的分割，難免會引起過度分割或不完全分割，同時檢測結(jié)果嚴(yán)重依賴于選取的閾值。變換法主要是多尺度變換法方法，如小波變換，這類方法可以進(jìn)行多尺度分析，能較好地克服了傳感器噪聲、配準(zhǔn)誤差等因素的影響，但沒有考慮像素間的相關(guān)性，檢測邊緣粗糙，誤檢像素較多。圖像建模法是通過建立模型來模擬空間上下文相關(guān)的分布信息，如馬爾科夫隨機(jī)場(markovrandomfield，mrf)充分考慮了鄰域內(nèi)像素之間的相互影響，較好地克服了孤立像素問題，但突變和細(xì)節(jié)信息表達(dá)較弱，偽變化檢測較多。為了解決上述方法存在的問題，已經(jīng)提出了將小波變換與mrf相結(jié)合的遙感圖像變化檢測方法，這些方法利用了小波變換的時頻分析能力并且考慮了像素的空間相關(guān)性，較好地克服了單一像素孤立性、噪聲和配準(zhǔn)誤差等因素引起的偽變化，但是由于在融合過程中直接舍棄了高頻信息，只取低頻信息進(jìn)行試驗，難免造成信息的損失，加之小波變換存在著方向選擇性和平移敏感性，使得變化檢測細(xì)節(jié)仍表達(dá)不足。

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

本發(fā)明旨在提供一種結(jié)合雙樹復(fù)小波變換(dual-treecomplexwavelettransform，dt-cwt)和馬爾科夫隨機(jī)場(mrf)的遙感圖像變化檢測方法與裝置，其不僅能有效平衡高頻信息的保留和噪聲的去除，而且邊緣檢測更加平滑，變化檢測結(jié)果具有較好的區(qū)域一致性，提高遙感圖像變化檢測的精度，可充分地消除由于現(xiàn)有技術(shù)的限制和缺陷導(dǎo)致的一個或多個問題。

本發(fā)明另外的優(yōu)點(diǎn)、目的和特性，一部分將在下面的說明書中得到闡明，而另一部分對于本領(lǐng)域的普通技術(shù)人員通過對下面的說明的考察將是明顯的或從本發(fā)明的實施中學(xué)到。通過在文字的說明書和權(quán)利要求書及附圖中特別地指出的結(jié)構(gòu)可實現(xiàn)和獲得本發(fā)明目的和優(yōu)點(diǎn)。

本發(fā)明提供了一種結(jié)合dt-cwt和mrf的遙感圖像變化檢測方法，其特征在于，所述方法具體包括以下步驟：

步驟1，輸入大小相同的經(jīng)過精確配準(zhǔn)的兩期同區(qū)域的遙感影像對應(yīng)的兩幅灰度圖像x1、x2；

步驟2，求所述兩幅灰度圖像x1、x2的差值圖像；

步驟3，對所述差值圖像進(jìn)行3層dt-cwt分解；

步驟4，利用基于fcm的mrf分割算法，對經(jīng)dt-cwt分解后得到的高頻子帶進(jìn)行變化特征提取；

步驟5，對通過步驟3的dt-cwt分解得到的低頻分量和通過步驟4的變化特征提取得到的高頻分量進(jìn)行重構(gòu)；

步驟6，利用所述基于fcm的mrf分割算法，對通過步驟5重構(gòu)后的結(jié)果進(jìn)行變化檢測；

步驟7，基于相似度融合規(guī)則，對通過步驟6得到的檢測結(jié)果進(jìn)行融合；

步驟8，輸出融合后的變化檢測結(jié)果。

優(yōu)選的，步驟2具體包括：

根據(jù)公式d(i,j)＝|x2(i,j)-x1(i,j)|求所述差值圖像d，其中，d(i,j)為差值圖像d在位置(i，j)處的像素值，x1(i,j)和x2(i,j)分別是灰度圖像x1和x2的在位置(i，j)處的像素值，i，j分別為1≤i≤p，1≤j≤q的整數(shù)，p為所述輸入圖像的行數(shù)，q為所述輸入圖像的列數(shù)。

優(yōu)選的，步驟3具體包括：采用兩組并行的基于正交變換的q-shift(quartersampleshift)實數(shù)濾波器組來對步驟2中得到的差值圖像d進(jìn)行3層dt-cwt分解。第一層分解采用相差一個單位延遲的奇數(shù)長度的濾波器組，分解結(jié)果的一樹作為第一層分解復(fù)小波系數(shù)的實部，另一樹作為第一層分解復(fù)小波系數(shù)的虛部；第二、三層分解均采用遲近似為1/4采樣間隔的偶數(shù)長度的濾波器組，分解結(jié)果一樹作為第二、三層分解復(fù)小波系數(shù)的實部，另一樹作為第二、三層分解復(fù)小波系數(shù)的虛部。

優(yōu)選的，步驟4具體包括以下子步驟：

步驟4.1，對所述經(jīng)過3層dt-cwt分解得到的高頻子帶hk,l，(k＝1,2,3；l＝±15°,±45°,±75°)中的第一層高頻分量h1,l(l＝±15°,±45°,±75°)置零，即得到第一層高頻變化特征提取結(jié)果，記為

步驟4.2，利用基于bayes閾值的fcm算法對所述經(jīng)過3層dt-cwt分解得到的高頻子帶hk,l(k＝1,2,3；l＝±15°,±45°,±75°)中的第二、三層高頻分量hk,l(k＝2,3；l＝±15°,±45°,±75°)進(jìn)行初始聚類分割，得到初始分割結(jié)果

步驟4.3，基于mrf與grf的等價性，利用條件迭代模型icm對通過步驟4.2獲得的初始分割結(jié)果進(jìn)行迭代更新，得到最終的最優(yōu)分割，記為

步驟4.4，對通過步驟4.3得到的最終分割結(jié)果進(jìn)行二進(jìn)制掩膜，得到掩膜圖bk,l(k＝2,3；l＝±15°,±45°,±75°)；

步驟4.5，將所述掩膜圖bk,l(k＝2,3；l＝±15°,±45°,±75°)矩陣像素值為1的記為變化類c，值為0的記為未變化類u；

步驟4.6，將c類像素對應(yīng)位置的二、三層高頻子帶hk,l(k＝2,3；l＝±15°,±45°,±75°)的像素值保持不變，u類像素對應(yīng)位置的二、三層高頻子帶的像素值置零，得到提取變化特征后的二、三層高頻子帶

優(yōu)選的，步驟4.2具體包括以下子步驟：

步驟4.2.1，根據(jù)最小錯誤概率準(zhǔn)則，求自適應(yīng)bayes閾值level；

步驟4.2.2，根據(jù)v1＝(min(hk,l(i,j))+level)/2、v2＝(level+max(hk,l(i,j))/2，分別對所述經(jīng)過3層dt-cwt分解得到的第二、三層高頻分量hk,l，(k＝2,3；l＝±15°,±45°,±75°)求初始聚類中心v1、v2，記vm＝[v1,v1],其中m為迭代次數(shù)，level為自適應(yīng)bayes閾值，max(hk,l(i,j))和min(hk,l(i,j))分別表示每個子帶hk,l，(k＝2,3；l＝±15°,±45°,±75°)中的像素最大和最小值，hk,l(i,j)表示位置(i,j)處的像素值，1≤i≤p，1≤j≤q的整數(shù)，p為hk,l(i,j)對應(yīng)的行數(shù)，q為hk,l(i,j)對應(yīng)的列數(shù)；

步驟4.2.3，根據(jù)式

對每個子帶hk,l(k＝2,3；l＝±15°,±45°,±75°)計算隸屬度矩陣u，其中，c為聚類中心數(shù)，n為權(quán)值，||hj-cj||為樣本hj到第i類聚類中心的距離度量，||hj-ck||為樣本hj到第k類聚類中心的距離度量；

步驟4.2.4，根據(jù)步驟4.2.3中得到的隸屬度矩陣u，求初始標(biāo)記場，并更新步驟4.2.2的聚類中心vm，m為迭代次數(shù)；

步驟4.2.5，重復(fù)步驟4.2.3、4.2.4不斷對所述標(biāo)記場進(jìn)行迭代更新，直到達(dá)到迭代終止條件(優(yōu)選的，迭代終止條件為15次)；

步驟4.2.6，輸出迭代最終的標(biāo)記場，也即是初始分割結(jié)果

優(yōu)選的，步驟4.3具體包括以下子步驟：

步驟4.3.1，根據(jù)步驟4.2的fcm初始分割結(jié)果計算特征場參數(shù)：均值μm和方差

步驟4.3.2，利用步驟4.3.1計算的特征場參數(shù)μm和方差來計算當(dāng)前分割的每個像素的特征場能量，并計算當(dāng)前分割的每個像素的標(biāo)記場能量；

步驟4.3.3，根據(jù)特征場和標(biāo)記場總能量之和最小值所在變化或未變化的類別來更新當(dāng)前分割結(jié)果，并記為m為迭代次數(shù)；

步驟4.3.4，重復(fù)步驟4.3.1-4.3.3，不斷對進(jìn)行迭代更新，直到達(dá)到循環(huán)終止條件(優(yōu)選的，循環(huán)終止條件設(shè)置為30次)，獲得最終的mrf-map最優(yōu)分割結(jié)果m為迭代次數(shù)。

步驟4.3.5，輸出所述最終的mrf-map最優(yōu)分割結(jié)果

優(yōu)選的，在步驟5中，對步驟3中dt-cwt分解得到的低頻分量{l1,l2,l3}和步驟4中變化特征提取得到的高頻分量進(jìn)行重構(gòu)，其中，步驟5具體包括：采用與步驟3中相同的兩組并行的基于正交變換的q-shift(quartersampleshift)實數(shù)濾波器組對步驟3中dt-cwt分解得到的低頻{l1,l2,l3}和步驟4經(jīng)過變化特征提取的高頻分量進(jìn)行重構(gòu)(即逆變換)。通過對l1和進(jìn)行逆變換得重構(gòu)后的第一層結(jié)果i1，通過對l2和進(jìn)行逆變換得重構(gòu)后的第一層結(jié)果i2，通過對l3和進(jìn)行逆變換得重構(gòu)后的第一層結(jié)果i3。

優(yōu)選的，步驟6具體包括以下子步驟：

步驟6.1，根據(jù)v1＝(min(ia(i,j))+level)/2、v2＝(level+max(ia(i,j))/2，對重構(gòu)后的結(jié)果i1，i2，i3分別求初始聚類中心v1、v2，記vm＝[v1,v1]，其中m表示迭代次數(shù)，level為自適應(yīng)bayes閾值，a＝1,2,3，max(ia(i,j))和min(ia(i,j))分別表示i1，i2或i3的像素最大和最小值，1≤i≤p，1≤j≤q的整數(shù)，p為ia(i,j)對應(yīng)的行數(shù)，q為ia(i,j)對應(yīng)的列數(shù)；

步驟6.2，根據(jù)式

對i1，i2和i3計算隸屬度矩陣u，其中c為聚類中心數(shù)，n為權(quán)值，||hj-cj||為樣本hj到第i類聚類中心的距離度量，||hj-ck||為樣本hj到第k類聚類中心的距離度量；

步驟6.3，根據(jù)步驟6.2中得到的隸屬度矩陣u，分別求初始標(biāo)記場并更新步驟6.2的聚類中心vm，m表示迭代次數(shù)；

步驟6.4，判斷迭代終止條件(優(yōu)選的，迭代終止條件設(shè)置為15次)，重復(fù)步驟6.2、6.3不斷對標(biāo)記場進(jìn)行迭代更新，得到最終迭代的標(biāo)記場，也即是初始分割結(jié)果w1^m、m表示迭代次數(shù)；

步驟6.5，對初始分割結(jié)果w1^m、由和分別求均值μm和方差即為特征場的參數(shù)，其中a＝1,2,3，p為輸入圖像的行數(shù)，q為輸入圖像的列數(shù)，i、j分別為1≤i≤p，1≤j≤q的整數(shù)，m表示迭代次數(shù)；

步驟6.6，采用條件迭代模型icm算法，由能量函數(shù)(只考慮二階鄰域系統(tǒng)，其中c為勢團(tuán)，c為二階勢團(tuán)集合，vc為勢函數(shù)，(w1,w2)表示二階鄰域系統(tǒng)，m為迭代次數(shù))計算當(dāng)前分割每個像素的標(biāo)記場能量um(w)，由能量函數(shù)計算當(dāng)前分割每個像素的特征場能量um(w,f)；

步驟6.7，根據(jù)能量最小原則em＝min(um(w)+um(w,f))，即根據(jù)迭代終止后特征場和標(biāo)記場總能量之和最小值所在變化或未變化類來更新當(dāng)前分割結(jié)果，更新后的分割結(jié)果記為m表示迭代次數(shù)；

步驟6.8，判斷循環(huán)終止條件(優(yōu)選的，迭代終止條件設(shè)置為30次)，重復(fù)步驟6.5、6.6和6.7不斷對當(dāng)前分割進(jìn)行迭代更新，達(dá)到最終的mrf-map最優(yōu)分割結(jié)果也即是三層變化檢測結(jié)果，并記為z1，z2，z3，其中m表示迭代次數(shù)。

優(yōu)選的，步驟7具體包括以下子步驟：

步驟7.1，統(tǒng)計比較通過步驟6獲得的變化檢測結(jié)果z2、z3與z1中對應(yīng)位置像素值是否相同，①若z1(i,j)≠z2(i,j)且z1(i,j)≠z3(i,j)，則記u＝0；②若z1(i,j)＝z2(i,j)≠z3(i,j)或z1(i,j)＝z3(i,j)≠z2(i,j)，則記u＝0.5；③若z1(i,j)＝z2(i,j)＝z3(i,j)，則記u＝1，其中z1(i,j)、z2(i,j)、z3(i,j)表示z1、z2、z3對應(yīng)(i,j)位置的像素值，u表示相似度；

步驟7.2，判斷相似度u是否滿足u≥0.5，若滿足，保持z1(i,j)值不變，若不滿足，取z2像素值z2(i,j)對z1像素值z1(i,j)進(jìn)行更新，更新后的結(jié)果即為融合后的變化檢測結(jié)果。

本發(fā)明還提供了一種結(jié)合dt-cwt和mrf的遙感圖像變化檢測的裝置，所述裝置包括計算機(jī)存儲介質(zhì)和處理器，其中，所述計算機(jī)存儲介質(zhì)中存有計算機(jī)程序，所述處理器被配置為執(zhí)行所述計算機(jī)程序，以執(zhí)行根據(jù)上述任一實施例所述的方法。

本發(fā)明有助于圖像信息的多方向表達(dá)和多尺度分析，充分利用了像素間的相關(guān)性，不僅能有效平衡高頻信息的保留和噪聲的去除，而且邊緣檢測更加平滑，變化檢測結(jié)果具有較好的區(qū)域一致性，大大降低了誤檢率和去除配準(zhǔn)誤差的影響，極大提高了遙感圖像變化檢測的精度。

附圖說明

圖1為根據(jù)本發(fā)明實施例的、結(jié)合dt-cwt和mrf的遙感圖像變化檢測方法的流程圖。

圖2為根據(jù)本發(fā)明實施例的、3層dt-cwt分解樹狀圖。

圖3為根據(jù)本發(fā)明實施例的、高頻變化特征提取流程圖。

圖4為根據(jù)本發(fā)明實施例的、基于自適應(yīng)bayes閾值的fcm分割流程圖。

圖5為根據(jù)本發(fā)明實施例的、icm分割流程圖。

圖6為根據(jù)本發(fā)明實施例的、基于相似度的變化檢測結(jié)果融合流程圖。

具體實施方式

下面對本發(fā)明進(jìn)行更全面的描述，其中說明本發(fā)明的示例性實施例。

如圖1所示，本發(fā)明提出的結(jié)合雙樹復(fù)小波變換(dual-treecomplexwavelettransform，dt-cwt)和馬爾科夫隨機(jī)場(mrf)的遙感圖像變化檢測方法包括以下步驟：

步驟1，輸入大小相同的經(jīng)過精確配準(zhǔn)的兩期同區(qū)域的遙感影像對應(yīng)的兩幅灰度圖像x1、x2；

其中，精確配準(zhǔn)可由軟件enviclassic下registration模塊并通過選擇同名控制點(diǎn)來完成，配準(zhǔn)精度小于一個像素，本發(fā)明進(jìn)行的變化檢測是基于像素的兩期的遙感圖像變化檢測，故在進(jìn)行變化檢測之前要對圖像進(jìn)行精確的配準(zhǔn)，避免因配準(zhǔn)引起的誤檢率上升。根據(jù)本發(fā)明的優(yōu)選實施例，所述兩幅灰度圖像x1、x2的大小均為p×q，其中，p為所述輸入圖像(即，所述灰度圖像)的行數(shù)，q為所述輸入圖像(即，所述灰度圖像)的列數(shù)。

步驟2，求所述兩幅灰度圖像x1、x2的差值圖像d；

具體的，采用差值法，根據(jù)公式d(i,j)＝|x2(i,j)-x1(i,j)|求所述差值圖像d，其中，d(i,j)為差值圖像d在位置(i，j)處的像素值，x1(i,j)和x2(i,j)分別是灰度圖像x1和x2的在位置(i，j)處的像素值，i，j分別為1≤i≤p，1≤j≤q的整數(shù)，p為所述輸入圖像的行數(shù)，q為所述輸入圖像的列數(shù)。

發(fā)明人利用差值法、比值法、對數(shù)法求取差值圖像進(jìn)行了對比實驗，發(fā)現(xiàn)差值法求取差異圖像效果較好，因此，本發(fā)明優(yōu)選使用了差值法，通過求取差值圖像來為后續(xù)變化檢測提供數(shù)據(jù)準(zhǔn)備。

步驟3，對所述差值圖像進(jìn)行dt-cwt分解；

步驟3具體包括：采用兩組并行的基于正交變換的q-shift(quartersampleshift)實數(shù)濾波器組來對步驟2中得到的差值圖像d進(jìn)行3層dt-cwt分解。第一層分解采用相差一個單位延遲的奇數(shù)長度的濾波器組，分解結(jié)果的一樹作為第一層分解復(fù)小波系數(shù)的實部，另一樹作為第一層分解復(fù)小波系數(shù)的虛部；第二、三層分解均采用遲近似為1/4采樣間隔的偶數(shù)長度的濾波器組，分解結(jié)果一樹作為第二、三層分解復(fù)小波系數(shù)的實部，另一樹作為第二、三層分解復(fù)小波系數(shù)的虛部。

圖2示出了3層dt-cwt分解樹狀圖，x為原始圖像，x(n)(n＝1，2，3)為第n層分解的高、低頻子帶，其中，陰影部分為低頻子帶，空白部分為六個方向的高頻子帶(±15°,±45°,±75°)。

dt-cwt的分解原理具體可描述為：第一層分解兩樹(a樹和b樹)均使用相差一個單位延遲的奇數(shù)長度的濾波器，在高層分解時使用群延遲近似為1/4采樣(優(yōu)選的采用群延遲近為1/4采樣)的偶數(shù)濾波器組，并且設(shè)定a樹1/4時延，b樹為3/4時延，以便保證兩樹之間所要求的1/2采樣周期的時延，這樣做的好處是：b樹的采樣點(diǎn)剛好位于a樹中因二抽取隔點(diǎn)采樣而舍棄的采樣值的點(diǎn)位(這樣就等價于沒進(jìn)行二抽取，而且非常易于實現(xiàn))。對于上述3層dt-cwt的分解，具體可以包括：分解的第一層兩樹均使用相差一個單位延遲的奇數(shù)長度的雙正交濾波器，在第二、三層分解時使用群延遲近似為1/4采樣(優(yōu)選的采用群延遲近為1/4采樣)的偶數(shù)濾波器組，并且設(shè)定一樹為1/4時延，一樹為3/4時延，以便保證兩樹之間所要求的1/2采樣周期的時延。

另外，dt-cwt是有一定冗余的，對于一維信號，產(chǎn)生了2:1的冗余度，對于n維信號，將產(chǎn)生2n:1的冗余度。為了克服奇/偶(odd/even)濾波器因嚴(yán)格線性相位引起的采樣結(jié)構(gòu)差異的缺點(diǎn)，本發(fā)明采用基于正交變換的q-shift(quartersampleshift)濾波器組對圖像進(jìn)行二維dt-cwt分解，每層分解可以得到一個低頻子帶和6個方向(±15°,±45°,±75°)的高頻子帶，所述低頻子帶記為{l1,l2,l3}，所述6個方向的高頻子帶記為hk,l，(k＝1,2,3；l＝±15°,±45°,±75°)，其中，k表示分解層數(shù)，l表示子帶的方向。由于dt-cwt不僅保持了傳統(tǒng)離散小波變換dwt(discretewavelettransform)良好的時頻分析能力，而且還具有近似的平移不變性、良好的方向選擇性、有限的數(shù)據(jù)冗余等特點(diǎn)，解決了復(fù)小波不能完全重構(gòu)的問題，故能夠較好的表達(dá)圖像的細(xì)節(jié)信息。

步驟4，利用基于fcm的mrf分割算法，對經(jīng)dt-cwt分解后得到的高頻子帶進(jìn)行變化特征提取；

具體來說，步驟4是對步驟3中經(jīng)過3層dt-cwt分解得到的高頻子帶hk,l，(k＝1,2,3；l＝±15°,±45°,±75°)進(jìn)行變化特征提取，具體實現(xiàn)方法見圖3詳解。經(jīng)過多尺度分解后的遙感圖像的低頻分量包含圖像的主要信息，高頻分量包含主要的突變信息(包括細(xì)節(jié)信息和噪聲)，由于噪聲絕大部分存在于分解的第一層高頻中，且第二層和第三層高頻信息同樣包含了圖像的大部分高頻信息，本發(fā)明通過舍去第一層高頻信息并利用mrf模型分割算法預(yù)先提取出其他層高頻域的變化特征，來達(dá)到在去除噪聲的同時較好地保留高頻細(xì)節(jié)信息的目的；其中，如圖3所示，步驟4具體包括以下步驟：

步驟4.1，對所述經(jīng)過3層dt-cwt分解得到的高頻子帶hk,l，(k＝1,2,3；l＝±15°,±45°,±75°)中的第一層高頻分量h1,l(l＝±15°,±45°,±75°)置零，作為第一層高頻特征提取結(jié)果其中，k表示分解層數(shù)，l表示子帶的方向；

步驟4.2，利用基于bayes閾值的fcm算法(具體實現(xiàn)方法見圖4詳解)對所述經(jīng)過3層dt-cwt分解得到的高頻子帶hk,l，(k＝1,2,3；l＝±15°,±45°,±75°)中的第二、三層高頻分量hk,l(k＝2,3；l＝±15°,±45°,±75°)進(jìn)行初始聚類分割，得到初始分割結(jié)果其中m表示迭代次數(shù)；

如圖4所示，步驟4.2具體包括以下子步驟：

步驟4.2.1，根據(jù)最小錯誤概率準(zhǔn)則，求自適應(yīng)bayes閾值level；

步驟4.2.2，根據(jù)v1＝(min(hk,l(i,j))+level)/2、v2＝(level+max(hk,l(i,j))/2，對所述經(jīng)過3層dt-cwt分解得到的第二、三層高頻分量hk,l，(k＝2,3；l＝±15°,±45°,±75°)分別求初始聚類中心v1、v2，記vm＝[v1,v1],其中m為迭代次數(shù)，level為自適應(yīng)bayes閾值，max(hk,l(i,j))和min(hk,l(i,j))分別表示每個子帶hk,l，(k＝2,3；l＝±15°,±45°,±75°)中的像素最大和最小值；

步驟4.2.3，根據(jù)式

對每個子帶hk,l(k＝2,3；l＝±15°,±45°,±75°)計算隸屬度矩陣u，其中，c為聚類中心數(shù)，n為權(quán)值，||hj-cj||為樣本hj到第i類聚類中心的距離度量，||hj-ck||為樣本hj到第k類聚類中心的距離度量，m表示迭代次數(shù)；

步驟4.2.4，根據(jù)步驟4.2.3中得到的隸屬度矩陣u，求初始標(biāo)記場并更新步驟4.2.2的聚類中心vm,其中m為迭代次數(shù)；

步驟4.2.5，重復(fù)步驟4.2.3、4.2.4不斷對標(biāo)記場進(jìn)行迭代更新，直到達(dá)到迭代終止條件(優(yōu)選的，迭代終止條件為15次)；

步驟4.2.6，輸出最終迭代的標(biāo)記場，也即為初始分割結(jié)果

步驟4.3，基于mrf與grf的等價性，利用條件迭代模型(icm)算法對4.2.6中的初始分割結(jié)果進(jìn)行迭代更新(優(yōu)選的，icm最大迭代次數(shù)設(shè)置為30次)，得到最終的最優(yōu)分割，記為

如圖5所示，步驟4.3具體包括以下子步驟：

步驟4.3.1，根據(jù)步驟4.2的fcm分割結(jié)果計算特征場參數(shù)：均值μm和方差

具體為，對基于bayes閾值的fcm分割結(jié)果由和分別求特征場參數(shù)μm和其中，m為迭代次數(shù)，μm為每個子帶變化與未變化類的均值，為每個子帶變化與未變化類的方差，p為所述輸入圖像的行數(shù)，q為所述輸入圖像的列數(shù)，i、j分別為1≤i≤p，1≤j≤q的整數(shù)；

步驟4.3.2，利用步驟4.3.1計算的特征場參數(shù)μm和來計算當(dāng)前分割的每個像素的特征場能量，并計算當(dāng)前分割的每個像素的標(biāo)記場能量；

具體為，對基于bayes閾值的fcm分割結(jié)果采用條件迭代模型icm算法，由能量函數(shù)(只考慮二階鄰域系統(tǒng)，其中，c為勢團(tuán)，c為二階勢團(tuán)集合，vc為勢函數(shù)，(w1,w2)表示二階鄰域系統(tǒng))計算當(dāng)前分割的每個像素的標(biāo)記場能量um(w)，由能量函數(shù)(其中，s為所有像素集合，f為特征場，m為迭代次數(shù))計算當(dāng)前分割的每個像素的特征場能量um(w,f)；

步驟4.3.3，根據(jù)能量最小原則來更新當(dāng)前分割結(jié)果，更新后的結(jié)果記為其中m為迭代次數(shù)；

具體的，根據(jù)能量最小原則em＝min(um(w)+um(w,f))，即,根據(jù)特征場和標(biāo)記場總能量之和最小值所在變化或未變化類來更新當(dāng)前分割結(jié)果，并記為m為迭代次數(shù)；

步驟4.3.4，重復(fù)步驟4.3.1-4.3.3，不斷對當(dāng)前分割進(jìn)行迭代更新，直到達(dá)到循環(huán)終止條件(迭代終止條件設(shè)置為30次)，獲得最終的mrf-map最優(yōu)分割結(jié)果m為迭代次數(shù)。

步驟4.3.5，輸出所述最終的mrf-map最優(yōu)分割結(jié)果m為迭代次數(shù)。

步驟4.4，對通過步驟4.3得到的最終分割結(jié)果進(jìn)行二進(jìn)制掩膜，得到掩膜圖bk,l(k＝2,3；l＝±15°,±45°,±75°)，其中m為迭代次數(shù)；

步驟4.5，將所述掩膜圖bk,l(k＝2,3；l＝±15°,±45°,±75°)矩陣像素值為1的記為變化類c，值為0的記為未變化類u；

步驟4.6，將c類像素對應(yīng)位置的二、三層高頻子帶hk,l(k＝2,3；l＝±15°,±45°,±75°)的像素值保持不變，u類像素對應(yīng)位置的二、三層高頻子帶的像素值置零，得到提取特征后的二、三層高頻分量

步驟5，對通過步驟3的dt-cwt分解得到的低頻分量{l1,l2,l3}和通過步驟4的變化特征提取得到的高頻分量進(jìn)行重構(gòu)；

具體的，采用與步驟3中相同的兩組并行的基于正交變換的q-shift(quartersampleshift)實數(shù)濾波器組，對步驟3中dt-cwt分解得到的低頻{l1,l2,l3}和圖3步驟1，圖4步驟6經(jīng)過變化特征提取的高頻分量進(jìn)行重構(gòu)(即逆變換)得到重構(gòu)后的三層結(jié)果i1，i2，i3，這里運(yùn)用重構(gòu)的思想來使得各層分量具有和原圖同樣地大小，避免了因尺度間插值引起的誤差，為后續(xù)進(jìn)行不同尺度間各分量的融合做準(zhǔn)備。

步驟6，利用所述基于fcm的mrf分割算法，對通過步驟5重構(gòu)后的三層結(jié)果i1，i2，i3進(jìn)行變化檢測；

本發(fā)明選擇基于dt-cwt與mrf模型相結(jié)合的變化檢測方法(同步驟4中的基于fcm初始分割的icm分割方法，即對重構(gòu)后的三層結(jié)果，先進(jìn)行基于自適應(yīng)bayes閾值的fcm初始分割，再進(jìn)行基于fcm的icm最終分割)，利用dt-cwt的多方向表達(dá)性、各向異性和多尺度特性，增強(qiáng)圖像信息的表達(dá)與分析，提高了變化檢測的精度，利用mrf能對像素間的相關(guān)性較好地描述，來降低檢測的誤檢率和消除配準(zhǔn)誤差的影響，提高魯棒性。

其中，步驟6具體包括：

步驟6.1，根據(jù)v1＝(min(ia(i,j))+level)/2、v2＝(level+max(ia(i,j))/2，分別對重構(gòu)后的結(jié)果i1，i2，i3求初始聚類中心v1、v2，記vm＝[v1,v1],其中m為迭代次數(shù)，level為自適應(yīng)bayes閾值，a＝1,2,3，max(ia(i,j))和min(ia(i,j))分別表示i1，i2或i3的像素最大和最小值，i、j分別為1≤i≤p，1≤j≤q的整數(shù)，p為ia(i,j)的行數(shù)，q為ia(i,j)的列數(shù)；

步驟6.2，根據(jù)式

對i1，i2和i3分別計算隸屬度矩陣u，其中c為聚類中心數(shù)，n為權(quán)值，||hj-cj||為樣本hj到第i類聚類中心的距離度量，||hj-ck||為樣本hj到第k類聚類中心的距離度量,m表示迭代次數(shù)；

步驟6.3，根據(jù)步驟6.2中得到的隸屬度矩陣u，分別求初始標(biāo)記場并更新步驟6.2的聚類中心vm，m表示迭代次數(shù)；

步驟6.4，判斷迭代終止條件(迭代終止條件設(shè)置為15次)，重復(fù)步驟6.2、6.3不斷對標(biāo)記場進(jìn)行迭代更新，得到最終迭代的標(biāo)記場，也即是初始分割結(jié)果w1^m、m表示迭代次數(shù)；

步驟6.5，對初始分割結(jié)果w1^m、由和分別求均值μm和方差即為特征場參數(shù)，其中a＝1,2,3；p為輸入圖像的行數(shù)，q為輸入圖像的列數(shù)，i、j分別為1≤i≤p，1≤j≤q的整數(shù),m表示迭代次數(shù)；

步驟6.6，采用條件迭代模型icm算法，由能量函數(shù)(只考慮二階鄰域系統(tǒng)，其中c為勢團(tuán)，c為二階勢團(tuán)集合，vc為勢函數(shù)，(w1,w2)表示二階鄰域系統(tǒng))計算當(dāng)前分割每個像素的標(biāo)記場能量um(w)，由能量函數(shù)(其中s為所有像素集合，s屬于s，f為特征場，m為迭代次數(shù))計算當(dāng)前分割每個像素的特征場能量um(w,f)；

步驟6.7，根據(jù)能量最小原則em＝min(um(w)+um(w,f))，即根據(jù)迭代終止后特征場和標(biāo)記場總能量之和最小值所在變化或未變化類來更新當(dāng)前分割結(jié)果，更新后的分割結(jié)果記為m表示迭代次數(shù)；

步驟6.8，判斷循環(huán)終止條件(迭代終止條件設(shè)置為30次)，重復(fù)步驟6.5、6.6和6.7不斷對當(dāng)前分割進(jìn)行迭代更新，達(dá)到最終的mrf-map最優(yōu)分割結(jié)果也即是三層變化檢測結(jié)果，并記為z1，z2，z3，其中m表示迭代次數(shù)。

步驟7，基于相似度融合規(guī)則，對通過步驟6得到的檢測結(jié)果進(jìn)行融合；

具體的，基于相似度融合規(guī)則，對步驟6中的三層變化檢測結(jié)果z1，z2，z3進(jìn)行融合，來進(jìn)一步提高變化檢測精度，經(jīng)試驗顯示，融合后的結(jié)果邊緣檢測更加平滑細(xì)膩，檢測雜點(diǎn)較少，具有較高的變化檢測精度和較強(qiáng)的魯棒性；其中，步驟7具體包括：

圖6示出了基于相似度的變化檢測結(jié)果融合過程，包括：

步驟7.1，統(tǒng)計比較通過步驟6獲得的變化檢測結(jié)果z2、z3與z1中對應(yīng)位置像素值是否相同，①若z1(i,j)≠z2(i,j)且z1(i,j)≠z3(i,j)，則記u＝0；②若z1(i,j)＝z2(i,j)≠z3(i,j)或z1(i,j)＝z3(i,j)≠z2(i,j)，則記u＝0.5；③若z1(i,j)＝z2(i,j)＝z3(i,j)，則記u＝1，其中u表示相似度；

步驟7.2，判斷相似度u是否滿足u≥0.5，若滿足，保持z1(i,j)值不變，若不滿足，取z2像素值z2(i,j)對z1像素值z1(i,j)進(jìn)行更新，更新后的結(jié)果即為融合后的變化檢測結(jié)果。

步驟8，輸出融合后的變化檢測結(jié)果，即為最終的變化檢測結(jié)果。

本發(fā)明還提出了一種結(jié)合dt-cwt和mrf的遙感圖像變化檢測的裝置，所述裝置包括計算機(jī)存儲介質(zhì)和處理器，其中，所述計算機(jī)存儲介質(zhì)中存有計算機(jī)程序，所述處理器被配置為執(zhí)行所述計算機(jī)程序，以執(zhí)行根據(jù)上述任一實施例所述的方法。

以上內(nèi)容僅為本發(fā)明的較佳實施例，對于本領(lǐng)域的普通技術(shù)人員，依據(jù)本發(fā)明的思想，在具體實施方式及應(yīng)用范圍上均會有改變之處，本說明書內(nèi)容不應(yīng)理解為對本發(fā)明的限制。