專利名稱：基于Treelet曲波域去噪的遙感圖像變化檢測(cè)方法
技術(shù)領(lǐng)域：
本發(fā)明屬于圖像處理技術(shù)領(lǐng)域，更進(jìn)一步涉及一種基于Treelet曲波域去噪的遙感圖像變化檢測(cè)方法。該方法可應(yīng)用于自然災(zāi)害中的災(zāi)情監(jiān)測(cè)和評(píng)估、森 林覆蓋率的監(jiān)測(cè)和評(píng)估、城市規(guī)劃等領(lǐng)域。
背景技術(shù)：
變化檢測(cè)是通過(guò)觀察不同時(shí)刻的一個(gè)物體或現(xiàn)象，鑒定其狀態(tài)差異的過(guò)程。隨著遙感技術(shù)和信息技術(shù)的發(fā)展，遙感圖像變化檢測(cè)已經(jīng)成為當(dāng)前遙感圖像分析研究的一個(gè)重要方向。遙感圖像變化檢測(cè)方法可以分為兩類先分類后比較法和先比較后分類法。先分類后比較法的優(yōu)點(diǎn)在于可以克服由于大氣、傳感器、季節(jié)、地面環(huán)境、多時(shí)相遙感圖像的分辨率等因素帶來(lái)的不同時(shí)相圖像間的差異。但該類方法存在誤差累計(jì)問(wèn)題，并且分類本身存在病態(tài)分割問(wèn)題，因此會(huì)導(dǎo)致最終的變化檢測(cè)精度不高。先比較后分類法的優(yōu)勢(shì)在于直觀、簡(jiǎn)單、易行，沒(méi)有先分類后比較法所存在的誤差累計(jì)問(wèn)題，是現(xiàn)在普遍應(yīng)用的一類多時(shí)相遙感圖像變化檢測(cè)方法?；诙喑叨葞缀畏治龅倪b感圖像變化檢測(cè)方法能夠在不同尺度上進(jìn)行檢測(cè)，并通過(guò)最優(yōu)尺度或者多尺度的融合確定最終的變化檢測(cè)結(jié)果，在一定程度上避免了單一尺度下對(duì)邊緣信息的保持和變化區(qū)域信息的檢測(cè)難以兼顧的缺點(diǎn)。西安電子科技大學(xué)在其專利申請(qǐng)“基于多尺度積和主成分分析的SAR圖像變化檢測(cè)方法”(專利申請(qǐng)?zhí)?200910023637. 1，公開(kāi)號(hào)CN101634709)中提出了一種多尺度積去噪和主成分分析融合的SAR圖像變化檢測(cè)方法。該方法雖然能夠減弱圖像誤配準(zhǔn)的影響， 但仍存在的不足是，該方法僅利用當(dāng)前尺度層和其下一個(gè)尺度層進(jìn)行去噪，沒(méi)有考慮當(dāng)前尺度層和其上一個(gè)尺度層的關(guān)系，使得該方法的去噪結(jié)果不夠理想，對(duì)噪聲信號(hào)較為敏感。 此外，由于SAR圖像并非完全服從廣義高斯模型，因此基于廣義高斯模型的最小錯(cuò)誤率閾值法會(huì)使得檢測(cè)結(jié)果中存在較多的偽變化信息，降低了遙感圖像變化檢測(cè)的精度。2010年黃世奇等在文獻(xiàn)“基于小波變換的多時(shí)相SAR圖像變化檢測(cè)技術(shù)”(測(cè)繪學(xué)報(bào)，2010，39(2) 180-186.)中提出了一種基于可靠尺度融合的遙感圖像變化檢測(cè)方法。 該方法通過(guò)局部變化系數(shù)和全局變化系數(shù)的關(guān)系確定可靠尺度，并采用加入權(quán)重的可靠尺度特征級(jí)融合策略進(jìn)行變化區(qū)域檢測(cè)。該方法雖然能夠?qū)ψ兓瘏^(qū)域類進(jìn)行了進(jìn)一步的類型分析，得到變化區(qū)域減弱類和變化區(qū)域增強(qiáng)類，但仍存在的不足是，該方法在進(jìn)行可靠尺度特征級(jí)融合時(shí)僅考慮了圖像的低頻信息，因此變化區(qū)域的邊緣信息不能夠得到較好的保持。此外，由于該方法使用EM算法進(jìn)行分類，使得檢測(cè)結(jié)果中存在較多的偽變化信息，影響了遙感圖像變化檢測(cè)的精度。

發(fā)明內(nèi)容
本發(fā)明針對(duì)上述現(xiàn)有技術(shù)存在的不足，提出了一種基于Treelet曲波域去噪的遙感圖像變化檢測(cè)方法。本發(fā)明對(duì)噪聲信號(hào)具有較好的魯棒性，能夠較好的保持變化區(qū)域的邊緣信息，減少偽變化信息，具有較高的變化檢測(cè)精度。本發(fā)明實(shí)現(xiàn)上述目的的思路是在對(duì)讀入的數(shù)據(jù)進(jìn)行中值濾波、構(gòu)造差異圖像、快速離散曲波分解和曲波變換系數(shù)分類處理之后，使用Treelet變換自適應(yīng)的對(duì)Detail尺度層進(jìn)行去噪，再利用曲波域低頻變化圖像及其邊緣圖計(jì)算變化比例閾值，最后進(jìn)行圖像分類和圖像融合。本發(fā)明的具體步驟包括如下 (1)讀入同一地區(qū)不同時(shí)刻獲取的兩幅遙感圖像。(2)中值濾波2a)確定正方形窗口 選取步驟(1)中的一幅遙感圖像，以該圖像中的某一像素點(diǎn)為中心，選取一個(gè)大小為NuXNu的正方形窗口，其中，Nu為奇數(shù)；2b)確定濾波值將正方形窗口中全部像素點(diǎn)的灰度值按照由大到小的順序排列，組成一個(gè)灰度序列，選取位于灰度序列中間位置的灰度值作為濾波值；2c)濾波用濾波值替代步驟2a)中像素點(diǎn)的灰度值；2d)重復(fù)步驟2a)至步驟2c)，直至處理完圖像中的全部像素點(diǎn)；2e)按照步驟2a)至步驟2d)，對(duì)步驟(1)中的另外一幅遙感圖像進(jìn)行處理，得到濾波后的兩幅圖像。(3)構(gòu)造對(duì)數(shù)差異圖像3a)對(duì)步驟(2)中濾波后的兩幅圖像分別進(jìn)行對(duì)數(shù)運(yùn)算，得到兩幅對(duì)數(shù)圖像；3b)對(duì)兩幅對(duì)數(shù)圖像對(duì)應(yīng)像素點(diǎn)進(jìn)行減法運(yùn)算，并對(duì)減法運(yùn)算的結(jié)果取絕對(duì)值，得到一幅對(duì)數(shù)差異圖像。(4)構(gòu)造絕對(duì)值差異圖像對(duì)步驟(2)中濾波后的兩幅圖像對(duì)應(yīng)像素點(diǎn)進(jìn)行減法運(yùn)算，并對(duì)減法運(yùn)算的結(jié)果取絕對(duì)值，得到一幅絕對(duì)值差異圖像。(5)對(duì)步驟(3)中的對(duì)數(shù)差異圖像進(jìn)行Sc層快速離散曲波分解。(6)根據(jù)曲波(Curvelet)變換域頻率的大小，對(duì)步驟(5)中的曲波變換系數(shù)進(jìn)行分類。(7)將Fine尺度層的曲波變換系數(shù)全部賦值為零。(8)Detail尺度層去噪8a)水平方向去噪對(duì)Detail尺度層的任一方向子帶，應(yīng)用Treelet變換自適應(yīng)的計(jì)算水平方向去噪閾值，并進(jìn)行軟閾值去噪，得到該方向子帶的水平方向去噪結(jié)果；8b)垂直方向去噪對(duì)步驟8a)中選取的方向子帶，應(yīng)用Treelet變換自適應(yīng)的計(jì)算垂直方向去噪閾值，并進(jìn)行軟閾值去噪，得到該方向子帶的垂直方向去噪結(jié)果；8c)平均加權(quán)求和將水平方向去噪結(jié)果和垂直方向去噪結(jié)果進(jìn)行平均加權(quán)求和，得到該方向子帶的最終去噪結(jié)果；8d)重復(fù)步驟8a)、8b)和8c)，直至處理完Detail尺度層的所有方向子帶。(9)曲波逆變換對(duì)Coarse尺度層的Curvelet變換系數(shù)、經(jīng)步驟(7)處理后的Fine尺度層的 Curvelet變換系數(shù)和經(jīng)步驟(8)處理后的Detail尺度層的Curvelet變換系數(shù)進(jìn)行逆Curvelet變換，得到重構(gòu)圖像。(10)計(jì) 算變化比例閾值IOa)低頻差異圖像將Fine尺度層和Detail尺度層的Curvelet變換系數(shù)全部置零，并與Coarse尺度層的Curvelet變換系數(shù)進(jìn)行逆Curvelet變換，得到僅包含低頻信息的低頻差異圖像；IOb)低頻變化圖像對(duì)低頻差異圖像采用最大類間方差(Otsu)法分類，得到低頻變化圖像；IOc)最終邊緣圖分別采用Canny、S0bel、PreWitt邊緣算子提取低頻變化圖像的邊緣，得到三幅邊緣圖，通過(guò)邏輯或運(yùn)算將此三幅邊緣圖合成一幅最終邊緣圖；IOd)計(jì)算變化比例閾值統(tǒng)計(jì)最終邊緣圖中的邊緣像素點(diǎn)個(gè)數(shù)和低頻變化圖像中變化像素點(diǎn)的個(gè)數(shù)，計(jì)算變化比例閾值。(11)分類Ila)通過(guò)變化比例閾值對(duì)步驟(4)獲得的絕對(duì)值差異圖像進(jìn)行分類，得到一幅分類圖像；lib)通過(guò)變化比例閾值對(duì)步驟(9)獲得的重構(gòu)圖像進(jìn)行分類，得到另一幅分類圖像。(12)獲得變化檢測(cè)結(jié)果圖對(duì)步驟Ila)獲得的分類圖像和步驟lib)獲得的分類圖像進(jìn)行融合，得到變化檢測(cè)結(jié)果圖。本發(fā)明與現(xiàn)有技術(shù)相比具有以下優(yōu)點(diǎn)第一，本發(fā)明由于采用了 Treelet變換自適應(yīng)的對(duì)Detail尺度層進(jìn)行去噪，克服了現(xiàn)有技術(shù)對(duì)噪聲較為敏感的缺點(diǎn)，使得本發(fā)明的檢測(cè)結(jié)果不容易受噪聲信號(hào)的干擾，對(duì)噪聲具有較好的魯棒性。第二，本發(fā)明根據(jù)Curvelet域低頻變化圖像及其邊緣圖計(jì)算變化比例閾值，克服了現(xiàn)有技術(shù)存在較多偽變換信息的缺點(diǎn)，使得本發(fā)明的檢測(cè)精度得到了提高。第三，本發(fā)明通過(guò)圖像融合獲得變化檢測(cè)結(jié)果圖像，克服了現(xiàn)有技術(shù)對(duì)變化區(qū)域邊緣信息保持不理想的缺點(diǎn)，使得本發(fā)明的變化區(qū)域邊緣信息得到了更加準(zhǔn)確、全面的保持。


圖1為本發(fā)明的流程圖；圖2為本發(fā)明的仿真效果圖。
具體實(shí)施例方式下面結(jié)合附圖1對(duì)本發(fā)明的步驟做進(jìn)一步的詳細(xì)描述。步驟1，讀入同一地區(qū)不同時(shí)刻獲取的兩幅遙感圖像。步驟2，中值濾波。2a)確定正方形窗口 選取步驟1中的一幅遙感圖像，以該圖像中的某一像素點(diǎn)為中心，選取一個(gè)NuXNu的正方形窗口，其中，Nu為奇數(shù)，本發(fā)明實(shí)施例中選取一個(gè)3X3的正方形窗口。2b)確定濾波值將正方形窗口中全部像素點(diǎn)的灰度值按照由大到小的順序排列，組成一個(gè)灰度序列，選取位于灰度序列中間位置的灰度值作為濾波值。2c)濾波將濾波值替代步驟2a)中像素點(diǎn)的灰度值。2d)重復(fù)步驟2a)至步驟2c)，直至處理完圖像中的全部像素點(diǎn)。2e)按照步驟2a)至步驟2d)，對(duì)步驟(1)中的另一幅圖像進(jìn)行處理，得到濾波后的兩幅圖像。步驟3，構(gòu)造對(duì)數(shù)差異圖像。3a)對(duì)步驟2中濾波后的兩幅圖像分別進(jìn)行對(duì)數(shù)運(yùn)算，得到兩幅對(duì)數(shù)圖像。其對(duì)數(shù)運(yùn)算公式為I3 (m, n) = log (I1 (m, n)+l)I4 (m, η) = log (I2 (m, η)+1)其中，I3和I4為兩幅對(duì)數(shù)圖像(大小為MXN)，I1和I2為濾波后的兩幅圖像，m和 η為圖像的行序號(hào)和列序號(hào)，m = 0，1，2，…，M_l，η = 0，1，2，…，Ν-1。3b)對(duì)兩幅對(duì)數(shù)圖像對(duì)應(yīng)像素點(diǎn)進(jìn)行減法運(yùn)算，并對(duì)減法運(yùn)算的結(jié)果取絕對(duì)值，得到一幅對(duì)數(shù)差異圖像L。L (m, n) = 113 (m, η)-I4 (m, η) |其中，L為對(duì)數(shù)差異圖像，I3和I4為兩幅對(duì)數(shù)圖像(大小為MXN)，m和η為圖像的行序號(hào)和列序號(hào)，m = 0，1，2，…，Μ-1，η = 0，1，2，…，Ν-1。步驟4，構(gòu)造絕對(duì)值差異圖像。對(duì)步驟2中濾波后的兩幅圖像對(duì)應(yīng)像素點(diǎn)進(jìn)行減法運(yùn)算，并對(duì)減法運(yùn)算的結(jié)果取絕對(duì)值，得到一幅絕對(duì)值差異圖像Α。A (m, n) = 111 (m, η)-I2 (m, η) |其中，A為絕對(duì)值差異圖像，I1和I2為濾波后的兩幅圖像(大小為MXN)，m和η 為圖像的行序號(hào)和列序號(hào)，m = 0，1，2，-,M-Ln = 0,1,2,…，N-1。步驟5，快速離散曲波分解。對(duì)步驟3中的對(duì)數(shù)差異圖像L進(jìn)行Sc層快速離散曲波分解。
權(quán)利要求
1.基于Treelet曲波域去噪的遙感圖像變化檢測(cè)方法，包括如下步驟(1)讀入同一地區(qū)不同時(shí)刻獲取的兩幅遙感圖像；(2)中值濾波2a)確定正方形窗口 選取步驟(1)中的一幅遙感圖像，以該圖像中的某一像素點(diǎn)為中心，選取一個(gè)大小為NuXNu的正方形窗口，其中，Nu為奇數(shù)；2b)確定濾波值將正方形窗口中全部像素點(diǎn)的灰度值按照由大到小的順序排列，組成一個(gè)灰度序列，選取位于灰度序列中間位置的灰度值作為濾波值； 2c)濾波用濾波值替代步驟2a)中像素點(diǎn)的灰度值； 2d)重復(fù)步驟2a)至步驟2c)，直至處理完圖像中的全部像素點(diǎn)； 2e)按照步驟2a)至步驟2d)，對(duì)步驟(1)中的另外一幅遙感圖像進(jìn)行處理，得到濾波后的兩幅圖像；(3)構(gòu)造對(duì)數(shù)差異圖像3a)對(duì)步驟(2)中濾波后的兩幅圖像分別進(jìn)行對(duì)數(shù)運(yùn)算，得到兩幅對(duì)數(shù)圖像； 3b)對(duì)兩幅對(duì)數(shù)圖像對(duì)應(yīng)像素點(diǎn)進(jìn)行減法運(yùn)算，并對(duì)減法運(yùn)算的結(jié)果取絕對(duì)值，得到一幅對(duì)數(shù)差異圖像；(4)構(gòu)造絕對(duì)值差異圖像對(duì)步驟(2)中濾波后的兩幅圖像對(duì)應(yīng)像素點(diǎn)進(jìn)行減法運(yùn)算，并對(duì)減法運(yùn)算的結(jié)果取絕對(duì)值，得到一幅絕對(duì)值差異圖像；(5)對(duì)步驟(3)中的對(duì)數(shù)差異圖像進(jìn)行Sc層快速離散曲波分解；(6)根據(jù)曲波變換域頻率的大小，對(duì)步驟(5)中的Curvelet變換系數(shù)進(jìn)行分類；(7)將Fine尺度層的Curvelet變換系數(shù)全部賦值為零； (S)Detail尺度層去噪8a)水平方向去噪對(duì)Detail尺度層的任一方向子帶，應(yīng)用Treelet變換自適應(yīng)的計(jì)算水平方向去噪閾值，并進(jìn)行軟閾值去噪，得到該方向子帶的水平方向去噪結(jié)果；8b)垂直方向去噪對(duì)步驟8a)中選取的方向子帶，應(yīng)用Treelet變換自適應(yīng)的計(jì)算垂直方向去噪閾值，并進(jìn)行軟閾值去噪，得到該方向子帶的垂直方向去噪結(jié)果；8c)平均加權(quán)求和將水平方向去噪結(jié)果和垂直方向去噪結(jié)果進(jìn)行平均加權(quán)求和，得到該方向子帶的最終去噪結(jié)果；8d)重復(fù)步驟8a)、8b)和8c)，直至處理完Detail尺度層的所有方向子帶；(9)曲波逆變換對(duì)Coarse尺度層的Curvelet變換系數(shù)、經(jīng)步驟(7)處理后的Fine尺度層的Curvelet 變換系數(shù)和經(jīng)步驟(8)處理后的Detail尺度層的Curvelet變換系數(shù)進(jìn)行逆Curvelet變換，得到重構(gòu)圖像；(10)計(jì)算變化比例閾值IOa)低頻差異圖像將Fine尺度層和Detail尺度層的Curvelet變換系數(shù)全部置零， 并與Coarse尺度層的Curvelet變換系數(shù)進(jìn)行逆Curvelet變換，得到僅包含低頻信息的低頻差異圖像；IOb)低頻變化圖像對(duì)低頻差異圖像采用最大類間方差(Otsu)法分類，得到低頻變化圖像；IOc)最終邊緣圖分別采用Canny、Sobel, Prewitt邊緣算子提取低頻變化圖像的邊緣，得到三幅邊緣圖，通過(guò)邏輯或運(yùn)算將此三幅邊緣圖合成一幅最終邊緣圖；IOd)計(jì)算變化比例閾值統(tǒng)計(jì)最終邊緣圖中的邊緣像素點(diǎn)個(gè)數(shù)和低頻變化圖像中變化像素點(diǎn)的個(gè)數(shù)，計(jì)算變化比例閾值；(11)分類Ila)通過(guò)變化比例閾值對(duì)步驟(4)獲得的絕對(duì)值差異圖像進(jìn)行分類，得到一幅分類圖像；lib)通過(guò)變化比例閾值對(duì)步驟(9)獲得的重構(gòu)圖像進(jìn)行分類，得到另一幅分類圖像；(12)獲得變化檢測(cè)結(jié)果圖對(duì)步驟Ila)獲得的分類圖像和步驟lib)獲得的分類圖像進(jìn)行融合，得到變化檢測(cè)結(jié)果圖。
2.根據(jù)權(quán)利要求1所述的基于Treelet曲波域去噪的遙感圖像變化檢測(cè)方法，其特征在于步驟3a)所述的對(duì)數(shù)運(yùn)算公式為13(m, n) = Iogd1 (m, n)+l)14(m, η) = log (I2 (m, n)+l)其中，I3和I4為兩幅對(duì)數(shù)圖像(大小為MXN)，I1和I2為濾波后的兩幅圖像，m和η為圖像的行序號(hào)和列序號(hào)，m = 0，1，2，…，Μ-1，η = 0，1，2，…，Ν-1。
3.根據(jù)權(quán)利要求1所述的基于Treelet曲波域去噪的遙感圖像變化檢測(cè)方法，其特征在于步驟(5)所述的快速離散曲波分解的分解層數(shù)為5。
4.根據(jù)權(quán)利要求1所述的基于Treelet曲波域去噪的遙感圖像變化檢測(cè)方法，其特征在于步驟(6)所述的曲波變換系數(shù)分類的步驟如下將第Sc分解層的低頻系數(shù)歸類為 Coarse尺度層，第Sc-l、Sc-2、…、2分解層的高頻系數(shù)歸類為Detail尺度層，第1分解層的高頻系數(shù)歸類為Fine尺度層，其中，Sc為快速離散曲波分解的分解層數(shù)。
5.根據(jù)權(quán)利要求1所述的基于Treelet曲波域去噪的遙感圖像變化檢測(cè)方法，其特征在于步驟8a)所述的計(jì)算水平方向去噪閾值的步驟如下計(jì)算第1分解層第I1個(gè)尺度向量的歸一化能量
6.根據(jù)權(quán)利要求1所述的基于Treelet曲波域去噪的遙感圖像變化檢測(cè)方法，其特征在于步驟8c)所述的平均加權(quán)求和的計(jì)算公式為
7.根據(jù)權(quán)利要求1所述的基于Treelet曲波域去噪的遙感圖像變化檢測(cè)方法，其特征在于步驟(10)所述的變化比例閾值的計(jì)算公式如下P = P^P2其中，P為變化比例閾值，P1為最終邊緣圖中的邊緣像素點(diǎn)總個(gè)數(shù)，P2為低頻變化檢測(cè)圖像中的變化像素點(diǎn)總個(gè)數(shù)。
8.根據(jù)權(quán)利要求1所述的基于Treelet曲波域去噪的遙感圖像變化檢測(cè)方法，其特征在于步驟Ila)所述的分類的步驟如下將絕對(duì)值差異圖像中所有像素點(diǎn)的灰度值按照由大到小的順序排列，將前P個(gè)最大灰度值對(duì)應(yīng)的像素點(diǎn)視為變化類，賦值為1 ；其他像素點(diǎn)視為非變化類，賦值為0，得到分類圖像，其中，P為變化比例閾值。
9.根據(jù)權(quán)利要求1所述的基于Treelet曲波域去噪的遙感圖像變化檢測(cè)方法，其特征在于步驟(12)所述的融合方式如下Mp (m, n) = X1 (m, η) Π X2 (m, η)其中，Mp(m, η)為變化檢測(cè)結(jié)果圖，X1和X2為兩幅分類圖像(大小為MXN)，m和η為圖像的行序號(hào)和列序號(hào)，m = 0，1，2，-,M-Ln = 0,1,2,…，N-I，Π表示邏輯與運(yùn)算。
全文摘要
本發(fā)明公開(kāi)了一種基于Treelet曲波域去噪的遙感圖像變化檢測(cè)方法，其實(shí)現(xiàn)步驟為(1)讀入數(shù)據(jù)；(2)中值濾波；(3)構(gòu)造對(duì)數(shù)差異圖像；(4)構(gòu)造絕對(duì)值差異圖像；(5)快速離散曲波分解；(6)曲波變換系數(shù)分類；(7)Fine尺度層置零；(8)Detail尺度層去噪；(9)曲波逆變換；(10)計(jì)算變化比例閾值；(11)分類；(12)獲得變化檢測(cè)結(jié)果圖。本發(fā)明對(duì)噪聲具有較好的魯棒性，能夠較好的保持變化區(qū)域的邊緣信息，減少偽變化信息，具有較高的檢測(cè)精度，可用于災(zāi)情監(jiān)測(cè)、森林覆蓋率評(píng)估、城市規(guī)劃等領(lǐng)域。
文檔編號(hào)G06T7/00GK102360500SQ20111019215
公開(kāi)日2012年2月22日 申請(qǐng)日期2011年7月8日 優(yōu)先權(quán)日2011年7月8日
發(fā)明者萬(wàn)義萍, 侯彪, 公茂果, 張小華, 焦李成, 王桂婷, 王爽, 田小林, 鐘樺 申請(qǐng)人:西安電子科技大學(xué)