利用高分三號SAR影像進行雙側變化檢測

蔡宣宣，張永紅，崔斌，2，3

(1.中國測繪科學研究院，北京 100039；2.武漢大學 測繪學院，武漢430079；3.城市空間信息工程北京市重點實驗室，北京 100089)

0 引言

利用遙感影像進行地表覆蓋和變化檢測已經越來越受到廣泛的關注和研究，其應用方向也在不斷增加，如民用領域的城市擴張研究[1-3]、農業及森林的調查監測[4]，軍事領域中的打擊效果評估、艦船位置變化監測等[5]。與傳統光學遙感不同，SAR衛星作為主動式微波成像衛星，成像不受天氣和光照條件限制，也正因這一特性，使利用SAR衛星進行地表覆蓋變化檢測有著更加廣闊的利用前景和能力[6]。高分三號作為我國首顆C波段合成孔徑雷達衛星，于2016年8月發射完成，是我國高分對地觀測系統的重要組成部分，擁有聚束、條帶成像、掃描成像等5大類12種成像模式，成像模式數量為世界上同類衛星之最，最高分辨率可達1 m[7-8]，其具體衛星參數如表1所示。在海洋、減災、水利和測繪等領域都有著重要的研究意義和應用前景。利用高分三號進行SAR變化檢測，達到了利用國產衛星進行全天時，全天候進行變化檢測的目的，極大改變了我國民用星載SAR圖像依靠進口的狀態，為國內各行業用戶提供高質量、高精度對地觀測數據[9]。對國產衛星應用和SAR變化檢測研究都有著極其重要的意義。

表1 高分三號屬性信息

單極化SAR影像變化檢測主要分為3個步驟，其中第二和第三步為當前的主要研究方向[10]。分別為：①圖像預處理。包括影像的配準、定標、裁剪、濾波等，為下步生成差異圖提供2個時相影像。②差異圖生成。這步主要是通過對2個時相影像的差異運算，反應2時相影像間對應像素的差異度，主要有LR、領域比(neighborhood based ratio，NR)[11]、均值比(mean-ratio,MR)[12]、似然比(likelihood ratio，LLI)[10]等。③差異圖分割。其目的是通過分割算法將差異圖中的變化圖斑與非變化圖斑進行分離，被廣泛使用的算法包括大津法(OTSU)[13]、最小誤差法(KI)[14]等。

廣義高斯模型及KI閾值分割(GKIT)閾值分割算法在分割LR差異影像上有著較為良好的應用效果，但由于GKIT算法假設差異圖中只存在變化和非變化兩類，就導致其只能發現單側變化，極大削弱了其變化檢測的能力[15]。針對這一缺陷，有學者[16]提出了利用DGKIT的方法，其主要思想是假設差異圖中存在的后向散射增強，減弱和非變化類像素均服從廣義高斯分布，在滿足KI準則的條件下，確定出分割三類像素的閾值。但由于差異圖的隨機性，2種變化類像素之間不管是差異程度還是數量上都可能存在較大差距，導致在擬合和求取準則函數最小值的過程中，無法穩定獲取正確的分割閾值，可能會導致在直方圖中一側確定出2個閾值而另一側被忽略的情況，極大削弱了檢測結果的可靠性。

本文正是針對進行雙閾值變化檢測這一目標，提出了在傳統針對LR差異圖進行單側GKIT分割變化檢測的基礎上，改進為利用互為相反數的2幅LR差異影像，分別進行GKIT分割，得到雙側變化的初始分割閾值，這樣可以穩定獲取2個初始的雙側閾值，然后利用馬爾科夫隨機場(MRF)分割分別進行分割精化[17]，最后再通過差值閾值去除上步分割結果中偽變化區域，融合后得到最終變化檢測結果。這樣就避免了出現在直方圖一側的閾值被忽略的情況。文章首次以高分三號影像作為變化檢測實驗對象，對文中算法進行實驗并驗證了精度，最后簡要分析了研究區發生變化的原因，證明了高分三號在變化檢測尤其是艦船及近海區域進行變化檢測的能力。

1 算法介紹

文中采用算法流程圖如圖1所示。

圖1 高分三號數據進行變化檢測算法流程圖

1.1 圖像預處理

進行變化檢測之前，需要將2個時相原始單視復數據(single looking complex，SLC)影像經過配準、裁剪、強度提取、定標、多視處理等流程，之后為了降低SAR影像固有斑點噪聲的影響，對影像進行濾波，用于生成差異圖的強度影像[18]。

1.2 差異圖生成

假設上步中用于生成差異圖的2時相影像分別為I1和I2，由于SAR圖像自身的噪聲特性，應用LR不但能將乘性噪聲轉化為加性噪聲，而且能較好地反映出2時相影像間的差異程度，所以本文中采用LR的方法。但與之前所不同，本文采用生成雙向差異圖的方法生成2幅對應像素互為相反數的差異影像Ds和Db：

(1)

和

(2)

假設Db中，圖像中灰度值較大的像素代表后向散射可能發生增強的像素，且灰度值越大發生后向散射增強的程度就越大，而灰度較小的像素代表后向散射可能減弱的像素，且灰度值越小發生后向散射減弱的程度就越小；相對地，在Ds中，灰度較大的像素代表可能發生后向散射減弱的像素而灰度較小的像素代表可能后向散射增強的像素。

1.3 差異圖分割

本文采用以2次單側GKIT閾值作為初始分割閾值，再利用MRF分割進行精化的分割方法，并通過約束條件去除了對數比值較大但差值較小的偽變化區域，得到最終分割結果。

1)GKIT。假設在差異影像中，未發生變化的區域要遠大于發生變化了的區域。GKIT的主要思想是假設廣義高斯分布分布更適合描述LR差異圖的概率分布，且差異影像中的變化和非變化像素均服從廣義高斯分布。如公式(3)所示。

i∈(u，c)

(3)

式中：p(x|Ci)表示差異圖像x中某一類Ci的概率密度函數；ai和bi是常數，可以通過與均值參數mi、方差參數σi和形狀參數βi的關系求出，而均值、方差和形狀參數都可以根據擬合的概率密度函數中求出相關參數[15]。

(4)

實際上由于假設中符合廣義高斯分布的差異圖直方圖本身具有對稱性，且尖峰拖尾形態中尖峰部分對應差異圖中廣大未變化區域的灰度值，而兩側拖尾形態對應的分別是差異圖中后向散射系數增強與減弱區的灰度值范圍。所以其在檢測雙側變化(后向散射增強與減弱2種變化)時有很大缺陷。其擬合結果都將只是除擬合出最高尖峰外，準則函數較小的那一側的閾值，而另一側的變化閾值將被忽略。

所以具有近似對稱性的廣義高斯分布直方圖在以下兩種情況會有較好的效果，一是利用雙閾值選取方法同時將差異圖中的后向散射增強和減弱閾值確定，如DGKIT算法；二是差異影像中發生變化的變化類單一，均為單一后向散射增強或減弱變化。

本文在假設可擬合出閾值一側的變化像素的均值大于非變化像素均值(可優先擬合出直方圖尖峰右側的閾值)的前提下，嘗試通過對雙向差異圖分割來解決這一問題。對差異圖Db應用GKIT得到閾值Td(>0)確定后向散射增強區域的初始分割結果；由于此時無法獲得后向散射減弱區域的閾值，故只能大致獲取初始閾值之后進行精化處理。本文就以對Ds應用GKIT獲取的閾值Ts(<0)作為分割后向散射減弱區域的閾值，選取的閾值Td和Ts均為分割差異影像Db，最終得到分割后向散射增強區域的初始影像Gb和后向散射減弱區域的初始影像Gs。

2)MRF分割。由于GKIT是通過全局直方圖選取閾值的分割結果，而全局閾值的缺點在于沒有利用到局部領域信息，選取結果會受到局部噪聲影響，且在無法擬合的尖峰一側閾值選取結果不精確。故文章采用將上步GKIT分割結果作為初始分割，再利用MRF分割對2幅初始分割影像進行精化，一方面有效利用了領域信息，抑制了圖像的噪聲，減小孤立的變化像元與孤立非變化像元對最終分割結果的影響，另一方面也能精化無法擬合單側閾值的一側閾值選取不精確的問題。

MRF模型與Gibbs分布存在等效關系，其能較好描述圖像空間的局部領域關系，使其在圖像分割領域得到了廣泛應用[17，19]。其目的是通過優化算法對圖像進行遞歸求解系統能量函數最小值及涉及到的統計參數，獲取對應的標號場，得到最終精化后的分割結果圖Is和Ib。

1.4 分割后處理

文中實驗數據采用的是未經拉伸的原始數據，灰度跨度較大，所以LR方法不可避免會遇到像素比值結果虛高的情況，比如研究區中出現的2組像素，第一組出現在陸地上且實際發生變化的2個對應像素灰度值為20和46.4，其對數比值為0.98，差值為26.4；另一組出現在海面上實際沒有發生變化的2個對應像素分別為0.5和1.09，其對數比值同樣為0.98，而差值只有0.59，這也從一個側面反應了應用經典LR方法生成差異圖的一個缺點[20]。為解決這一問題，在得到分割結果之后，將前后時相影像差值的絕對值影像S中閾值小于t，但分割結果是變化的像素分別從變化結果Fs和Fb上去除，只保留差值大于閾值t的部分，得到變化結果Fs和Fb，融合后得到最終的變化結果圖F。

2 數據準備

2.1 前后時相影像數據

試驗區原始數據基本信息如表2所示，本文中格式均為方位向在前，距離向在后，且未經過影像地理編碼。研究區原始影像如圖2所示。

表2 高分三號試驗區數據基本信息

圖2 寧波地區前后時相原始影像

2.2 實驗區域概況

如圖3所示，研究區為浙江寧波市區及周邊海域。寧波位于浙江省東部，長江三角洲南翼，為浙江省經濟中心，經濟發展十分迅速。研究區的主要地物類型有城區、海域、灘涂，根據2016年數據，附近舟山港穩居全球第一大港口，海上船舶往來密集，可以用來測試高分三號的船舶檢測能力。

圖3 研究區范圍

2.3 潛在可能發生真實地表變化

SAR影像變化可分為后向散射增強變化和減弱變化，如果根據地物后向散射強度將地物分為強、中、弱三類，則水面和裸地可認為是后向散射較弱和正常區域，船舶和人工建筑物則具有較強的后向散射強度。故反應在2時相SAR強度影像上有后向散射增強與減弱的變化之分。影像獲取時間分別為2016-11-14和2017-03-10，期間由于植被變化和船舶移動等影響，引起影像上后向散射強度變化。

3 實驗與分析

3.1 數據處理

1)預處理。預處理分為配準、裁剪、強度提取、多視處理和濾波等過程。多視處理的視數比為5∶4，裁剪后影像大小為4 071像素×3 753像素，分辨率為14 m左右。生成差異圖時對影像進行了3×3增強LEE濾波處理，以減少斑點噪聲對圖像影響。濾波后結果如圖4所示。

圖4 濾波后前后時相強度影像

2)差異圖生成。與傳統的SAR變化檢測中只生成一副差異圖不同，本文不再以后時相除以前時相影像后取對數來確定差異影像，而是生成2幅互為相反數的差異影像Db和Ds，如圖5所示。

圖5 對數比差異圖像

3)差異圖分割。對經過中值濾波的差異圖進行GKIT分割，假設Db中后向散射增強區域均值大于非變化區域均值，對Db進行分割得到Tb(=1.002 9)分割后向散射增強區域閾值，Db中大于Tb的像素是后向散射增強像素。對Ds進行分割得到Ts(=-1.097 9)分割后向散射減弱區域，差異圖Db中小于Ts的像素被認為是后向散射減弱像素。初始分割結果如圖6所示。

圖6 GKIT初始分割結果

如圖7所示，GKIT在直方圖中只擬合出了尖峰右側的閾值，而尖峰左側無法擬合出閾值。故在選取初始閾值后采用MRF分割進行精化處理。

圖7 差異圖Ds各類型概率密度分布

經過MRF精化后的結果如圖8所示。

圖8 MRF精化分割結果

未經歸一化原始影像生成LR差異圖會帶來比值較大而差值較小的偽變化區域，如圖8(a)中大片白色海域變化區域，考慮通過設定差值差異閾值(本文t=5)，將偽變化區域去除，并將上步中后向散射增強變化結果和減弱變化結果進行融合，得到最終變化結果圖。其中后向散射增強結果和減弱結果如圖9所示，整體變化結果如圖10。后向散射增強像素、非變化像素、減弱像素分別為255、128、0。

圖9 偽變化區域去除后結果

圖10 最終變化結果

3.2 典型變化與非變化圖斑分析

圖11是圖10中區域1的2時相信息。變化主要原因是不同獲取時間艦船位置及沿海區域發生變化。沿海地區發生了散射強度增強，主要原因可能是由于海產品養殖增加了反射地物和沿海含水量的變化；城鎮地區發生此變化的原因主要是城市擴張使裸地變成城市建筑區。由于沒有同時獲取的光學影像，文中沒有給出2時相對應的光學影像。圖12是圖10區域2的2時相SAR影像，由于地表植被及土壤含水量不同引起地表變化。圖13所示為鄰近時間光學影像，可以看出地表覆蓋情況發生變化情況。

圖11 區域1 2時相影像

圖12 區域2前后時相影像

圖13 區域2 2時相鄰近時間光學影像

為計算方便，精度評價時，將后向散射增強和減弱像素均標記為變化像素(灰度值255)。其他情況下，仍然分別顯示為后向散射增強和減弱區域。通過計算錯檢率(false alarms，FA)和漏檢率(missed alarms，MA)，總錯誤率(overall error，OE)及Kappa系數4項作為評價指標進行精度評價。分別使用本文方法和只進行了去除偽變化區域而沒有利用MRF分割精化的GKIT方法。由于缺少真實地表覆蓋變化圖，本文選取2個區域作為樣本代替真實地表變化圖。區域1，18 391個非變化像素，900個變化像素；區域2，2 284個非變化像素，912個變化像素。人工勾繪的變化圖斑如圖14所示。

圖14 人工勾繪的變化圖斑

為了驗證本文方法的實用性，與DGKIT算法的檢測結果進行了對比，DGKIT算法的檢測結果如圖15所示，其中2個分割閾值分別為-1.985 7和-1.093 7，可見DGKIT只在單側確定了2個閾值而忽略了另一側的變化情況，檢測結果也反應出了這一結果，圖中大部分沒有用發生變化的區域都被錯分為減弱區域。與真實地表變化相差較大，故沒有對其進行精度評價，只是給出了其檢測結果。

圖15 DGKIT變化結果

區域1 2種方法分割結果和精度評價結果如表3和圖16所示。可見經過MRF分割精化之后的結果，FA和MA分別減少了0.15%和0.16%，在本身局部已經具有較高精度的檢測結果的基礎上，考慮到鄰域關系后，精度又有了部分提升，Kappa系數達到了0.889 1。而DGKIT算法的檢測結果，如圖16(c)中所示，與真實結果相差較大。而由本文算法檢測結果圖16(a)和圖16(b)可見，利用適合的算法對高分三號影像進行處理，能在近海艦船檢測(圖16(a)中大部分黑白斑塊)和沿海地區及海岸線檢測方面產生較為良好的效果。

表3 區域1精度評價

圖16 區域1變化檢測結果

區域2分割結果和精度評價結果如表4和圖17所示。在經過MRF精化分割結果后，盡管MA有小幅下降，但OE仍然有了1.15%的提升。最終的Kappa系數也達到了0.851 1。但DGKIT檢測結果，檢測結果全部為錯誤的黑色減弱像素，與真實地表變化相差較大。

表4 區域2精度評價

圖17 區域2變化檢測結果

3.3 變化結果分析

通過算法檢測，提取了整個影像覆蓋地區的變化情況，變化像素約占影像對應研究區大小的1.71%。從圖中可以看出，變化主要分布在2個主要區域，一是集中在市內的鄞州區，二是北侖區和鎮海區沿海。鄞州區由于區域內存在較多農業用地，相隔4個月后，土地植被種植情況及土壤含水量都發生了不同程度變化，但由城市建設引起的包含人工地物變化的情況則發生較少。相比之下，北侖、鎮海區的沿海地區發生變化的主要原因是因為沿海灘涂，包括海帶紫菜等近海養殖業及海面船舶位置變化引起的后向散射變化是造成這一區域變化的主要原因。

4 結束語

本文針對在SAR變化檢測中應用DGKIT算法在增強和減弱像素變化存在較大差異的情況下可能會產生單側出現2個閾值的情況，提出了利用雙向差異圖確定雙側初始閾值，并利用MRF分割精化，再去除結果中存在的偽變化區域，最終生成雙側變化檢測結果的方法，確定了2時相影像中發生后向散射增強和減弱的變化區域。

通過對高分三號浙江寧波地區2016年11月14日和2017年3月10日2景的實驗檢測結果的分析，一方面證明了本文算法是一種有效的SAR變化檢測方法，另一方面，首次利用高分三號SAR影像作為變化檢測實驗對象也顯示出國產高分三號SAR影像在變化檢測尤其是沿海灘涂及近海海產養殖、海岸線變化和艦船檢測上有著較好的檢測效果。