《測繪學報》
構建與學術的橋樑 拉近與權威的距離
虢英傑, 朱蘭豔
昆明理工大學國土資源工程學院, 雲南 昆明 650093
收稿日期：2016-11-16； 修回日期：2017-05-19
基金專案：國家自然科學基金（41261043）
第一作者簡介：虢英傑(1993—), 男, 碩士, 研究方向為現從事定量遙感、變化資訊識別和提取。 E-mail:884833945@qq.com
通信作者：朱蘭豔, E-mail：KMOY2222@sina.com
摘要：目前， 大多數高解析度遙感影像變化檢測方法的效率不高， 精度也難以達到實際應用要求， 本文引入頻域分析的方法進行變化檢測預處理， 達到提高變化資訊檢測精度和效率的目標。 首先， 對影像頻域能量的分佈特徵進行分析， 並以此為基礎， 利用累加距離匹配函數， 完成試驗影像的合理分塊；其次， 通過變化性紋理規則度的定義和計算， 實現對同能量下變化資訊的精確描述， 縮小變化區域搜索範圍， 有效抑制背景資訊對目標資訊的干擾。 最後， 通過3種常用的變化檢測方法的驗證，
A Preprocess for Change Detection Method from High Resolution Remotely Sensed Imagery in Frequency Domain
GUO Yingjie, ZHU Lanyan
Abstract: At present, the efficiency of change detection from high resolution remotely sensed imagery is generally not very high, and accuracy couldn't meet the practical requirements. A method which introduced frequency domain analysis is put forward in this paper. Firstly, it studies the distribution of image frequency energy. Using it as a basis, the image division is accomplished by summed-up distance matching function. Secondly, the definition and computation of texture regularity of change information(TRCI) clearly and accurately expresses the change information according to equal energy rule of image blocks, greatly narrowing the search area. Thus, the background information could be greatly suppressed so as to improve the efficiency and precision of detection. Finally, through measuring and comparing of the accuracy and efficiency of three different detection algorithms before and after introducing the preprocess, we find that it could provides more accurate results on the existed basis.
Key words: preprocess for change detection frequency domain analysis TRCI optimal image blocking
隨著社會與科技的發展， 人類開發資源與改造自然的能力不斷增強， 自然界的變化和人類的各種活動每天都在改變著地表景觀[1]。 人口的快速增長及城市化的發展， 更是加快了這種變化的速度[2]。 近年來， 隨著高解析度遙感影像的普及化使用和相關研究的不斷深入， 各種變化檢測方法逐漸被發掘和提出， 但由於在大範圍影像中， 不同尺度的地物光譜相互重疊， 降低了光譜域的統計可分性， 使同一地物在不同時相的影像中呈現光譜、空間特徵的差異， 嚴重影響變化檢測的效果， 從而限制了有較高自動化水準的高解析度遙感影像變化檢測方法在實際工程中的應用[3-5]。
根據以上分析， 本文嘗試從影像地物的視覺特徵出發， 在保證地表景觀完整性的前提下， 確定影像中主要地物的紋理週期性和方向性， 實現影像分塊。 以此為基礎， 結合子影像塊的變化性紋理規則度， 完成變化資訊所在影像塊的衡量和標識， 從而縮小變化區域的搜索範圍，
1 研究方法1.1 研究思路
本文以高空間解析度衛星圖像為試驗資料， 經過圖像抑噪等預處理後， 首先對圖像進行二維離散傅裡葉變換， 通過簡單的圖像變換和數學推導進行線狀特徵的頻譜分析， 證明線狀特徵和頻域變換後的譜線之間的方向關係；然後根據頻域能量的徑向分佈特徵和SDMF所表徵的週期性， 得到圖像分塊的分塊參數；最後， 結合論文所提出的變化性紋理規則度， 實現有效變化資訊的精確表示。 完成確定變化資訊所在區域的任務， 為檢測演算法在現有的檢測精度和效率的基礎上做進一步提高提供了保障，
圖 1 引入頻域分析的變化檢測預處理過程Fig. 1 The preprocess of change detection based on frequency spectrum analysis
圖選項
1.2 影像頻域能量的資訊表徵
為了定量分析遙感影像的空間頻率及其對應能量值所反映的地物紋理、細部結構和邊緣特徵， 通過二維離散傅裡葉變換對影像進行空域向頻域的轉換[8]。 設影像的大小為M×N， 其二維空域表達為f(x,y)， 對應的頻域表達為F(u,v)， 則二維離散傅裡葉變換的數學定義式為
(1)
根據傅裡葉變化位移定理， 設影像的變數x,y在二維空間記憶體在位移x-x0和y-y0， 那麼二維離散傅裡葉變化滿足以下關係
(2)
從式(2) 可以看出， 時(空)域信號的平移並不能改變其頻率和幅度成分， 發生改變的僅僅是相位譜部分。 由於空間域內單個條帶圖像必定帶有直流電平的屬性[9]， 故其所對應的譜線必然會通過頻域平面的中心點，並集中反映在垂直于原特徵線方向的譜線上，這一現象稱之為傅氏變換幅度譜的自配准性質[10]。
為了進一步探討影像中地物特徵與能量譜線之間的關係，需要將傅裡葉頻譜能量分析和統計方法結合起來考慮。本文採用文獻[11]中所提出的環狀採樣和楔狀採樣的方法，假設(r,θ)為(u,v)平面上的極座標，以採樣間隔為θi≤arctan(v/u)＜θi+1，在以頻譜圖中心為原點出發的扇面上對r進行求和，其關係為
(3)
則稱Er(θi)為影像頻譜能量的角向分佈特徵(楔特徵)，如圖 2(a)所示。
圖 2 能量譜的環狀採樣和楔狀採樣Fig. 2 Radius sampling and angle sampling of the frequency spectrum
圖選項
在高解析度遙感影像中，地物特徵線會隨機分佈在很多方向上[12]，但根據傅裡葉變換能量譜的自配准性，相同方向上的特徵線無論分佈如何，最終其能量總會疊加在一起[13]。而楔狀採樣作為對影像頻域能量分佈方向性的描述方法，其角向分佈特徵曲線的峰值將也必然會出現在與其能量譜的譜線相垂直的方向上。
統計分析兩期影像頻域能量的角向分佈後，即可得到兩期試驗影像的角向分佈特徵曲線[14]，從而表徵研究區域主要地物的排列方向，進一步可確定影像分塊的旋轉參數α。
若以採樣間隔為ri2≤u2+v2＜ri+12，在以頻譜圖中心為圓心的一個環面上對θ求和，其關係為
(4)
則稱Eθ(ri)為影像頻譜能量的徑向分佈特徵(環特徵)，如圖 2(b)所示。
應該注意的是，在一幅遙感影像中往往包含有各種週期性成分、非週期成分以及雜訊等資訊[15]。在徑向採樣中，除了直流分量外，如果有某一波峰特別突出，往往說明原圖像中包含以該波峰頻率為主頻率的週期成分存在，而且峰值的大小表徵了週期性成分的強弱程度[16-18]。
1.3 影像分塊參數的確定
大多數變化檢測方法需要對地表景觀進行完整的描述，這就要求在對影像進行分塊的同時，儘量保證影像中主要地物的完整性，不會因影像分塊導致地物破碎，從而避免影響地物提取精度以及變化檢測精度[19]。傳統影像分塊方式及其效果見圖 3。
圖 3 傳統影像分塊方式及分塊效果Fig. 3 The tradition way of image division and effectiveness
圖選項
在進行影像分塊的時候，為避免分塊線穿過主要地物，應根據影像的頻譜能量角向分佈曲線的峰值θ，對影像進行順時針旋轉，使影像中主要地物儘量按照經緯方向進行排列，以期在合適的步長下完成對地物的完整分塊。
為了使圖像既有足夠的紋理基元組合，又盡可能形成單一紋理，應針對不同影像中主要地物的紋理基元的週期性，設計影像分塊步長。論文選取累加距離匹配函數(summed-up distance matching function, SDMF)對圖像紋理基元的週期性進行描述[20]，通過分析影像中紋理基元的週期性，完成分塊步長的選取。
文獻[20]所提出的累加距離匹配函數，是在犧牲少量計算精度的代價下，提高傳統單視窗灰度共生矩陣(cray level co-occurrence matrix, GLCM)對紋理基元週期性的執行效率，其旨在計算影像不同方向間不同步長的亮度差值的累積量，從而判定紋理週期。行方向和列方向上的SDMF分別可以表達為
(5)
(6)
式中，M、N分別表示圖像的尺寸大小；r和c分別表示行數和列數；δ表示圖元間的空間距離。通常，可以將δ的範圍取作0~N-1，以獲取步長為1~N下的SDMF值，從而使SDMF的紋理週期計算具有全域性，進一步確定最佳分塊步長。
為了規範化探討SDMF在具備規則性紋理的影像中的特性，通過10個大小為43×344像元的矩形構建人工規則紋理圖像(見圖 4(a))，其中，設黑色亮度為0，白色亮度為1。對規則紋理圖像的SDMF做進一步的分析後發現，SDMF會產生大小、尺寸不一的波峰(見圖 4(b))，借助SDMF的一階、二階導數，計算各次波峰的峰值和穀值，通過設定閾值，判斷各次峰值和穀值的顯著性，獲取各顯著峰的尺寸，進而得到紋理圖像的週期。為了使顯著性判定閾值具有良好的自我調整性，經過規則紋理圖像實驗驗證，論文採用文獻[21]中的觀點，使用公式(7) 計算顯著性判定閾值
圖 4 人工規則紋理圖像及其SDMF分析圖Fig. 4 The artificial image and statistical analysis of SDMF in the grid row
圖選項
(7)
式中，SDMF的第一顯著峰值是全域最優的平均穩定週期[22]，第一顯著峰的尺寸即可表徵紋理基元的尺寸。但對於規則性紋理圖像，紋理週期是指基元特徵重複出現的尺寸大小[23]。因此，SDMF的所有顯著峰均可以表示紋理圖像的週期，但對於步長的選取而言，需要綜合考慮各種紋理基元的多樣性，又要照顧到地物的相對獨立性[24]。進而，論文規定步長的選取公式為
(8)
式中，S為影像分塊步長；μ為單位步長，即SDMF所反映的紋理週期；m為步長倍數。
應該格外注意的是，經過反復試驗，在研究區域中變化資訊對整體資訊的貢獻率低於52%的情況下，由兩期SDMF所確定的步長大小基本可以保證同一，此時將其最或是值作為單位步長即可滿足分析需求；而當影像中地物發生劇烈變化，導致變化資訊占比超過52%時，將難以保證兩期影像的SDMF所反映出的地物週期性一致。此時，應將兩期影像中較大的地物週期作為單位步長，並進行適當擴大步長倍數m，這樣雖然了犧牲了少量運算效率，但總體上可以保證變化資訊不會因影像分塊而丟失，並可以較好地反映地物週期性。
1.4 有效變化資訊所在影像塊的確定
紋理感知的精神心理學認為：人類識別紋理特徵的能力主要來自對紋理基元重複模式的判斷，其中3個重要的特徵是結構元素的規則度、重複性和方向性[25]。當研究區域地物發生變化時，必然會導致其紋理資訊和邊緣特徵的變化，從而帶來紋理規則度的改變，借此可以衡量圖像變化的程度，進而縮小變化區域的搜索範圍。
首先，利用紋理圖像頻譜的徑向特徵曲線E(r)統計其能量分佈情況，對徑向特徵曲線進行歸一化處理；其次，由其下包絡曲線B(r)對徑向特徵曲線進行擬合，計算E(r)和B(r)的能量比值e；再次，將B(r)中各處的點值乘以e，得到同能量下包絡曲線T(r)；最後，求取E(r)對T(r)的方差σEi，即確定圖像紋理規則度[26]。
通過兩期影像的紋理規則度，計算變化性紋理規則度ΔσE為
(9)
該演算法可以確保下包絡線和徑向特徵曲線所映射的幅度譜具有相同的能量，而且具有單調遞減的特性，從而有效減少多期影像間因光譜能量、輻照度等條件的不同對變化檢測所產生的干擾。
根據以上指標的定義過程，可以得到結論：圖像中結構性紋理所占的比例越大，圖像的紋理規則度就越高；圖像中包含的變化資訊越大，變化性紋理規則度則越高。這也就是說，變化性紋理規則度ΔσE可以反映影像塊產生變化的劇烈程度，通過設定對變化資訊分割閾值η，即可實現對變化資訊所在影像塊的準確判斷。
2 試驗和結果分析2.1 試驗
為了驗證本文所提出方法的有效性，選擇2012年5月17日和2015年8月11日獲取的雲南某地區QuickBird衛星遙感影像作為研究物件。QuickBird影像包含了4個多光譜波段和一個全色波段，其中全色波段波譜範圍為0.45~0.90 μm，全色波段空間解析度0.61 m。截取全色波段中2800×2800圖元大小的圖像作為研究資料。
首先，進行圖像預處理。論文採用文獻[27]中所提出的適合於去除高斯雜訊和椒鹽雜訊的偏微分方程模型進行圖像雜訊抑制處理。處理後發現，遙感資料亮度回應不明顯，為使圖像更好地呈現地物的細節結構特徵，將試驗資料進行長條圖拉伸，儘量使拉伸後的圖像長條圖分佈於0~255全部灰度級，以增加對比度，有利於地物類型的識別。
然後，計算影像分塊的兩個參數：旋轉參數α和分塊步長s。通過對兩期影像頻譜能量的角向分佈特徵進行統計分析，得其角向分佈特徵曲線，見圖 5和圖 6。
圖 5 2012年試驗影像的角向分佈特徵曲線Fig. 5 The image distribution of wedges feature curve in 2012
圖選項
圖 6 2015年試驗影像的角向分佈特徵曲線Fig. 6 The image distribution of wedges feature curve in 2015
圖選項
從影像角向分佈特徵曲線分析，雖然兩期影像的頻域能量分佈形式略有不同，但可以清晰地表徵出2012年和2015年試驗影像塊的頻域能量明顯集中於47°和137°，其餘的方向並未出現明顯的峰值，表現為一條較為平穩的曲線。這表明原影像中包含大量位於該能量峰值方向上的邊緣或直線資訊，而主要地物的優勢邊緣是按照47°的方向進行排列，而其餘地物特徵未表現出明顯方向性。
為了使主要地物按照經緯方向排列，儘量避免由於分塊不當而帶來的景觀破碎，取試驗影像角向特徵分佈曲線峰值方向的餘角α作為旋轉參數，將影像順時針旋轉47°，得到適合分塊的影像分佈方式，見圖 7。
圖 7 旋轉處理後的兩期試驗影像Fig. 7 Two experiment images of rotation in 2012 and 2015
圖選項
對2012年試驗影像的SDMF計算值進行統計，計算其SDMF的一階導數f′(x)、二階導數f″(x)，分析各次顯著峰的波峰與波谷的位置，從而準確反映出影像的紋理週期單位，該SDMF統計圖及其分析見圖 8。
圖 8 2012年試驗影像的SDMF統計分析結果Fig. 8 The SDMF statistical results of experiment image in 2012
圖選項
通過對一階、二階導數的判斷，得到2012年試驗影像SDMF計算結果中所有峰值和穀值。經計算，確定本實驗的閾值為4.606×107，統計各次顯著波峰值、穀值及其所對應步長，見表 1。
表 1 2012年SDMF的各次顯著波峰值、穀值統計表Tab. 1 The peak and valley value of SDMF and the associated significant wave statistical list in 2012
表選項
為了進一步考察在變化資訊所占比例相對較小的兩期影像中SDMF的分佈情況，統計2015年影像SDMF計算結果中各次顯著波峰值、穀值及其所對應步長，見表 2。
表 2 2015年SDMF的各次顯著波峰值、穀值統計表Tab. 2 The peak and valley value of SDMF and the associated significant wave statistical list in 2015
表選項
結合表 1和表 2可以看出，兩期影像中各次特徵點所對應的步長值僅出現微小擾動，但並沒有固定的偏離方向和大小，呈現出偶然誤差的特性。通過計算整條顯著波內各週期步長大小的最或是值，發現兩期遙感影像中SDMF所反映的單位步長分別為：219.655、219.725。這進一步說明在變化區域占比相對較小的情況下，兩期影像的SDMF所反映出的週期性規律基本具備同一性。
為了保證分塊後的子影像中包含儘量多的地物種類，結合試驗驗證，以兩倍顯著峰尺寸作為步長，即以440像元為步長，將影像分成81幅子圖像，其結果可基本保證影像塊中地物的完整性。
最後，在完成影像分塊後，分別統計各子圖像塊的頻譜能量徑向分佈，計算其同能量下包絡線，獲取各影像塊的紋理規則度，最終得到對應影像塊的變化性紋理規則度，並按其行列分佈統計為表 3。
表 3 試驗圖像的變化性紋理規則度Tab. 3 The value of TRCI about exemplificative images
表選項
通過對各影像塊的變化性紋理規則度進行統計分析後發現，除影像資訊空白區的變化性紋理規則度為0外，其餘54幅子影像塊均存在大小不一的變化資訊；進一步的，通過人工判斷表中不同紋理規則度在提取變化資訊能力上的強弱，確定變化資訊分割閾值η=0.03，通過對變化資訊的描述，有效變化資訊所在的影像塊已被準確標識(見表中加粗部分)，以子圖像塊(6, 2) 為例，驗證變化性紋理規則度對變化資訊描述的效果，見圖 9和圖 10。
圖 9 2012年試驗影像塊的歸一化徑向分佈曲線及同能量下包絡曲線Fig. 9 The normalized radial distribution curve and co-energy lower enveloping curve in image block in 2012
圖選項
圖 10 2015年試驗影像塊的歸一化徑向分佈曲線及同能量下包絡曲線Fig. 10 The normalized radial distribution curve and co-energy lower enveloping curve in image block in 2015
圖選項
結合兩期驗證子影像塊發現，在2012年試驗影像塊的歸一化徑向分佈曲線中，徑向分佈曲線和下包絡曲線的頻譜能量差之和為1.696，紋理規則度σE1=0.014；2015年試驗影像塊的徑向分佈曲線中存在明顯波峰，其徑向分佈曲線和下包絡曲線的頻譜能量差之和為2.432，紋理規則度達到σE2=0.062。
兩期試驗影像塊的變化性紋理規則度ΔσE=0.048>0.03。反映出兩期影像塊的結構性紋理發生較劇烈的變化，通過對試驗影像塊的目視解譯，可以發現該區域確實產生大量的地物變化，驗證了此方法的可靠性和有效性。
2.2 評價與分析
目前變化檢測的主流方法大致可以分為面向像元級、物件導向級、面向目標級3類，試驗將採用3類檢測手段中較為常用的方法，通過對比使用預處理方法前後的檢測精度及效率的變化情況，完成對變化預處理方法的評價。
對於像元級變化檢測方法，採用影像比值運算法，是指將兩期遙感影像中對應的像元值做比值運算，獲取差異影像，通過選取適當的閾值，實現對變化區域的提取。
對於物件級變化檢測方法，採用MeanShift-CVA法，是指通過MeanShift圖像分割演算法進行地理物件的提取，然後通過變化向量分析法(change vector analysis, CVA)計算兩個不同時期的影像差異，最後通過閾值判斷，完成變化資訊的提取。
對於目標級變化檢測方法，以建築物為例，通過傅裡葉變換實現影像從空域向頻域的變換，統計其頻譜能量分佈，獲取建築物的紋理區域和邊緣特徵，構建相適應的濾波器，完成建築物的二值化提取，最後通過對二值化資訊進行圖元級的變化檢測即實現變化資訊的提取。
通過以上3種代表方法在同一硬體環境、同一運算平臺下，對試驗影像進行變化檢測，統計其執行變化檢測預處理前後各演算法的檢測精度和效率(計入變化檢測預處理時間)，見表 4。
表 4 3種變化資訊檢測方法執行變化檢測預處理前後檢測精度和效率變化表Tab. 4 The comparisons of accuracy and efficiency before and after the preparing preprocess by three different change detection methods
表選項
針對統計結果，從宏觀層面看，3個級別的變化檢測方法在進行預處理後，檢測精度和運行效率均有所提高，這充分地反映了本文所提出的變化檢測預處理方法對提高檢測精度及其執行效率的有效性。
從微觀層面看，檢測預處理對3種檢測級別起到的作用有所差異。其中，圖元級比值法精度較未進行預處理前僅提高3.39%，執行時間減少5.79 s，運算效率提高10.80%，是3種方法中收效最差的，其根本原因在於圖元級變化檢測方法對高解析度遙感影像的不適用性，導致檢測預處理對其效果較差。物件級MeanShift-CVA法檢測精度提高9.54%，節省執行時間45.38 s，運算效率提高53.57%；在此基礎上，作為對檢測結果紋理資訊同質性的描述，其緊質度也有一定的提高。頻域下目標級變化檢測法精度提高12.50%，節省執行時間69.69 s，運算效率提高54.53%，在3種方法中，是收效最大的，這是因為在目標級變化檢測中，需要對兩期影像分別進行目標地物的提取，其往往會受到背景地物在光譜和頻域上的影響，通過對影像分塊，很大程度的抑制了目標地物的提取誤差，從而提高變化檢測精度。
進一步的，從演算法性質上來講，由於全自動方法在很多實際應用中會產生精度不足、演算法不夠健全和穩定等問題，尚無法適用於多種複雜情況[28]；故本文提出了一種人工干預的半自動化處理方法，借助複雜情況下的人工判斷，對檢測過程加以修正，最終達到提高檢測精度的要求，而恰當的干預次數和干預位置就成為半自動化演算法實現高精度、高效率的關鍵所在。論文在選取步長倍數和變化資訊分割閾值兩處操作上進行人工干預。經過試驗驗證，僅使人工干預次數在原演算法的基礎上平均增加34.12%，這意味著，所進行的人工干預不但可以對複雜情況下的參數進行彈性設置，而且不會對演算法原有的自動化程度造成過多的影響，從而側面驗證了本實驗的科學性。
應該注意的是，目標級變化檢測的精度很大程度上是由目標地物的提取精度決定，變化檢測預處理對提高目標級變化檢測精度的有效性，進一步從側面反映出，論文中的影像分塊方法並沒有引起大規模的地物破碎，而造成變化資訊丟失、檢測精度下降，再次證明了影像分塊方法的科學性和合理性。
3 結語
本文以頻域能量守恆定律和傅立葉轉換幅度譜的自配准性質為依據，通過合理的影像分塊，結合對應子影像塊的變化性紋理規則度，確定有效變化資訊的範圍；研究表明，高解析度遙感影像變化檢測的頻域分析預處理方法可以説明原檢測演算法實現更高層次的檢測目標，對於減少背景資訊的干擾、縮小變化區域搜索範圍具有一定的參考價值。但應該注意到的是，檢測精度及效率所能夠提高的程度主要由原檢測演算法的本質及其性質決定，如何弱化其對預處理方法的較大影響，使之更具普適性，是一個有待於進一步研究和解決的問題。
【引文格式】虢英傑，朱蘭豔。 高解析度遙感影像變化檢測的頻域分析預處理方法[J]. 測繪學報，2017，46(6)：743-752. DOI: 10.11947/j.AGCS.2017.20160570
更多精彩內容：
故其所對應的譜線必然會通過頻域平面的中心點，並集中反映在垂直于原特徵線方向的譜線上，這一現象稱之為傅氏變換幅度譜的自配准性質[10]。為了進一步探討影像中地物特徵與能量譜線之間的關係，需要將傅裡葉頻譜能量分析和統計方法結合起來考慮。本文採用文獻[11]中所提出的環狀採樣和楔狀採樣的方法，假設(r,θ)為(u,v)平面上的極座標，以採樣間隔為θi≤arctan(v/u)＜θi+1，在以頻譜圖中心為原點出發的扇面上對r進行求和，其關係為
(3)
則稱Er(θi)為影像頻譜能量的角向分佈特徵(楔特徵)，如圖 2(a)所示。
圖 2 能量譜的環狀採樣和楔狀採樣Fig. 2 Radius sampling and angle sampling of the frequency spectrum
圖選項
在高解析度遙感影像中，地物特徵線會隨機分佈在很多方向上[12]，但根據傅裡葉變換能量譜的自配准性，相同方向上的特徵線無論分佈如何，最終其能量總會疊加在一起[13]。而楔狀採樣作為對影像頻域能量分佈方向性的描述方法，其角向分佈特徵曲線的峰值將也必然會出現在與其能量譜的譜線相垂直的方向上。
統計分析兩期影像頻域能量的角向分佈後，即可得到兩期試驗影像的角向分佈特徵曲線[14]，從而表徵研究區域主要地物的排列方向，進一步可確定影像分塊的旋轉參數α。
若以採樣間隔為ri2≤u2+v2＜ri+12，在以頻譜圖中心為圓心的一個環面上對θ求和，其關係為
(4)
則稱Eθ(ri)為影像頻譜能量的徑向分佈特徵(環特徵)，如圖 2(b)所示。
應該注意的是，在一幅遙感影像中往往包含有各種週期性成分、非週期成分以及雜訊等資訊[15]。在徑向採樣中，除了直流分量外，如果有某一波峰特別突出，往往說明原圖像中包含以該波峰頻率為主頻率的週期成分存在，而且峰值的大小表徵了週期性成分的強弱程度[16-18]。
1.3 影像分塊參數的確定
大多數變化檢測方法需要對地表景觀進行完整的描述，這就要求在對影像進行分塊的同時，儘量保證影像中主要地物的完整性，不會因影像分塊導致地物破碎，從而避免影響地物提取精度以及變化檢測精度[19]。傳統影像分塊方式及其效果見圖 3。
圖 3 傳統影像分塊方式及分塊效果Fig. 3 The tradition way of image division and effectiveness
圖選項
在進行影像分塊的時候，為避免分塊線穿過主要地物，應根據影像的頻譜能量角向分佈曲線的峰值θ，對影像進行順時針旋轉，使影像中主要地物儘量按照經緯方向進行排列，以期在合適的步長下完成對地物的完整分塊。
為了使圖像既有足夠的紋理基元組合，又盡可能形成單一紋理，應針對不同影像中主要地物的紋理基元的週期性，設計影像分塊步長。論文選取累加距離匹配函數(summed-up distance matching function, SDMF)對圖像紋理基元的週期性進行描述[20]，通過分析影像中紋理基元的週期性，完成分塊步長的選取。
文獻[20]所提出的累加距離匹配函數，是在犧牲少量計算精度的代價下，提高傳統單視窗灰度共生矩陣(cray level co-occurrence matrix, GLCM)對紋理基元週期性的執行效率，其旨在計算影像不同方向間不同步長的亮度差值的累積量，從而判定紋理週期。行方向和列方向上的SDMF分別可以表達為
(5)
(6)
式中，M、N分別表示圖像的尺寸大小；r和c分別表示行數和列數；δ表示圖元間的空間距離。通常，可以將δ的範圍取作0~N-1，以獲取步長為1~N下的SDMF值，從而使SDMF的紋理週期計算具有全域性，進一步確定最佳分塊步長。
為了規範化探討SDMF在具備規則性紋理的影像中的特性，通過10個大小為43×344像元的矩形構建人工規則紋理圖像(見圖 4(a))，其中，設黑色亮度為0，白色亮度為1。對規則紋理圖像的SDMF做進一步的分析後發現，SDMF會產生大小、尺寸不一的波峰(見圖 4(b))，借助SDMF的一階、二階導數，計算各次波峰的峰值和穀值，通過設定閾值，判斷各次峰值和穀值的顯著性，獲取各顯著峰的尺寸，進而得到紋理圖像的週期。為了使顯著性判定閾值具有良好的自我調整性，經過規則紋理圖像實驗驗證，論文採用文獻[21]中的觀點，使用公式(7) 計算顯著性判定閾值
圖 4 人工規則紋理圖像及其SDMF分析圖Fig. 4 The artificial image and statistical analysis of SDMF in the grid row
圖選項
(7)
式中，SDMF的第一顯著峰值是全域最優的平均穩定週期[22]，第一顯著峰的尺寸即可表徵紋理基元的尺寸。但對於規則性紋理圖像，紋理週期是指基元特徵重複出現的尺寸大小[23]。因此，SDMF的所有顯著峰均可以表示紋理圖像的週期，但對於步長的選取而言，需要綜合考慮各種紋理基元的多樣性，又要照顧到地物的相對獨立性[24]。進而，論文規定步長的選取公式為
(8)
式中，S為影像分塊步長；μ為單位步長，即SDMF所反映的紋理週期；m為步長倍數。
應該格外注意的是，經過反復試驗，在研究區域中變化資訊對整體資訊的貢獻率低於52%的情況下，由兩期SDMF所確定的步長大小基本可以保證同一，此時將其最或是值作為單位步長即可滿足分析需求；而當影像中地物發生劇烈變化，導致變化資訊占比超過52%時，將難以保證兩期影像的SDMF所反映出的地物週期性一致。此時，應將兩期影像中較大的地物週期作為單位步長，並進行適當擴大步長倍數m，這樣雖然了犧牲了少量運算效率，但總體上可以保證變化資訊不會因影像分塊而丟失，並可以較好地反映地物週期性。
1.4 有效變化資訊所在影像塊的確定
紋理感知的精神心理學認為：人類識別紋理特徵的能力主要來自對紋理基元重複模式的判斷，其中3個重要的特徵是結構元素的規則度、重複性和方向性[25]。當研究區域地物發生變化時，必然會導致其紋理資訊和邊緣特徵的變化，從而帶來紋理規則度的改變，借此可以衡量圖像變化的程度，進而縮小變化區域的搜索範圍。
首先，利用紋理圖像頻譜的徑向特徵曲線E(r)統計其能量分佈情況，對徑向特徵曲線進行歸一化處理；其次，由其下包絡曲線B(r)對徑向特徵曲線進行擬合，計算E(r)和B(r)的能量比值e；再次，將B(r)中各處的點值乘以e，得到同能量下包絡曲線T(r)；最後，求取E(r)對T(r)的方差σEi，即確定圖像紋理規則度[26]。
通過兩期影像的紋理規則度，計算變化性紋理規則度ΔσE為
(9)
該演算法可以確保下包絡線和徑向特徵曲線所映射的幅度譜具有相同的能量，而且具有單調遞減的特性，從而有效減少多期影像間因光譜能量、輻照度等條件的不同對變化檢測所產生的干擾。
根據以上指標的定義過程，可以得到結論：圖像中結構性紋理所占的比例越大，圖像的紋理規則度就越高；圖像中包含的變化資訊越大，變化性紋理規則度則越高。這也就是說，變化性紋理規則度ΔσE可以反映影像塊產生變化的劇烈程度，通過設定對變化資訊分割閾值η，即可實現對變化資訊所在影像塊的準確判斷。
2 試驗和結果分析2.1 試驗
為了驗證本文所提出方法的有效性，選擇2012年5月17日和2015年8月11日獲取的雲南某地區QuickBird衛星遙感影像作為研究物件。QuickBird影像包含了4個多光譜波段和一個全色波段，其中全色波段波譜範圍為0.45~0.90 μm，全色波段空間解析度0.61 m。截取全色波段中2800×2800圖元大小的圖像作為研究資料。
首先，進行圖像預處理。論文採用文獻[27]中所提出的適合於去除高斯雜訊和椒鹽雜訊的偏微分方程模型進行圖像雜訊抑制處理。處理後發現，遙感資料亮度回應不明顯，為使圖像更好地呈現地物的細節結構特徵，將試驗資料進行長條圖拉伸，儘量使拉伸後的圖像長條圖分佈於0~255全部灰度級，以增加對比度，有利於地物類型的識別。
然後，計算影像分塊的兩個參數：旋轉參數α和分塊步長s。通過對兩期影像頻譜能量的角向分佈特徵進行統計分析，得其角向分佈特徵曲線，見圖 5和圖 6。
圖 5 2012年試驗影像的角向分佈特徵曲線Fig. 5 The image distribution of wedges feature curve in 2012
圖選項
圖 6 2015年試驗影像的角向分佈特徵曲線Fig. 6 The image distribution of wedges feature curve in 2015
圖選項
從影像角向分佈特徵曲線分析，雖然兩期影像的頻域能量分佈形式略有不同，但可以清晰地表徵出2012年和2015年試驗影像塊的頻域能量明顯集中於47°和137°，其餘的方向並未出現明顯的峰值，表現為一條較為平穩的曲線。這表明原影像中包含大量位於該能量峰值方向上的邊緣或直線資訊，而主要地物的優勢邊緣是按照47°的方向進行排列，而其餘地物特徵未表現出明顯方向性。
為了使主要地物按照經緯方向排列，儘量避免由於分塊不當而帶來的景觀破碎，取試驗影像角向特徵分佈曲線峰值方向的餘角α作為旋轉參數，將影像順時針旋轉47°，得到適合分塊的影像分佈方式，見圖 7。
圖 7 旋轉處理後的兩期試驗影像Fig. 7 Two experiment images of rotation in 2012 and 2015
圖選項
對2012年試驗影像的SDMF計算值進行統計，計算其SDMF的一階導數f′(x)、二階導數f″(x)，分析各次顯著峰的波峰與波谷的位置，從而準確反映出影像的紋理週期單位，該SDMF統計圖及其分析見圖 8。
圖 8 2012年試驗影像的SDMF統計分析結果Fig. 8 The SDMF statistical results of experiment image in 2012
圖選項
通過對一階、二階導數的判斷，得到2012年試驗影像SDMF計算結果中所有峰值和穀值。經計算，確定本實驗的閾值為4.606×107，統計各次顯著波峰值、穀值及其所對應步長，見表 1。
表 1 2012年SDMF的各次顯著波峰值、穀值統計表Tab. 1 The peak and valley value of SDMF and the associated significant wave statistical list in 2012
表選項
為了進一步考察在變化資訊所占比例相對較小的兩期影像中SDMF的分佈情況，統計2015年影像SDMF計算結果中各次顯著波峰值、穀值及其所對應步長，見表 2。
表 2 2015年SDMF的各次顯著波峰值、穀值統計表Tab. 2 The peak and valley value of SDMF and the associated significant wave statistical list in 2015
表選項
結合表 1和表 2可以看出，兩期影像中各次特徵點所對應的步長值僅出現微小擾動，但並沒有固定的偏離方向和大小，呈現出偶然誤差的特性。通過計算整條顯著波內各週期步長大小的最或是值，發現兩期遙感影像中SDMF所反映的單位步長分別為：219.655、219.725。這進一步說明在變化區域占比相對較小的情況下，兩期影像的SDMF所反映出的週期性規律基本具備同一性。
為了保證分塊後的子影像中包含儘量多的地物種類，結合試驗驗證，以兩倍顯著峰尺寸作為步長，即以440像元為步長，將影像分成81幅子圖像，其結果可基本保證影像塊中地物的完整性。
最後，在完成影像分塊後，分別統計各子圖像塊的頻譜能量徑向分佈，計算其同能量下包絡線，獲取各影像塊的紋理規則度，最終得到對應影像塊的變化性紋理規則度，並按其行列分佈統計為表 3。
表 3 試驗圖像的變化性紋理規則度Tab. 3 The value of TRCI about exemplificative images
表選項
通過對各影像塊的變化性紋理規則度進行統計分析後發現，除影像資訊空白區的變化性紋理規則度為0外，其餘54幅子影像塊均存在大小不一的變化資訊；進一步的，通過人工判斷表中不同紋理規則度在提取變化資訊能力上的強弱，確定變化資訊分割閾值η=0.03，通過對變化資訊的描述，有效變化資訊所在的影像塊已被準確標識(見表中加粗部分)，以子圖像塊(6, 2) 為例，驗證變化性紋理規則度對變化資訊描述的效果，見圖 9和圖 10。
圖 9 2012年試驗影像塊的歸一化徑向分佈曲線及同能量下包絡曲線Fig. 9 The normalized radial distribution curve and co-energy lower enveloping curve in image block in 2012
圖選項
圖 10 2015年試驗影像塊的歸一化徑向分佈曲線及同能量下包絡曲線Fig. 10 The normalized radial distribution curve and co-energy lower enveloping curve in image block in 2015
圖選項
結合兩期驗證子影像塊發現，在2012年試驗影像塊的歸一化徑向分佈曲線中，徑向分佈曲線和下包絡曲線的頻譜能量差之和為1.696，紋理規則度σE1=0.014；2015年試驗影像塊的徑向分佈曲線中存在明顯波峰，其徑向分佈曲線和下包絡曲線的頻譜能量差之和為2.432，紋理規則度達到σE2=0.062。
兩期試驗影像塊的變化性紋理規則度ΔσE=0.048>0.03。反映出兩期影像塊的結構性紋理發生較劇烈的變化，通過對試驗影像塊的目視解譯，可以發現該區域確實產生大量的地物變化，驗證了此方法的可靠性和有效性。
2.2 評價與分析
目前變化檢測的主流方法大致可以分為面向像元級、物件導向級、面向目標級3類，試驗將採用3類檢測手段中較為常用的方法，通過對比使用預處理方法前後的檢測精度及效率的變化情況，完成對變化預處理方法的評價。
對於像元級變化檢測方法，採用影像比值運算法，是指將兩期遙感影像中對應的像元值做比值運算，獲取差異影像，通過選取適當的閾值，實現對變化區域的提取。
對於物件級變化檢測方法，採用MeanShift-CVA法，是指通過MeanShift圖像分割演算法進行地理物件的提取，然後通過變化向量分析法(change vector analysis, CVA)計算兩個不同時期的影像差異，最後通過閾值判斷，完成變化資訊的提取。
對於目標級變化檢測方法，以建築物為例，通過傅裡葉變換實現影像從空域向頻域的變換，統計其頻譜能量分佈，獲取建築物的紋理區域和邊緣特徵，構建相適應的濾波器，完成建築物的二值化提取，最後通過對二值化資訊進行圖元級的變化檢測即實現變化資訊的提取。
通過以上3種代表方法在同一硬體環境、同一運算平臺下，對試驗影像進行變化檢測，統計其執行變化檢測預處理前後各演算法的檢測精度和效率(計入變化檢測預處理時間)，見表 4。
表 4 3種變化資訊檢測方法執行變化檢測預處理前後檢測精度和效率變化表Tab. 4 The comparisons of accuracy and efficiency before and after the preparing preprocess by three different change detection methods
表選項
針對統計結果，從宏觀層面看，3個級別的變化檢測方法在進行預處理後，檢測精度和運行效率均有所提高，這充分地反映了本文所提出的變化檢測預處理方法對提高檢測精度及其執行效率的有效性。
從微觀層面看，檢測預處理對3種檢測級別起到的作用有所差異。其中，圖元級比值法精度較未進行預處理前僅提高3.39%，執行時間減少5.79 s，運算效率提高10.80%，是3種方法中收效最差的，其根本原因在於圖元級變化檢測方法對高解析度遙感影像的不適用性，導致檢測預處理對其效果較差。物件級MeanShift-CVA法檢測精度提高9.54%，節省執行時間45.38 s，運算效率提高53.57%；在此基礎上，作為對檢測結果紋理資訊同質性的描述，其緊質度也有一定的提高。頻域下目標級變化檢測法精度提高12.50%，節省執行時間69.69 s，運算效率提高54.53%，在3種方法中，是收效最大的，這是因為在目標級變化檢測中，需要對兩期影像分別進行目標地物的提取，其往往會受到背景地物在光譜和頻域上的影響，通過對影像分塊，很大程度的抑制了目標地物的提取誤差，從而提高變化檢測精度。
進一步的，從演算法性質上來講，由於全自動方法在很多實際應用中會產生精度不足、演算法不夠健全和穩定等問題，尚無法適用於多種複雜情況[28]；故本文提出了一種人工干預的半自動化處理方法，借助複雜情況下的人工判斷，對檢測過程加以修正，最終達到提高檢測精度的要求，而恰當的干預次數和干預位置就成為半自動化演算法實現高精度、高效率的關鍵所在。論文在選取步長倍數和變化資訊分割閾值兩處操作上進行人工干預。經過試驗驗證，僅使人工干預次數在原演算法的基礎上平均增加34.12%，這意味著，所進行的人工干預不但可以對複雜情況下的參數進行彈性設置，而且不會對演算法原有的自動化程度造成過多的影響，從而側面驗證了本實驗的科學性。
應該注意的是，目標級變化檢測的精度很大程度上是由目標地物的提取精度決定，變化檢測預處理對提高目標級變化檢測精度的有效性，進一步從側面反映出，論文中的影像分塊方法並沒有引起大規模的地物破碎，而造成變化資訊丟失、檢測精度下降，再次證明了影像分塊方法的科學性和合理性。
3 結語
本文以頻域能量守恆定律和傅立葉轉換幅度譜的自配准性質為依據，通過合理的影像分塊，結合對應子影像塊的變化性紋理規則度，確定有效變化資訊的範圍；研究表明，高解析度遙感影像變化檢測的頻域分析預處理方法可以説明原檢測演算法實現更高層次的檢測目標，對於減少背景資訊的干擾、縮小變化區域搜索範圍具有一定的參考價值。但應該注意到的是，檢測精度及效率所能夠提高的程度主要由原檢測演算法的本質及其性質決定，如何弱化其對預處理方法的較大影響，使之更具普適性，是一個有待於進一步研究和解決的問題。
【引文格式】虢英傑，朱蘭豔。 高解析度遙感影像變化檢測的頻域分析預處理方法[J]. 測繪學報，2017，46(6)：743-752. DOI: 10.11947/j.AGCS.2017.20160570
更多精彩內容：