基于目標相干性表征差異的多波段SAR相干變化檢測方法

冀廣宇 董勇偉 卜運成 李焱磊 周良將 梁興東

(中國科學院電子學研究所微波成像技術重點實驗室 北京 100190)

(中國科學院大學 北京 100049)

1 引言

相干變化檢測(Coherent Change Detection,CCD)是一種可以檢測場景中發生的微小變化的技術，它利用相位信息形成的相干性觀測度量可以探測出亞波長量級的變化，并以相干性變化差異的形式將其表現出來[1]。這是利用相干信息相對于非相干信息用于變化檢測的優勢之處。通過CCD處理形成的相干變化差異圖具體表現為：場景中發生變化的區域呈現低相干特性(相干系數趨于0)，未發生變化的區域呈現高相干特性(相干系數趨于1)。然而，場景中并非所有呈現低相干特性的區域均為我們所關注的變化區域(稱為目標變化區域)，也可能是低相干干擾造成的虛警。這些低相干干擾區域包括回波散射強度弱形成的低信噪比區域，極易受風吹雨淋造成體散射失相干的植被區域等。形成這些低相干干擾的根本原因是：CCD方法采用相位信息進行的相干性度量過于敏感，許多不是我們關注的區域也因發生足以用CCD方法檢測出的相位變化而被SAR傳感器敏銳地探測出來。因此，如何從低相干區域中排除引起虛警的干擾，是改善CCD性能的關鍵技術，最直觀的解決方法為對低相干區域分類。傳統的CCD方法只能區分高、低相干區域，不具備進一步區分能力。因此需要增加新的觀測量實現進一步的分類。

目前，國內外的學者針對這一問題開展了許多研究工作。文獻[2,3]從時間維度擴展，分別建立以時間為參量的模型對場景中的變化區域分類，提取出僅在特定時間發生變化的低相干區域；文獻[4]考慮全極化情況，利用極化相干矩陣建立檢驗統計量，對不同極化散射類型的變化場景分類，提取出屬于表面散射的冰雪消融變化；文獻[5]通過估計場景中的噪聲功率建立統計量，將真實變化區域與低信噪比區域區分，去除場景中低信噪比引起的低相干虛警干擾；還有的文獻利用圖像中目標的散射特征，分別采用聚類，半監督甚至深度學習的方式進行場景分類，實現變化檢測[6–8]。以上各方法或是側重于觀測方式的選取，或是關注于場景本身的散射屬性，沒有將二者有機地結合起來考慮，因而在適用范圍上有一定的局限性。CCD方法如同使用一把具有精細刻度的尺子，能夠丈量微小的變化，傳統CCD方法之所以產生低相干虛警干擾，是因為這把尺子的尺度雖精細但單一，能夠同時觀測到目標變化與虛警干擾但無法對其進一步區分。將CCD方法在頻率維度上擴展，采用不同頻率電磁波對場景觀測，利用場景中目標對多波段雷達信號的散射差異形成多尺度觀測，使得對目標變化區域與虛警干擾區域形成不同的觀測結果，從而將其有效區分。目前進行多波段CCD研究的文獻較少。文獻[9,10]采用極化變化檢測的方法處理C, L兩個波段的農作物場景SAR圖像，文中所采用的方法將兩個波段的數據結合起來，可檢測出的變化區域為兩個單獨波段變化檢測結果的并集。文獻[11]使用X, C, S, L不同波段的SAR數據對北極地區分別進行CCD與非相干變化檢測(NCCD)處理，通過對比實驗結果得出結論：(1)CCD比NCCD能檢測出更微小的變化；(2)NCCD處理結果與所使用的波段關系不大；(3)CCD中，觀測電磁波波長越短(頻率越高)，變化區域去相干程度越高。該文獻只是對實驗結果做初步定性分析，并指出：需要進一步對檢測出的變化進行分類才具有應用意義。

在多波段觀測模式下，不同波段的電磁波對場景探測，獲得的目標散射特性存在較大差異，這與電磁波在不同頻率下的穿透性能，可探測目標最小截面積，觀測目標變化去相干程度等特性有關。不同類型的目標發生不同尺度的變化在各個波段探測下形成不同的相干變化檢測結果。我們依此對各種變化分類，并從中提取我們感興趣的目標變化區域。基于上述原理，本文充分分析不同波段電磁波的目標探測特性及其對不同尺度變化的去相干程度，據此提出一種多波段CCD方法。該方法首先獲取場景在不同波段下的相干變化差異圖像，然后根據目標在不同波段下的相干變化差異表征用改進的期望最大化(Expectation-Maximization, EM)算法進行分類，根據先驗知識確定目標變化區域所屬的類別，最后用Dempster-Shafer(DS)證據理論實現多波段相干變化差異圖融合，提取目標變化區域。由于自然界中絕大多數場景對不同波段電磁波照射下呈現明顯不同的散射差異[12]，故采用本文提出的多波段CCD方法可適用于軍事、農業、災害監測等更廣闊領域的場景中。通過實驗數據處理結果表明，本文方法可以有效減小除目標變化區域外低相干干擾區域造成的影響。

本文的結構安排如下：第2節介紹不同波段電磁波探測目標特性的差異，第3節為多波段CCD算法原理與處理流程，第4節為實驗數據驗證，第5節對本文內容作總結。

2 目標的多波段相干性表征

場景中目標受電磁波照射，產生的響應為場景中各個目標散射中心響應之和，即目標反映在SAR圖像上的復投影結果。我們參考文獻[13,14]建立的參數模型，獲取的SAR圖像I表示為：

其中，E為場景中目標總響應，Nnoi為噪聲項。Ei為場景中每個目標的散射中心響應，包括復散射系數項和相位項。其中，前者與電磁波頻率f，目標觀測角度信息(俯仰角與偏航角)及目標結構參數T相關，后者與電磁波頻率f及目標散射中心距離參數有關。為了避免空間去相干對CCD結果的干擾，在數據采集過程中，我們保持場景變化前后的兩次觀測俯仰角θ與偏航角φ固定，使兩次觀測之間的復相干系數γ不包含空間去相干。根據噪聲去相干與變化去相干理論可知，影響復相干系數γ的重要因素為SAR圖像信噪比與發生變化的程度[15]，這兩點分別對應于場景中每個目標的復散射系數項與相位項，因此對于每個目標，其相干系數統計值可表示為。

關于場景中的不同結構目標在不同頻率電磁波探測下對相干系數的影響，已有一些文獻從不同角度直接或間接地通過理論分析與實驗驗證得出結論[15–19]，可總結如下：

(1) 高波段雷達可通過CCD方法檢測小尺寸或表面粗糙程度較弱的目標造成的微小尺度變化，同時對大尺寸或表面粗糙程度強的目標造成的變化更具有良好檢測性能；

(2) 低波段雷達通過CCD方法可檢測大尺寸或表面粗糙程度較強的目標造成的較大尺度變化，對于小尺寸或表面粗糙程度弱目標的幅值響應弱，使得圖像信噪比低，導致相干性差，易引起虛警；

(3) 中間波段雷達對目標的檢測效果介于高、低波段之間，使得在目標尺寸結構或表面粗糙程度處于一定范圍內的情況下，對于一定尺度目標變化的檢測結果具有區分度；

(4) 低波段雷達獨有穿透特性，可以檢測深層或植被冠層以下目標發生的較大尺度變化。

根據上述結論可知，在場景中存在尺寸和表面粗糙度有較大差異的目標，且發生不同尺度變化的時候，得到的CCD結果會有許多值，可據此把場景中的目標分成許多類；每種目標在不同波段電磁波探測下得到不同大小的CCD結果，分別分屬于每個波段的不同分類。我們利用這種分類及不同波段結果中分類的差異性設計多波段CCD方法，具體在第3節中闡述。

3 多波段CCD方法

根據上一節分析的目標在多波段探測下的相干性表征特性，本文提出一種多波段CCD方法。該方法能夠充分利用不同波段電磁波探測目標及目標變化造成的散射差異進行相干變化檢測，通過對多波段相干變化差異圖分類，結合目標在不同波段探測下的散射特性構成的先驗知識提取目標變化區域，排除CCD結果中的虛警。多波段CCD方法的主要步驟如下：

步驟1 對場景多波段SAR圖像進行基于SIFT特征的圖像配準與輻射校正，使之滿足CCD處理需求[20]。

步驟2 分別計算配準后各個波段變化前后SAR圖像的相干變化統計量(一般為相干系數)，形成各波段的相干變化差異圖。

步驟3 對每個波段分別處理：以該波段的相干變化差異圖為數據，建立基于相干變化檢驗統計量概率密度函數(probability density function,pdf)的有限混合模型(Finite Mixture Modem,FMM)，用改進的EM算法自適應求解FMM中分模型數量n及各個分模型pdf參數，據此在該波段對檢測場景按照變化情況分類，形成該波段觀測下的分類結果，并確定場景在該波段觀測下對每種分類的響應度。

步驟4 在各個波段中，根據場景待檢測區域中感興趣地物的少量監督樣本在該波段相干變化差異圖的直方統計結果確定目標變化區域與非變化區域所屬分類，然后根據響應度結果采用DS證據理論進行多波段相干變化差異圖像融合，獲取最終的多波段CCD結果。

基于上述處理步驟形成的多波段CCD方法流程圖如圖1所示。

上述處理過程中，步驟1和步驟2已有較多文獻進行相關方法研究，可以分別參考文獻[21]與文獻[1]中的算法進行處理。下面針對能夠體現多波段CCD方法特點的步驟3和步驟4進行詳細討論。

3.1 變化檢測場景分類

根據第2節的理論，受雷達頻率，目標尺寸，表面粗糙程度，變化尺度等因素影響，使用不同波段探測場景中的目標獲取的CCD檢測結果中存在相干程度的差異。因此，需要根據每個波段獲取的相干變化差異圖對檢測場景的變化情況依相干度大小進行分類，為下一步根據分類結果實現多波段融合提供可靠的數據。

在每個波段數據中，我們需要借助步驟2采用相干變化檢驗統計量的pdf進行場景分類。這里我們不妨采用傳統的相干系數檢驗統計量，其pdf表示為[22]：

式(2)中，|γ|表示相干系數的真實值，N表示像元統計個數(統計窗口大小)，為超幾何分布，(|γ|，N)可作為式(2)所示pdf的一組參數。傳統CCD只是根據檢驗結果相干性的高低將場景分成變化/非變化兩類，本文描述的多波段CCD中，分類情況更為復雜：受目標結構影響，場景未發生變化時|γ|值并非固定在1附近，可能因信噪比差異而在 [0，1]區間內取值；受目標變化程度影響，場景發生變化時|γ|值并非固定趨近于0，也可能在[0，1]區間內取值。因此根據不同目標散射特性與變化程度，依據相干統計結果可將場景分為許多類。每一類受目標紋理信息影響，統計像元數據之間的相關性可導致N的實際值與名義值不符，且場景中不同分類的N值也可能不同。據此，我們將某一波段觀測下檢測場景的相干變化差異圖看作是基于式(2)pdf分布的FMM，將每個單一分布看作一個分類(分模型)，總體分布表示為：

其中，ai為第i分類的加權系數，滿足，與第i個分布的形狀參數(|γi|，Ni)一起構成第i個分布的分模型參數。n為FMM分類個數。本文以EM算法為基礎實現基于FMM的變化場景分類。

EM算法是一種極大似然估計算法，它先構建初始分類，通過迭代的方式用最大似然準則調整樣本的類別歸屬，其中E步計算場景中每個像元對當前模型參數中的每個分模型的響應度，M步根據E步計算的響應度估計每個分模型新的模型參數。這樣循環迭代，直至模型參數滿足收斂條件，實現混合分布模型中各個分模型的參數估計。將傳統EM算法用于本文應用中存在以下兩個問題：一是分模型個數n固定，而本文應用中，分模型的個數未知；二是在EM算法的M步對分模型參數估計時，無法像高斯混合分布模型那樣建立關于分布參數最大似然估計的顯式表達式，也無法獲取本文分布的對數累積量顯式表達式[23]。對此，本文方法做出如下改進。

針對第1個問題，本文采用自適應估計方法。首先，建立一個足夠多分模型個數的混合分布模型。然后在迭代過程中，如果某一個分模型的加權系數ai小于某個閾值，則可認為該分類無意義予以去除，并重新計算剩下分模型的加權系數；如果有兩個分模型參數|γi|之差小于某個范圍時，可認為它們屬于同一分類，將其合并，具體操作為：將兩個分模型加權系數相加，取加權系數較大的參數|γi|，Ni共同作為合并后的分模型參數。這樣實現混合模型個數的自適應估計。

針對第2個問題，本文在M步中采用中心矩估計量估計分模型參數。對于式(2)所示的分布，其各階原點矩表達式[22]為：

對于統計數據，我們用樣本加權平均方法計算每個分模型的樣本k階中心矩，計算公式為：

其中，ξ(i，j)為第j個樣本像元在第i個分模型處的響應度。我們通過式(6)所示的極值計算方法獲取分模型參數。

在實際處理中，我們建立關于參數|γ|和N的k階中心矩2維索引圖代入式(6)，從而避免了大量復雜耗時的計算。

3.2 多波段融合

我們將3.1節獲取的變化檢測中第f個波段的第i個分類結果記作，為便于下文分析，在各個波段中，不妨將每個分類結果按照由小到大順序排列，即令。接下來采用圖像融合方法獲取多波段融合后的相干變化差異圖。

多波段圖像融合方法主要有基于DS證據理論[24]和基于模糊集理論[25]的融合方法，它們都具有不確定信息的處理能力。針對本文對多波段差異圖像融合的具體應用，我們采用沒有對分類狀態響應度函數作近似的DS證據理論融合方法。

在進行多波段融合之前，我們先根據DS證據理論框架構建變化狀態空間，并需要根據一定的先驗知識確定各個波段中每個分類隸屬的變化狀態。根據3.1節描述，每個分類所屬變化狀態與該波段探測下目標的信噪比以及變化程度有關，在這些參數難以預估的情況下，選擇少量監督樣本可以輔助確定分類所屬的變化類型[3]。監督樣本的選擇有兩個準則：一是選擇的樣本區域在場景中為具有一定代表性的地物，二是所選擇區域具有明確的變化狀態(變化/不變化)，且這兩種狀態均有與之相符的樣本。根據監督樣本所屬的變化狀態及其相干變化差異圖的直方圖統計結果確定變化類集合與非變化類集合。對于可能存在的分類集合，其|γ|值介于變化與非變化之間，它們可能屬于低信噪比的非變化情況，也可能屬于一定程度變化去相干的微小變化情況，在沒有先驗知識的情況下，令表示不確定情況。剩下其它類集合表示為 ? 。這樣構成了關于變化狀態的完備空間。令場景中的每個像元對狀態空間中每種狀態的概率分配函數等于對該狀態類別的響應度，這樣可保證每個像元中對所有狀態的概率分配函數之和為1。根據證據合成方法，變化狀態的多波段融合結果表示為：

其中，mf(·)表示第f個波段的概率分配函數，⊕表示DS證據融合操作。這樣對場景中逐個像元求解，即可獲取多波段融合相干變化差異圖。

需要說明的是，所有波段電磁波對水面等幾乎完全鏡面反射區域的散射強度都比較弱，一般均低于噪聲水平，因此在各個波段的相干變化差異圖中均呈現低相干特性。對這類區域的變化檢測超出了本文多波段CCD方法的檢測范圍，因此我們假定該區域沒有發生變化，將其看作虛警干擾從SAR幅度圖像上做掩模處理將其濾除，即：在所有波段變化前后SAR圖像上，幅度均低于某一閾值(噪聲水平)的像元，所在區域認為未發生變化，在最終的多波段CCD結果中將該位置數據置為1。

4 實驗數據驗證

采用ESAR機載SAR數據對本文方法進行實驗驗證。該數據變化前后采集時間間隔為1個月(2007年2月－2007年3月)，地點為瑞典Eggby小鎮郊外的一片場景。對該場景分別用L波段與P波段HH極化雷達進行變化前后數據采集，所得SAR圖像及對照光學圖像如圖2所示。在變化前后的兩次數據采集時間間隔內，場景中發生土壤翻修，草地修剪，冰雪消融等類型的目標變化，主要發生在場景中草地，裸地的部分區域以及湖面區域，在圖2(e)的光學圖像中將部分變化區域予以簡要標注。經過圖像配準后對每個波段SAR圖像進行CCD處理，得到的相干變化差異圖如圖3所示。比較圖3中兩幅圖可發現，場景中左側大面積植被區域在L波段差異圖中呈現低相干特性，在P波段差異圖中為高相干特性，這是因為相比于L波段，P波段電磁波具有良好的穿透特性，沒有受到植被冠層受風吹雨淋發生的變化去相干影響；場景中局部裸地、草地區域(紅色圓圈部分)在L波段中呈現高相干特性，而在P波段中呈現低相干特性，這是因為該處場景表面粗糙度較弱，導致波長較長的P波段對該處的目標散射強度較弱，造成低信噪比，從而受到噪聲去相干的影響。這兩部分區域均不是我們感興趣的目標變化區域，在各自波段的檢測結果中形成虛警干擾。

基于圖3所示兩個波段的相干變化差異圖，我們用3.1節提出的FMM及EM算法對場景分類。結合SAR圖像與差異圖信息，我們初始制定分模型數目均為5, M步迭代中取式(6)參數k=2。經過迭代后，L波段與P波段結果中分模型數目分別為4和3。圖4分別給出了兩個波段差異圖的統計直方圖與建立的FMM以及分模型的最終迭代結果，可見FMM對統計直方圖的擬合效果良好。

在圖3的兩幅相干變化差異圖中，我們分別選取黃色和綠色方框區域作為變化區域與非變化區域的樣本監督數據，并將各自波段的直方圖與各自的擬合分布模型比較，結果如圖5所示。為方便比較，圖5中將樣本監督數據直方圖乘以相應的倍數，使其與擬合分布分模型的概率密度大小相當。通過比較我們得知，L波段數據中，變化類與非變化類分別對應第1分模型與第3、第4分模型；P波段數據中，變化類與非變化類分別對應第1分模型與第3分模型。根據3.2節處理方法可得，兩個波段的第2分模型分別對應各自的不確定類。然后計算像元在各個分類中的響應度(即計算概率分配函數)，代入式(7)計算，求得多波段融合相干變化差異圖結果 (mP⊕mL)(C)如圖6(c)所示(場景中無因散射強度弱以致掩模去除區域)。圖6(a)與圖6(b)為使用模糊集理論的融合方法[25]，前者使用加權算術平均的方式融合，后者使用加權幾何平均方式融合。結合圖3與圖6可見，相比于單一波段CCD結果，多波段融合結果中，既排除了L波段中植被受風吹擾動造成的低相干干擾，又免受P波段中低散射強度導致低信噪比造成的低相干干擾影響，使得結果中只剩下土壤翻動，草地修剪，湖面冰雪消融造成的目標變化區域，更加接近真實值。在模糊融合結果與證據融合結果比較中，比較圖6(a)–圖6(c)3種多波段融合結果，可見加權算術平均的模糊融合結果在檢測閾值大于0.5情況下存在較多虛警，其它兩種融合結果均具有較大的閾值選取動態范圍，具體檢測性能需要通過詳細指標分析。

為了進一步說明本文使用的多波段CCD方法的性能，采用受試者工作特征(Receiver Operat-ing Characteristic, ROC)曲線與最優Kappa系數兩個指標進行分析。ROC曲線越靠近左上角，最優Kappa系數值越大，說明方法性能越好。圖7(a)，圖7(b)分別展示了實驗采用的兩個單波段CCD與多波段CCD方法的ROC曲線和選取不同檢測門限的Kappa系數。從圖中可以看出，本文使用的DS證據融合多波段CCD方法的ROC曲線更接近左上方，最優Kappa系數值為0.81，大于兩個單波段CCD方法(L波段：0.63，P波段：0.54)和兩種模糊融合方法(均為0.79)。因此，以上分析驗證了本文提出的多波段CCD方法對去除單波段CCD方法中低相干干擾，降低虛警率方面的正確性與有效性。

5 結束語

本文根據場景中目標及其發生的變化在多波段電磁波觀測情況下表征出的相干性差異，提出一種多波段CCD方法。該方法先將待檢測場景根據各波段的相干變化差異圖用改進EM算法進行各自的自適應分類處理，然后依據少量監督樣本分別在各個波段中確定分類的變化狀態，最后使用DS證據理論方法處理，獲取多波段融合CCD結果，去除了每個波段各自產生的低相干干擾，達到降低虛警率的目的。由于用于多波段CCD處理的重軌SAR數據較少，本文中僅用一組雙波段數據驗證該方法的有效性與正確性。接下來的研究應以更多波段對不同場景探測的重軌SAR數據為基礎，進一步分析波段種類、波段數目、場景地物類型等因素對多波段CCD方法實現效果的影響。

致謝本文作者感謝歐空局提供ESAR機載L波段與P波段重軌SAR數據。