本文內容來源於《測繪學報》2022年第12期（審圖號GS京（2022）1419號）

基於超像素分割與多信息融合的疊掩和陰影區域檢測法

劉鵬¹

1. 西安電子科技大學雷達信號處理國家重點實驗室，陝西 西安 710071;

2. 地理空間信息國家重點實驗室，陝西 西安 710054;

3. 西安測繪研究所，陝西 西安 710054

基金項目：國家自然科學基金(62031005)

摘要：InSAR信號處理中, 疊掩和陰影區域嚴重破壞了InSAR相位圖的連續性, 導致高程反演出現錯誤。本文提出一種基於多信息融合的超像素檢測算法, 以彌補傳統的基於幅度和相干係數閾值檢測、恆虛警檢測和局部頻率檢測等逐像素檢測方法的不足。由於缺少真實場景的數字高程模型, 導致無法驗證所提算法在實際應用中的準確性, 因此本文提出了基於先進積分方程近似模型的頻域回波模擬技術, 用於模擬給定場景的SAR回波數據以驗證本文算法。仿真數據和TH-2數據試驗表明, 本文算法能夠區分SAR圖像中98%以上的疊掩區和陰影區。最後, 通過將檢測出的疊掩和陰影區域在相位解纏過程中進行掩膜, 進一步提高了相位解纏速度和DEM產品精度。

關鍵詞：天繪二號 疊掩 陰影 簡單線性疊代聚類 快速回波仿真

引文格式：劉鵬, 李真芳, 樓良盛, 等. 基於超像素分割與多信息融合的疊掩和陰影區域檢測法[J]. 測繪學報，2022，51(12)：2517-2530. DOI: 10.11947/j.AGCS.2022.20220111

LIU Peng, LI Zhenfang, LOU Liangsheng, et al. Layover and shadow regions detection based on superpixel segmentation and multi-information fusion[J]. Acta Geodaetica et Cartographica Sinica, 2022, 51(12): 2517-2530. DOI: 10.11947/j.AGCS.2022.20220111

閱讀全文：http://xb.chinasmp.com/article/2022/1001-1595/20221214.htm

引 言

干涉合成孔徑雷達(interferometric synthetic aperture radar，InSAR)技術得益於其獲取目標高程信息的能力，被廣泛應用於地形測繪、海洋監測、災害監測與預警等領域^[1-4]。由於合成孔徑雷達屬於主動輻射測距成像模式的一種技術，使得疊掩和陰影成為SAR圖像中普遍存在的現象。在InSAR信號處理中，疊掩和陰影區域的相位奇異性導致相位無法正確濾波和相位解纏，錯誤的相位信息導致相鄰區域的相位展開錯誤，並造成高程反演出現錯誤。因此需要準確判斷並標記SAR圖像的陰影區和疊掩區，避免該區域影響整個相位解纏過程。

近年來，許多學者對InSAR中陰影和疊掩的檢測方法進行了相關研究，並取得了一定成果，如基於幅度的閾值分割方法^[5]、基於相干係數的閾值分割方法^[6]、基於邊緣特徵和模糊理論的方法^[7]、基於局部頻率估計的方法^[8]、基於特徵值分解的方法^[9]和基於恆虛警(CFAR)檢測的方法^[10]等。其中，局部頻率估計和特徵值分解更為有效，但其都有各自的局限性。一方面，上述方法都是逐像素進行判別，SAR圖像實際上反映的是場景對電磁波的雙站散射能力，目標的幾何和物理特性通過影響場景電磁雙站散射而被SAR系統捕獲並成像，SAR圖像上每個分辨單元是由完全隨機分布的許多散射點矢量疊加組成，在SAR圖像上形成相干斑噪聲，進而使得簡單閾值分割方法失效，逐像素的判別結果較為離散，即使做形態學處理，也很難反映SAR數據的紋理信息；另一方面，局部頻率估計需要對每個像素周圍指定窗口進行二維FFT操作，對SAR數據逐像素進行判別計算量較大，特徵值分解方法需要獲取同一區域的多幅SAR圖像，通過構造圖像陣列的方式判別疊掩和陰影，無法應用於大部分工程需求，恆虛警檢測(constant false alarm，CFAR)基於目標與背景對比強烈的準則，需要估計目標的分布函數，對於較大的疊掩和陰影區域檢測性能不佳。

針對逐像素檢測方法在準確度和速度上的不足，考慮到在SAR圖像中屬於同一地物同等坡度的相鄰像素具有相似信息，疊掩或陰影區域具有散射信息近似相等且距離上相近的特點，如果能夠將具有相鄰且相似信息的像素組合成一個基元塊，則可以大大降低後續圖像處理任務的複雜性。基於此，本文結合超像素分割方法提出一種基於多信息融合的超像素檢測算法，實現了快速且準確的疊掩和陰影區域檢測。近年來，常用的超像素分割方法是基於圖論和梯度下降聚類的方法。基於圖論的方法具體包括基於隨機遊走的熵率的方法^[11]、偽布爾優化的方法^[12]和基於圖的貪心聚類算法^[13]。基於梯度下降聚類的超像素分割包括簡單線性疊代聚類(SLIC)方法^[14]、疊代邊緣細化方法^[15]、多尺度分割方法^[16]和分水嶺算法^[17]。這些方法由於其處理速度快、邊界擬合度較好，被廣泛應用於光學圖像分割中。其中，SLIC方法由於在遵守圖像邊界的能力、速度、內存效率以及分割性能上的優越性而被廣泛使用^[18]。然而，由於SAR圖像固有的相干斑噪聲效應^[19]，以及相對於光學圖像而言SAR圖像具有更多像素數量，通常達到10⁶~10⁸數量級，使得傳統SLIC在超像素分割方面表現不佳。近年來多尺度分割和SLIC方法在SAR圖像分割領域得到了廣泛的應用。儘管過去已提出了許多的多尺度分割方法，但嵌入在eCognition軟體中的多尺度分割(multiscale segmentation，MS)方法仍然是目前最有效的方法^[20]。針對SAR圖像存在的相干斑噪聲問題以及SAR數據固有的數據特徵，文獻[21]基於概率密度函數的SLIC(PDF-based SLIC)方法，實現了L波段的SAR圖像超像素分割。文獻[22—23]通過改進SLIC超像素算法中的距離度量方法實現了PolSAR圖像中目標的分類。然而SLIC算法應用於SAR圖像分割中依舊有些問題沒有很好地解決。首先，大部分研究都是集中在改進建立有效的距離度量以區分地物類別，並未考慮到初始化步驟對於SAR圖像分割的重要性；其次，疊代過程中複雜的距離度量使得算法在大尺寸SAR圖像中效率低，不利於實際工程應用。因此，需要考慮SAR圖像中疊掩和陰影檢測的實際情況，設計一種快速且準確的SLIC處理流程，以及將SLIC嵌入到後續處理中實現疊掩和陰影的快速檢測。

本文基於改進後的SLIC算法，結合相干係數、閾值幅度分割方法和局部頻率估計方法，提出一種基於多信息融合的超像素檢測算法，識別SAR圖像中的疊掩和陰影區域，來保證InSAR後續處理的可靠性和有效性。由於缺少真實場景的數字高程模型，導致無法驗證本文算法在實際應用中的準確性，因此本文提出一種基於電磁散射模型的快速回波仿真法，將先進積分方程近似模型(advanced integral equation model，AIEM)^[24-25]和頻域回波仿真^[26]相結合，通過利用衛星工具包(satellite tool kit，STK)軟體設計仿真所需的軌道參數，實現星載全極化雙站SAR的原始回波數據仿真，回波能夠反映多站、多角度、全極化、複雜地物的回波特性以及反映SAR圖像的疊掩、陰影、相干斑噪聲等物理特徵，可用於檢驗本文所提出的基於多信息融合的超像素檢測算法有效性。結合天繪二號衛星^[27-28]的2019年9月獲取的河北赤城精度檢測場區域和2020年6月獲取的哥倫比亞區域數據，證明本文算法的有效性。

1 SAR圖像疊掩和陰影檢測理論研究

1.1 疊掩和陰影的幾何特性

工程上利用兩個在俯仰向觀測角度略有差別的接收天線，對同一目標進行各自成像處理，通過干涉的方法獲取目標高程信息，稱其為干涉合成孔徑雷達(interferometric synthetic aperture radar，InSAR)技術。如圖 1所示，S₁和S₂代表裝載雷達天線的兩部衛星，當雷達觀測具有陡峭地形的山區或城市時，具有相同斜距的目標點P₁和P₂在錄取回波時會疊加到同一個像素內，將這一現象稱為疊掩，將電磁波無法照射的區域成為陰影區域。由圖 1中的幾何關係可以得到局部入射角θ的具體表達式為

圖 1 幾何關係

Fig. 1 Geometric relationship

圖選項

式中，β為雷達下視角；α為目標點的坡度角。當θ=0°時，斜坡的坡面與雷達視線垂直，在該目標點兩側的點將具有相同的斜距，使得不同層次高程的目標被疊加到同一距離單元內形成疊掩；當θ< 0°時，此時α>β，底部的目標點到雷達的斜距大於頂部的目標點到雷達的斜距，使得SAR圖像上出現頂底倒置的現象，進而形成疊掩現象。陰影區域一部分是由於電磁波被其他目標遮擋而形成，另一部分是由於θ>π/2，即α<β－π/2時電磁波鏡面反射導致。

1.2 基於信號特徵的傳統逐像素檢測法

1.2.1 基於幅度閾值分割方法

SAR圖像中疊掩和陰影在強度上有明顯的區別，文獻[9]指出可以計算整個SAR圖像的平均值A，認為疊掩區域Al>2A，陰影區域As< 0.15A。實際上水域或局部入射角θ接近90°的區域，它們通常都具有回波能量低的特點，這些區域在SAR圖像上的強度一般較弱，導致所檢測出的陰影區域比圖 1幾何約束下的陰影區域更多，過檢區域對後續高程反演是有益的。

1.2.2 基於局部頻率估計方法

文獻[8]指出InSAR系統形成的圖像之間存在頻率偏移Δf，並給出了具體表達式

(2)

式中，f₀為載波的頻率(單位Hz)；R為雷達到目標點的斜距(單位m)；k為常數，對於雙星收發分置雷達系統取2；B為基線長度(單位m)；γ為基線傾角(單位rad)。由式(2)可以看出，當β－α< 0或β－α>π/2時Δf取值為負，此時對應疊掩和陰影區域。

本文採用文獻[8]中頻譜分析方法估計局部頻率，在指定窗口內做二維FFT求干涉相位圖頻譜，頻譜峰值出現在二維局部頻率對應的位置，通過判斷頻率符號確定可能的疊掩和陰影區域。

1.2.3 基於CFAR檢測方法

均值類CFAR算法是目前主流的CFAR檢測算法，其實現思路是以被檢測單元為中心的周圍單元作為參考窗，將窗口內的所有數據平均作為背景功率估值，將估計的背景功率與係數(門限因子)相乘得到決策閾值T，根據不同的背景功率和門限因子。文獻[10]採用SO-CFAR算法進行檢測，將參考窗分為對稱的兩個窗口，以較小的背景雜波數據來估計背景功率，其背景功率為

(3)

式中，N為檢測窗口長度；L_i為窗口內像素幅度值；虛警率Pf可表示為

(4)

式中，ε_so為門限因子。SO-CFAR算法適用於多個相鄰的目標檢測，目標周圍虛警率低，但是雜波邊緣區域虛警概率提升。

2 SAR圖像疊掩和陰影檢測算法實現

2.1 SLIC超像素分割算法的改進

SAR圖像反映的是場景反射電磁波的能力，因此SAR圖像是強度圖像，傳統的SLIC算法適用於RGB圖像的分割，鑑於SAR圖像像素數量通常達到10⁶~10⁸數量級，需要在保證分割精度的同時儘可能減少傳統SLIC算法運算量以達到工程上快速分割的目的，本節對SLIC算法進行改進，並給出了改進後的SLIC具體實現流程。預處理時需要將SAR圖像轉換成取值範圍為[0, 1]的灰度圖像，然後放大100倍作為圖像的亮度信息，得到亮度圖L。

2.1.1 初始化聚類中心

對於一個輸入為R×B=N的SAR圖像，其中R和B分別為圖像的行數和列數，N為SAR圖像像素點總數。初始化K個聚類中心，原始SLIC算法將輸入圖像均勻劃分為K個矩形網格，如圖 2所示通過將種子點錯位分布，使得生成的超像素具有六邊形網格的特徵。首先，假設初始相鄰聚類中心的間距為

；然後，計算圖像第一列聚類中心個數後重新更新S。最終，每個聚類中心可以表示為fCk=[lCk,xCk,yCk]^T，k∈[1,K]。同樣圖像中的每個像素點可以表示為fi=[l_i, x_i, y_i]^T，i∈[1,N]；l表示對應像素的亮度，取值[0, 100]；x、y表示該點在圖像中的位置。

圖 2 初始化種子點

Fig. 2 Initializing the seed point

圖選項

為避免將超像素放置在圖像邊緣或噪聲點上影響聚類性能，本文將聚類中心移動到與八連通域中的最低梯度位置相對應的種子位置。梯度度量方法為

(5)

式中，L(x,y)代表圖像中坐標為(x,y)的點像素亮度信息；·代表L2範數。

2.1.2 像素點的聚類

SLIC算法通過檢測每個像素點周圍2S×2S範圍內聚類中心之間的相似度，將相似度最高的聚類中心的標籤賦予該像素。SLIC算法定義了最大的空間距離Ns=S和最大的顏色差異Nc來對空間距離和顏色距離進行差異整合，由於Nc的值不容易直接定義，所以SLIC算法將Nc取值為常量m∈[1, 40]，最終給出了像素間的相似性公式為

(6)

d(i,Ck)描述了像素i和聚類中心Ck之間的相似度。可以將式(6)表示為向量的形式，特徵向量f的加權L2範數形式為

(7)

式中，對角矩陣W=diag(1,λ,λ)，λ=m²/S²；f=[l,x,y]^T。則d(i,Ck)可寫成

(8)

對於每一個需要和像素i計算相似度的候補集群C_k，||fi||W²是他們公共的一項，此項對於任意一個關於像素i的運算都相等，因此可以省略。同時去除根號後，則式(8)最終簡化為

(9)

進一步通過近似的方法將d(i,Ck)近似到ds(i,Ck)減少了計算量，將聚類問題轉化為argminds(i,Ck)。

2.1.3 疊代更新

初始化標籤圖像l={l(i)=－1i∈[1,N]}，距離圖像D={d(i)=∞i∈[1,N]}，循環k∈[1,K]計算Ck鄰域2S×2S範圍內ds(i,Ck)，若ds(i,Ck) <d(i)，則更新d(i)=ds，l(i)=k，並更新聚類中心，計算誤差項E=Dn－Dn－1，Dn表示第n次疊代。循環上述操作進行疊代，終止條件為E< threshold。試驗表明疊代10次左右，大部分圖像都可以得到滿意結果。

2.1.4 後處理增強連通性

經過上述疊代優化可能出現多連通情況、超像素尺寸過小以及出現孤立超像素等，這些情況可以通過增強連通性解決。本文具體處理步驟如下：

(1) 為確保不相鄰區域的標記是不同的，將具有不同區域但相同標籤的重新分配標籤。

(2) 對標籤圖像進行形態學開操作處理，去除孤立的小點和毛刺。

(3) 上述操作完成後，某些區域可能會被關聯，分成多個區域或被其他區域吸收，為確保不相鄰區域的標記是不同的，需要對區域重新編號，以便他們從1開始順序增加。

上述步驟可能導致超像素個數增多或減少，最終輸出標籤圖像l={l(i)i∈[1,N]}，其中l(i)表示SAR圖像上每個像素對應的超像素編號，對應的超像素特徵結構可表示為

(10)

式中，l′Ck表示超像素Ck所包含所有SAR圖像中原像素的平均強度信息；x′Ck、y′Ck表示超像素Ck的等效位置信息；cohCk表示超像素Ck的平均相干係數信息；K′表示最終分割的超像素個數。

2.2 基於多信息融合的超像素檢測算法

基於改進後的SLIC超像素分割算法，結合InSAR相干係數、傳統閾值分割和局部頻率估計方法，本文提出了一種疊掩陰影判別算法。具體流程如圖 3所示。

圖 3 算法流程

Fig. 3 Algorithm flowchart

圖選項

對於給定輸入為R×A=N的SAR圖像。首先，採用均值濾波的方式去除SAR圖像中的斑點噪聲；然後，對SAR圖像中的輪廓特徵進行精細化分割，獲得標籤圖像l={l(i)i∈[1,N]}和超像素特徵結構SPCk,k∈[1,K′]參數。初始化疊掩、陰影和疑似疊掩陰影的標籤圖像Il、Is和Ils。計算SAR圖像平均強度L_avr，遍歷k∈[1,K′]，對每個超像素做局部頻率估計得到疊掩和部分陰影區的區域集合Ils，通過強度和相干係數閾值分割法在區域Ils中識別疊掩區Il，則可得到部分陰影區Is=Ils－Il。估計出部分陰影區的平均強度avr(L>(Is))，將2×avr(L>(Is))和相干係數小於0.6作為最終判定陰影區的閾值進行陰影區域的檢測。陰影區域通常沒有回波能量，因此相干性較差，根據處理過程中的經驗，一般認為相干係數小於0.6的區域包含陰影區域。在相干係數差的區域中通過設置閾值提取陰影區域。閾值分割得到的陰影區域中同樣包含回波能量弱的區域，這是由於SAR圖像每個像素值反映的是電磁波與目標點作用後被接收機接收的回波能量值，鑑於散射係數受到入射角影響較大，對於入射角大於70°的面目標而言，回波能量相對較弱，使得在數據分布特性方面與陰影區域相似，造成這部分區域相干係數較差。對於高程反演而言，相干係數較差的區域也會影響相位生成精度，因此回波能量弱的區域也應該被有效檢測。

2.3 基於回波模擬技術的算法有效性檢驗

為驗證本文算法的正確性，需要明確疊掩和陰影的實際位置，由1.1節的分析可知，對於給定的觀測場景數字高程模型(DEM)，通過坡度和下視角的關係可判斷疊掩區域，通過射線追蹤的方法可判斷陰影區域。因此通過準確地模擬回波錄取過程以獲取真實場景的SAR數據，可以檢驗本文算法的有效性。SAR回波過程可被描述為發射信號與系統衝激響應的卷積^[29]，將其擴展到雙基，具體可表達為

(11)

式中，τ為快時間參數(單位s)；t為慢時間參數(單位s)；f為頻率參數(單位Hz)；接收信號SRx(τ,f,t)被描述為式(11)中STx(τ,f)與通道脈衝響應函數hc(τ,t)的卷積以及0均值高斯白噪聲n₀(τ,f,t)的形式。在快時間的相干窗口內假設發射機位置PosTx(t)和接收機位置PosRx(t)相對時間是不變的。此時目標點的通道脈衝響應可被寫為

(12)

式中，δ(τ－R(t,r)/c)為雙基距離延時的δ函數；σ(t,r)為在採樣時刻t對應的WGS-84坐標系位置r=(x，y，z)處目標點的反射率函數，反射率函數是由觀測幾何角、目標點材質等特性決定的；雙基距離函數R(t,r)為

(13)

在每個方位向慢時間採樣時刻，4 dB波束寬度內對應的照射區域內的脈衝響應hc(τ,t)可被描述為

(14)

式中，W(t,r)是由雷達發射功率、收發天線增益以及雷達工作過程中由熱損失等引起的系統損耗等因素共同決定的參數；St為慢時刻t時4 dB波束照射到的場景目標集合，去除其中的陰影區域。根據雙站SAR雷達方程^[29]可得到W(t,r)的具體表達為

(15)

式中，Pt為發射信號的功率；Gt(t,r)為慢時間t接收位置r處目標的發射天線方向圖增益；Gr(t,r)為接收天線方向圖增益。

SAR成像過程實質上是重構目標點反射率函數的過程，因此SAR數據仿真的準確性取決於反射率函數計算的準確性。文獻[24]將AIEM模型與實測數據進行對比，證明了AIEM模型在計算隨機粗糙面雙站散射係數的可靠性。AIEM模型考慮基爾霍夫項、交叉項和補償項，可以簡化寫成最終的簡單形式，散射係數計算公式為

(16)

式中，q、p在極化h和v中取值表示4種極化的狀態；k=2π/λ為波數，λ為波長(單位m)；δ為均方根高度(單位m)；kx=ksinθcosϕ,ky=ksinθsinϕ,kz=kcosθ，θ為局部入射角，ϕ目標點的入射方位角；ksx=ksinθscosϕs,ksy=ksinθssinϕs,ksz=kcosθs，θs為散射角，ϕs目標點的散射方位角；關鍵函數Iqpn的計算詳細見文獻[25]。W(n)(K)如式(17)所示用來描述指數相關函數隨機粗糙面的起伏情況

(17)

式中，K=sqrt((ksx－kx)+(ksy－ky)²)；ρ為相關長度(單位m)。

在AIEM模型的推導坐標系下ϕ=π，因此只要給出δ、ρ、θ、θs、ϕs這5個參數以及介電常數即可準確計算目標點的散射係數。通常對於目標場景採用三角面元建模，每個場景面元都有一個確定的觀測幾何角θ、θs、ϕs，若用A表示三角面元的面積，則每個面元接收電磁波的有效面積為Acosθ。由於SAR圖像存在相干斑噪聲，因此模擬的SAR圖像也必須具有斑點噪聲特性。文獻[19]指出斑點噪聲被建模為實部和虛部服從0均值高斯分布的乘性噪聲。因此最終的反射率函數表達為

(18)

式中，N(0, 1)表示服從均值為0方差為1的高斯分布特性的數。

根據上述建模分析，最終的雙站SAR回波方程可被描述為

(19)

3 試驗結果及分析

一方面，通過將疊掩和陰影區域檢測之後的結果作為掩膜圖屏蔽疊掩和陰影區域進行干涉相位解纏，可以大大縮短解纏時間並提高解纏精度。對檢測出的疊掩和陰影區域進行標註，以便後續通過升降軌融合方法或已有的DEM對疊掩和陰影區域進行修正，生成高精度的DEM產品。另一方面，TH-2致力於構建全球DEM庫，因此必須儘可能提高疊掩和陰影的檢測效率，而超像素分割過程作為疊掩和陰影檢測的重要步驟，需要提高其運算效率和分割性能。為了證明本文提出的改進SLIC算法的有效性，對TH-2獲取的SAR數據集進行了性能和運算效率分析，並將本文算法與目前已被開發能夠用於SRA圖像分割的方法(傳統SLIC^[18]、MR方法^[20]及PDF-based SLIC方法^[21])進行對比。其中，MR方法是集成在商業軟體eCognition v8.9中的方法，在該版本中MR方法並未使用多核並行計算技術進行加速。為對比4種算法的運算效率，其餘方法均採用C++程式語言實現且不進行並行加速。

為了評估不同分割方法的性能，本文採用3個指標來評估分割結果，分別為邊界契合度(boundary fit，BF)、平均變異係數(mean coefficient of variation，MCV)和運算效率。

BF參數用來衡量超像素分割邊界與真實邊界的一致性，定義為

(20)

式中，ES∩T用來描述通過算法分割的邊界和地物真實邊界的重合像素數；ET表示地物真實邊界的像素數。因此BF的值域為[0, 1]之間，BF的值越接近於1表示超像素塊與圖像邊緣的一致性更好。本文中BF參數選取局部分割區域進行統計求解。

通常不同分割算法分割出的超像素塊之間差異較大，因此每個超像素塊的平均值略有差異，採用方差難評估多種算法之間的性能，CV參數可以消除由於不同超像素塊的平均水平不同對離散程度比較的影響。MCV參數定義為

(21)

式中，N表示本次分割的超像素總數；超像素塊Ck表示當前處理的第k個超像素塊；M為當前的超像素塊中的像素個數；l(Ck,i)表示超像素Ck中第i個像素的強度；l′Ck表示超像素塊Ck的平均強度值。MCV值越小，則分割結果的波動程度越小分割算法的分割性能越好。

如圖 4所示，試驗數據選自TH-2號2020年6月獲取的哥倫比亞區域實測數據，大小為18 214×23 566像素，截取其中大小為3000×3000像素的SAR圖像作為輸入。在MR方法中，尺度參數設置為50，形狀參數設置為0.4，顏色參數設置為0.5。在SLIC方法中，初始超像素個數K取值均為10 000，顏色和空間重要性度量參數m取值15，其餘所需參數均參考對應文獻。由圖 4結果可以看出，SAR圖像中通過分割得到的超像素區域中的顏色是均勻的，即屬於同一類的像素在回波能量上具有相似的表現。從細節上可以看出不同算法在性能上有一定的差異，主要體現在邊界的定位以及分割效率上。對圖 4中矩形的區域進行評估，由圖 4中黃色框和表 1中BF和MCV參數結果可以看出MR方法和本文算法都能較好地分割出不同區域的邊界，而PDF-based SLIC方法和傳統SLIC方法分割性能相對較差。本次試驗中採用MR方法進行分割用時18.28 s，本文算法用時1.99 s，結果表明本文算法在更短的時間內分割出的超像素更加均勻和緊湊，這為後續疊掩和陰影區域檢測提供了重要的保障。

圖 4 TH-2數據超像素分割結果

Fig. 4 Superpixel segmentation results for the TH-2

圖選項

表 1 性能比較

Tab. 1 Performance comparison

s

表選項

在TH-2號2020年6月獲取的哥倫比亞區域數據集中截取[500, 1000, 2000, …, 10 000]²數據量大小的圖像作為輸入，在每次試驗中超像素個數設置為圖像行數的10倍，顏色和空間重要性度量參數m固定為15。具體結果如圖 5所示。

圖 5 運算效率

Fig. 5 Operating efficiency

圖選項

隨著輸入數據量的增大，MR方法需要根據相異性準則對空間上非相鄰的區域進行合併使得運算量增大，PDF-based SLIC需要估計像素周圍的概率密度函數進行聚類使得計算量增加，由於本文算法在計算上的有效近似以及後處理的優化使得計算效率明顯優於其他方法。

3.2 仿真驗證

本文利用STK軟體生成觀測河北赤城精度檢測場區域的軌道數據，通過構造500 m基線得到發射機和接收機對應的天線相位中心文件，其餘具體仿真參數見表 2。

表 2 仿真參數

Tab. 2 Simulation parameters table

表選項

觀測區域的DEM見圖 6(a)，本文採用西安電子科技大學InSAR課題組研發的RDSpace軟體對仿真回波數據進行成像，本文給出發射星的成像結果見圖 6(b)，對比圖 6(d)成像後地理編碼的結果和圖 6(a)可以證明本文回波仿真方法的有效性。

圖 6 回波仿真SAR成像結果

Fig. 6 SAR results of echo simulation

圖選項

對於3500×3588像素大小的SAR圖像，分別設置初始超像素個數K=^{[3, 10, 3]}×10 000，顏色和空間重要性度量參數m=^{[15, 15, 5]}。為方便呈現分割後的結果，本文選取圖 7(a)紅色框出來的區域進行分析，圖 7(a)為紅色框中的SAR圖像，3種分割參數均實現了對SAR圖像的分割，並很好地保證了SAR圖像中目標的輪廓信息。對比圖 7(b)、(c)可以看出當m取值相同時，超像素個數越多分割的結果越細，分割完成後的超像素個數接近初始的超像素個數，初始超像素越多適應目標邊界的能力越強。對比圖 7(b)、(d)可以看出當初始超像素取值相同時，m取值越小超像素形狀和尺寸越不規則，但是分割出的超像素更加均勻和緊湊，邊界粘附性更好。從圖 7(d)結果可以看出，生成的超像素個數比輸入的超像素多，這是因為參數m取值越小顏色權值越大，此時算法更傾向於像素強度值上的鄰近性，使得在初次分割時大量相同標籤的區域並未接壤，在後續處理中為確保不相鄰區域的標記是不同的，對具有不同區域但相同標籤的重新分配標籤，進而使得超像素個數增多。在疊掩和陰影檢測時，更關注不同類別區域地物能量的變化，因此通常m取值較小，為保證結果的可靠性，一般m取值範圍為1~40。

圖 7 SLIC超像素分割結果

Fig. 7 The results of SLIC superpixel segmentation

圖選項

超像素形狀和尺寸的規則性並不能提高疊掩和陰影檢測正確性，對於疊掩和陰影的檢測通常更關心目標能量上的變化。通常希望分割的結果能夠更緊密地黏附到目標邊界，因此設置初始超像素個數K=30 000，顏色和空間相對重要性度量指標m=5進行疊掩和陰影的檢測。並將本文算法與文獻[5—6, 10]以及根據幾何關係判斷的真實疊掩和陰影區域進行對比驗證，詳細結果見圖 8。圖 8(a)、(e)是根據DEM和軌道信息幾何關係準確判斷的疊掩和陰影位置，由於SAR圖像反映的是目標點的反射率函數被脈衝響應函數調製後的點散布函數形式，在該試驗中SAR圖像一個像素上的點是5 m×4 m範圍內的DEM點的疊加，因此首先需要將DEM反投影到SAR圖像上再進行插值，才能確保高程點與SAR像素對應。可以看出CFAR檢測方法出現很多漏檢區域和無效區域，疊掩的檢測正確率不足50%。本文算法幾乎檢測出全部的疊掩區域，未檢測的部分大多是離散的孤立點區域，這些區域通常只占1~2個像素在形態學處理時被去除，對後續InSAR處理影響可忽略。對於陰影區域而言，陰影區域多為噪聲，回波能量較弱，因此陰影區域幾乎都可以根據強度值檢測出來，由圖 8(h)中可以看出，檢測出的陰影區域比真實陰影區域較多，這是因為這些區域普遍集中在背坡面，具有較大的入射角，幾乎沒有回波信號，導致在圖像上呈現與陰影區域相似的特徵。

圖 8 不同方法檢測結果

Fig. 8 The results of different methods

圖選項

為驗證檢測結果的有效性，定義疊掩和陰影檢測正確率為檢測出的真實疊掩和陰影數占準確的疊掩陰影總數比例，定義疊掩檢錯率為檢測錯誤的總像素占檢測出的總像素的比例，為驗證本文算法對DEM產品誤差精度的提升效果，通過計算未掩膜處理時DEM生成的高程誤差均值以及經過疊掩和陰影區域掩膜處理後的高程誤差均值，並計算相位解纏速度提升指標和DEM產品誤差精度提升指標。具體指標結果見表 3。

表 3 疊掩陰影檢測正確率

Tab. 3 Detect the correct ratio of layover and shadow

(%)

表選項

由表 3可以看出，CFAR檢測方法對於複雜地形的SAR圖像表現效果較差，相干係數聯合幅度閾值檢測的方法檢測正確率相對於本文算法較低，具體是因為相干係數是由周圍信息估計得到的，幅度閾值是按照像素點去檢測的，所以實際的相干係數和像素點的幅度會有不一致性導致正確率降低，錯檢率增加，這一現象在邊界處更為明顯。通過將本文算法檢測出的疊掩和陰影區域作為掩膜圖，在相位解纏過程中進行掩膜使得相位解纏時間從16.64 s降低到13.13 s，相位解纏速度提升17.31%，生成的DEM產品高程誤差均值從2.32 m降低到1.91 m，DEM產品誤差精度提升17.58%，優於其他方法。圖 9給出了仿真數據InSAR處理後生成的DEM產品高程誤差值，誤差分析試驗區域為圖 9中紅色線段位置，圖 9(b)為未掩膜的DEM生成誤差，圖 9(c)給出了掩膜後的DEM生成誤差，DEM生成誤差是通過InSAR處理得到的DEM產品和原始回波仿真輸入的觀測場景DEM做差得到。圖 9中對比紅色方框中的結果可以明顯看出，由於解纏錯誤導致誤差較大的區域，通過被有效地掩膜使得高程反演誤差降低，證明本文算法的有效性。

圖 9 DEM產品高程誤差

Fig. 9 Elevation errors of DEM product

圖選項

3.3 天繪二號數據試驗

將本文所提算法應用在天繪二號衛星在2019年9月獲取的河北赤城精度檢測場SAR圖像中，對於6046×6971像素大小的SAR圖像，設置初始超像素個數K=30 000，顏色和空間重要性度量參數m=5。考慮到原始接收到的SAR數據中由於地物信息複雜存在強散射點數據，導致如果直接對SAR圖像量化後顯示，造成SAR圖像顯示較暗問題，若不進行處理會使得後續陰影檢測失效，因此需要去除SAR數據中後5%強度值數據再進行後續試驗。

如圖 10所示，圖 10(a)、(c)為本文算法檢測出的疊掩區域以及將疊掩區域邊界標記到SAR圖像上結果，圖 10(b)、(d)為本文算法檢測出的陰影區域以及將陰影區域邊界標記到SAR圖像上結果。為證明本文算法的有效性，選取圖 10中綠色框中的區域的SAR圖像以及對應的相干係數圖，分別採用本文算法、CFAR檢測方法和相干係數聯合幅度閾值檢測方法進行試驗。試驗結果如圖 11所示，從中可以看出相干係數聯合幅度閾值檢測的方法檢測出的疊掩區域較真實疊掩區域更大，無法適應疊掩的邊界，檢測出的陰影區域較小，且疊掩和陰影漏檢嚴重。相干係數聯合幅度閾值檢測方法受到SAR圖像固有的斑點噪聲影響檢測結果離散，檢測性能較差。實際SAR圖像相比於仿真的SAR圖像而言斑點噪聲更嚴重，使得相干係數聯合幅度閾值檢測方法失效。而本文算法檢測的疊掩和陰影區域邊界一致性較好，檢測性能更優。

圖 10 河北赤城精度檢測場SAR圖像試驗結果

Fig. 10 The results of the SAR Image of Hebei Chicheng calibration field

圖選項

圖 11 不同方法檢測結果

Fig. 11 The results of different methods

圖選項

為進一步驗證本文算法應用於實測數據的有效性，將檢測出的疊掩和陰影區域在相位解纏過程中進行掩膜。用於計算高程誤差的輔助DEM選自天繪中心測量得到的5 m精度的DEM數據，未掩膜時相位解纏用時23.77 s，DEM生成的高程誤差均值為3.79 m，掩膜後相位解纏用時16.93 s，DEM生成的高程誤差均值為2.74 m，相位解纏速度提升了28.78%，DEM產品精度提升了27.7%。選取圖 12中紅色線段區域進行分析，其對應的經緯度範圍從(115.623 683°E, 41.210 212°N)到(115.666 799°E, 41.179 376°N)，圖 12(a)為未掩膜時DEM生成後的高程誤差值，圖 12(b)為掩膜後DEM生成後的高程誤差值。為便於對比掩膜前後高程誤差變化結果，採用紅色框標記同一區域，可以看出經過掩膜處理後，圖 12中紅色框區域內的高程誤差明顯降低，證明了本文提出的基於超像素分割與多信息融合的疊掩和陰影區域檢測方法的有效性。

圖 12 TH-2實測數據DEM產品高程誤差

Fig. 12 Elevation errors of DEM products based on TH-2 measured data

圖選項

4 結論

試驗結果表明，本文提出的基於多信息融合的超像素檢測算法能夠快速準確地提取SAR圖像中的疊掩和陰影區域，在疊掩區域識別上具有顯著的優勢，在陰影區域識別上，由於低回波能量區域和實際地物遮擋形成的陰影區域在數據特徵上相似，在沒有實際精確DEM輔助的情況下無法準確區分低回波能量區域和實際地物遮擋形成的陰影區域，但是在大部分實際需求中這些區域都會對後續處理造成影響，因此陰影區域的準確區分並不是必要的。本文所提算法檢測的結果可作為相位解纏環節中的掩膜產品，提升相位解纏的速度和DEM產品的精度。同時提取的超像素信息還可用於SAR目標識別、升降軌融合等領域。本文提出的基於AIEM模型的頻域回波仿真算法，實現了星載全極化雙站SAR的原始回波數據仿真，模擬的回波能夠反映SAR圖像的疊掩、陰影、相干斑噪聲等物理特徵，回波仿真結果可用於成像算法驗證和系統設計參數驗證。

作者簡介

第一作者簡介：劉鵬(1995—)，男，碩士，助理研究員，研究方向為SAR/InSAR數據處理。E-mail：lp18392965216@163.com通信作者：李真芳, E-mail：lzf@xidian.edu.cn

初審：張艷玲

覆審：宋啟凡

終審：金 君

資訊