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

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

0 引言

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

表1 高分三號屬性信息

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

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

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

1 算法介紹

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

圖1 高分三號數(shù)據(jù)進(jìn)行變化檢測算法流程圖

1.1 圖像預(yù)處理

進(jìn)行變化檢測之前，需要將2個(gè)時(shí)相原始單視復(fù)數(shù)據(jù)(single looking complex，SLC)影像經(jīng)過配準(zhǔn)、裁剪、強(qiáng)度提取、定標(biāo)、多視處理等流程，之后為了降低SAR影像固有斑點(diǎn)噪聲的影響，對影像進(jìn)行濾波，用于生成差異圖的強(qiáng)度影像[18]。

1.2 差異圖生成

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

(1)

和

(2)

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

1.3 差異圖分割

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

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

i∈(u，c)

(3)

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

(4)

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

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

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

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

MRF模型與Gibbs分布存在等效關(guān)系，其能較好描述圖像空間的局部領(lǐng)域關(guān)系，使其在圖像分割領(lǐng)域得到了廣泛應(yīng)用[17，19]。其目的是通過優(yōu)化算法對圖像進(jìn)行遞歸求解系統(tǒng)能量函數(shù)最小值及涉及到的統(tǒng)計(jì)參數(shù)，獲取對應(yīng)的標(biāo)號場，得到最終精化后的分割結(jié)果圖Is和Ib。

1.4 分割后處理

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

2 數(shù)據(jù)準(zhǔn)備

2.1 前后時(shí)相影像數(shù)據(jù)

試驗(yàn)區(qū)原始數(shù)據(jù)基本信息如表2所示，本文中格式均為方位向在前，距離向在后，且未經(jīng)過影像地理編碼。研究區(qū)原始影像如圖2所示。

表2 高分三號試驗(yàn)區(qū)數(shù)據(jù)基本信息

圖2 寧波地區(qū)前后時(shí)相原始影像

2.2 實(shí)驗(yàn)區(qū)域概況

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

圖3 研究區(qū)范圍

2.3 潛在可能發(fā)生真實(shí)地表變化

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

3 實(shí)驗(yàn)與分析

3.1 數(shù)據(jù)處理

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

圖4 濾波后前后時(shí)相強(qiáng)度影像

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

圖5 對數(shù)比差異圖像

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

圖6 GKIT初始分割結(jié)果

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

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

經(jīng)過MRF精化后的結(jié)果如圖8所示。

圖8 MRF精化分割結(jié)果

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

圖9 偽變化區(qū)域去除后結(jié)果

圖10 最終變化結(jié)果

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

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

圖11 區(qū)域1 2時(shí)相影像

圖12 區(qū)域2前后時(shí)相影像

圖13 區(qū)域2 2時(shí)相鄰近時(shí)間光學(xué)影像

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

圖14 人工勾繪的變化圖斑

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

圖15 DGKIT變化結(jié)果

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

表3 區(qū)域1精度評價(jià)

圖16 區(qū)域1變化檢測結(jié)果

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

表4 區(qū)域2精度評價(jià)

圖17 區(qū)域2變化檢測結(jié)果

3.3 變化結(jié)果分析

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

4 結(jié)束語

本文針對在SAR變化檢測中應(yīng)用DGKIT算法在增強(qiáng)和減弱像素變化存在較大差異的情況下可能會產(chǎn)生單側(cè)出現(xiàn)2個(gè)閾值的情況，提出了利用雙向差異圖確定雙側(cè)初始閾值，并利用MRF分割精化，再去除結(jié)果中存在的偽變化區(qū)域，最終生成雙側(cè)變化檢測結(jié)果的方法，確定了2時(shí)相影像中發(fā)生后向散射增強(qiáng)和減弱的變化區(qū)域。

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