吳洪憲,鄭少蘭*
(1.廣東省國土資源測繪院,廣東 廣州 510500)
差異圖的構造和分析是影像變化檢測領域的主要方法。通常采用預處理、構造差異圖和閾值分割3個步驟進行無監督的影像變化檢測[1]。SAR影像變化檢測差異圖的構造方法主要是比值形式的,因其可有效抑制SAR影像內的乘性噪聲[2]。經典的比值形式差異圖構造方法包括比值法、對數比值法和區域均值比法[3]。Weken D V D[4]等將模糊數學中的貼近度概念引入影像變化檢測領域,將像素鄰域窗口作為模糊集合,采用比值形式的相似度公式來構建兩期影像的差異圖,進而評價兩個模糊集合的接近程度。該方法基于鄰域信息進行分析,對SAR影像中的乘性噪聲有較好的抑制作用,同時很好地保留了變化區域的邊緣細節。
過渡區域是圖像分割領域一個重要概念。ZHANG Y J[5]等率先將過渡區域應用到影像處理方面;LI Z Y[6]等分別利用均值濾波器和全通濾波器對差異圖進行處理,再通過計算二者間的灰度差值來提取過渡區域,雖然影像分割效果有所提升,但影像中的斑點噪聲也容易被誤分到過渡區域中。鑒于此,本文考慮將全通濾波器換成雙邊濾波器。雙邊濾波器不僅考慮了空間上像素間的接近程度,還增加了像素光度、色彩方面的權重。對中心像素的空間域和灰度域進行加權計算,既可有效保持邊緣細節,又可平滑噪聲,使圖面更純凈[7]。抽取過渡區域后,通過計算其像素均值確定影像的分割閾值。本文利用上述方法分別處理經典算子和模糊貼近度差異圖,并采用相關參數指標對變化檢測結果進行定量評價。結果表明,基于模糊貼近度差異圖的影像變化檢測效果更好。
本文采用的數據資料為變化檢測數據集,每種數據集包含同一區域的兩期變化影像和參考變化影像,其中參考變化影像用于變化檢測結果的定量指標評價。本文采用的變化檢測數據集包括加拿大Ottawa地區、西班牙Andasol地區和深圳市寶安區的實驗數據。其中,加拿大Ottawa地區實驗數據的兩期變化影像分別為1997年Ottawa地區雨季來臨前影像和雨季中影像,用于檢測洪水淹沒變化;西班牙Andasol地區實驗數據是該地區1987年和2013年兩期局部變化影像,用于檢測山區農村變化;深圳市寶安區實驗數據是該地區2018年和2020年兩期局部變化影像,用于檢測城市擴展變化。3個變化檢測數據集均為網絡下載,其中Ottawa地區數據為Radarsat SAR影像,Andasol地區數據為SeaSat影像,寶安區數據為Sentinal-1 IW_GRD級影像。
2.1.1 比值法差異圖
差異圖的構造通常都是基于預處理后的SAR影像。比值法差異圖的核心公式為;

式中,I1、I2為預處理后同一區域的兩期SAR影像。
2.1.2 對數比值法差異圖
比值法構造的差異圖對實際發生變化的區域檢測效果較好,但也會檢測出較多的偽變化區域。因此,Bovolo[1]等采用了對數比值法差異圖,構造公式為;

對數比值法可減少錯誤檢測數,但可能會造成變化區域的漏檢。
2.1.3 區域均值比法差異圖
上述兩種差異圖構造方法都是基于兩期影像對應位置上單一像素點的灰度信息進行操作的,無法較好地規避影像中的噪聲像點。因此,利用區域窗口像素信息構造差異圖的區域均值比法被提出。其構造公式為;

2.1.4 模糊貼近度差異圖
貼近度是模糊數學中提出的一個概念,常用于衡量兩個模糊集合的相似性程度。對于兩個時相SAR圖像S1和S2中對應位置(i,j|1≤i≤I,1≤j≤J)上的像素p1和p2,分別取其窗口內的鄰域N1和N2作為其對應的模糊集合。本文選取3像素×3像素的窗口,則Nk={Sk(i-1,j+Δj),Sk(i,j+Δj),Sk(i+1,j+Δj)},其中Δj={-1,0,1},k={1,2}。由于SAR影像中的噪聲主要是乘性的,比值法能更好地對其進行抑制,因此模糊貼近度差異圖的構造公式為;

式中,K為窗口內像素的個數。D(i,j)越大,表示兩個模糊集合的差異越大,貼近度越小,變化的概率越大;反之亦然。
分別利用經典算子和模糊貼近度構造差異圖后,對于差異圖位置(i,j)上的像素,其灰度級差異為B(i,j)與Dave(i,j)間的絕對差值,即

式中,B(i,j)為差異圖經雙邊濾波處理后的影像;Dave(i,j)為位置(i,j)上5像素×5像素窗口的像素均值。
Δg的大小與原圖相同,其值表示差異圖中像素灰度層面的差異。Δg均值和方差的計算公式為;

設置閾值為th=μ+a×σ,其中a為大于零的參數(本文設定a=0.8),抽取過渡區域。統計過渡區域中差異圖像素的灰度平均值T,再利用T對差異圖進行閾值分割,得到變化檢測結果。
變化檢測結果的可靠性需通過將檢測結果與實際發生變化區域的參考圖進行比較得到[2]。定量評價指標包括漏檢數、虛警數、總錯誤檢測數、漏檢率、虛警率等,還可通過Kappa系數的大小評估整體檢測性能。Kappa系數常用于評價分類的準確性,如表1所示,Kappa系數不同的值域代表不同級別的一致性,一致性越高分類越準[8]。

表1 Kappa系數的一致性分類表
不同方法下Ottawa地區實驗數據的變化檢測結果如圖1所示,可以看出,由于只針對單一像素進行運算,比值法和對數比值法的變化檢測結果視覺效果較差,虛警像素較多,對影像噪聲的抑制效果較差;區域均值比法和本文方法的結果較接近,畫面相對純凈,較好地抑制了影像噪聲,便于確定洪水淹沒區域的分布范圍。
不同方法下Andasol地區實驗數據的變化檢測結果如圖2所示,可以看出,對比參考變化影像,由于山體陰影的影響,4種方法的變化檢測結果中均存在一些虛警像素,畫面不夠純凈,其中本文方法的結果相對較好,抑制了大部分影像噪聲,更加接近該區域農村變化的真實情況。

圖2 Andasol地區實驗數據以及基于不同差異圖的變化檢測結果圖
不同方法下深圳寶安區實驗數據的變化檢測結果如圖3所示,可以看出,對照參考變化影像,4種方法均存在一定的虛警像素,其中比值法的虛警像素相對較多,區域均值比法和本文方法的變化檢測結果相對較好,與參考變化影像保持了很高的一致性。

圖3 寶安區實驗數據以及基于不同差異圖的變化檢測結果圖
Ottawa地區實驗數據變化檢測結果的定量性能評價數據如表2所示,可以看出,4種方法中比值法和對數比值法的總錯誤檢測數較多,Kappa系數相對較低;本文方法的結果與參考變化影像之間的Kappa系數為0.913,比比值法和對數比值法分別高0.226和0.365,僅比區域均值比法少0.006。

表2 Ottawa地區實驗數據的定量評價指標表
Andasol地區實驗數據變化檢測結果的定量性能評價數據如表3所示,可以看出,4種方法中本文方法的虛警率最低,Kappa系數最高,為0.714,比比值法、對數比值法和區域均值比法的Kappa系數分別高0.082、0.114和 0.024。

表3 Andasol地區實驗數據的定量評價指標表
寶安區實驗數據變化檢測結果的定量性能評價數據如表4所示,可以看出,4種方法的變化檢測效果均較好,Kappa系數均在0.840以上,與參考變化影像的一致性較高;其中本文方法的Kappa系數為0.905,比比值法、對數比值法和區域均值比法的結果分別高0.065、0.033和 0.012。

表4 寶安區實驗數據的定量評價指標表
綜上所述,針對不同的實驗數據,本文方法的各項指標在4種方法中是占優的;從變化檢測結果圖來看,由于本文方法利用窗口鄰域像素信息計算相似性測度,避免了單個像素運算時不好區分噪聲的缺陷,對SAR影像中乘性斑點噪聲的抑制效果更好,因此檢測結果的虛警率較低,圖面更加純凈,變化檢測結果與參考變化影像的一致性非常高,說明該方法準確率更高、可靠性更強。
本文分別利用3種經典算子和模糊貼近度構造兩個時相SAR影像的差異圖,在一定程度上抑制了乘性噪聲;再利用雙邊濾波器對差異圖進行處理,既可保持細節又能去除噪聲;然后通過計算差異圖中像素的灰度差異提取差異圖中像素的過渡區域,進而確定分割閾值;最后根據分割閾值分割差異圖,并利用相關指標定量評價不同差異圖的變化檢測效果。結果表明,基于模糊貼近度構造的差異圖得到的變化檢測結果影像與參考影像保持高度或幾乎完全的一致性,變化檢測準確率高、可靠性強,整體變化檢測效果優于其他3種經典算子。
利用本文方法得到的變化檢測結果仍存在一些雜點,這是由于雜點的灰度值接近變化區域內的像素灰度值,被誤分到灰度過渡區域,導致錯誤檢測。今后可研究變化區域的連通性,盡可能地清除誤分的雜點,從而提高方法的精度和可靠性。