基于對數域加性信號分解的時序SAR圖像相干斑抑制方法

康 健 童風雨 白雨松 丁 翔 冀騰宇 張 柘

①(蘇州大學電子信息學院 蘇州 215006)

②(西北工業大學數學與統計學院 西安 710072)

③(蘇州市微波成像處理與應用技術重點實驗室 蘇州 215123)

④(蘇州空天信息研究院 蘇州 215123)

⑤(中國科學院空天信息創新研究院 北京 100190)

1 引言

合成孔徑雷達(Synthetic Aperture Radar,SAR)通過主動發射、接收并處理電磁波信號，成為具有全天時、全天候對地觀測能力的遙感技術，在國防、自然資源等領域發揮著重要作用。近年來，新的衛星及信號處理技術的飛速發展極大地促進了SAR對地觀測技術的進步，使得同時滿足高分辨率、大幅寬以及高重訪頻率的SAR圖像獲取成為可能，從而提升了其在高精度、大范圍地表動態監測等方面的能力[1–3]。通過將不同時間獲取的SAR圖像進行配準，得到的SAR圖像時間序列能提供被觀測地區在時間維度的變化信息，進而能更好地服務于農業監測、災后評估以及城市規劃等方面[4]。

SAR成像過程中的相干處理特性使得圖像中像素幅度值呈現出隨機擾動的現象，這主要是雷達波與像素點內各個散射點相互作用產生的回波相參疊加所導致的。與常見的加性高斯噪聲不同，這種相干斑噪聲是一種乘性噪聲，其非平穩特性對SAR圖像的后續解譯構成了嚴重挑戰[5,6]。一直以來，SAR圖像的相干斑抑制是SAR圖像處理的重要任務之一，早期的去噪方法采用簡單的低通濾波器對單視SAR圖像進行濾波，在噪聲被抑制的同時，圖像細節信息也存在一定程度的丟失。為了解決此問題，文獻[7,8]均利用空間域自適應濾波的方法，在降噪的同時，可以充分保留圖像的細節特征。文獻[9]引入了全變分(Total Variation,TV)正則化項，且利用對數變換將噪聲的乘性模型轉化為滿足加性模型的數據擬合項，進而提出了相干斑抑制的優化模型并利用最優化方法得到去噪結果。近幾年，隨著深度學習方法在圖像處理領域的廣泛應用，一些工作也將其用來抑制SAR圖像相干斑，并取得了良好的效果[10–12]。其他單通道SAR圖像抑制方法詳見綜述類文獻[6,13,14]。

近年來，世界各國均著眼于SAR衛星星座的建設，如歐空局的哨兵1號(Sentinel-1)衛星星座、中國的高分3號(Gaofen-3)衛星星座，使得具有全球尺度、高重訪頻率的SAR圖像處理成為領域內研究的重點。空間域配準得到的具有高時間分辨率的SAR圖像時間序列能充分地反映出觀測時間段內地物變化信息。針對時序SAR圖像的相干斑抑制，國內外學者也做了很多相關工作。文獻[15–17]將單通道SAR圖像搜索同質像素點或者圖像塊進行去噪的思想拓展到時序SAR圖像中，即在時空方向上均進行同質像素點或圖像塊的搜索，進而再進行濾波。文獻[18]提出的多基線干涉SAR技術(SqueeSAR)能自適應地將同質分布式散射點進行選取并聯合濾波，進而同時實現對多時相SAR幅度及相位進行去噪。文獻[19]提出的多時相SAR圖像三維塊匹配(Multitemporal SAR-BM3D,MSAR-BM3D)方法利用時間維度上的冗余信息以及變換域協同濾波將單通道SAR圖像三維塊匹配(SAR Block Matching 3D algorithm,SAR-BM3D)方法[20]進行拓展，使其能實現對時序SAR圖像的相干斑抑制效果。文獻[21]提出的基于比值的多時相SAR圖像去噪(Ratio-based Multitemporal SAR,RABASAR)方法通過利用SAR時間序列圖像與其在時間維度平均得到的均值圖像計算對數比，使SAR圖像相干斑的空間非平穩性能被很好地抑制，從而能進一步提升時序SAR圖像的去噪效果。

然而，除了空間域SAR圖像中原有的相干斑噪聲帶來的挑戰之外，時間序列中存在的變化信息，如出現人造目標、建筑物的結構變化等，均能使得觀測區域在時間維度上出現異常的散射強度，如圖1所示。由于這些沿時間維度突變信號來源于觀測區域地物變化，特別是人造目標，現有的大部分針對SAR圖像時間序列進行相干斑抑制的方法并沒有對此信號成分進行建模，從而導致去噪結果中含有人為引入的噪點或者出現“虛假散射點”現象。

圖1 含有沿時間維度突變信號的SAR幅度時間序列Fig.1 SAR amplitude time series with outliers along the temporal dimension

另外，在一些應用領域，如目標檢測等，這些沿時間維度突變信號通常反映了一些人造目標在觀測時間內的變化情況，也能為后續的解譯工作提供有價值的信息。為此，Baier等人[22]提出了一種非局部低秩信號分解的方法對時序SAR圖像進行相干斑抑制(Robust Nonlocal Low-Rank SAR Time Series Despeckling Considering Speckle Correlation by Total Variation Regularization,DespecKS-NLLRTV)，分別將觀測數據分解為低秩信號成分、時間相關的相干成分以及稀疏的沿時間維度突變信號成分，在相干斑抑制的同時也能降低突變信號對去噪結果產生的影響。雖然此方法在實際數據中的效果明顯，但是由于其加性信號分解是在原始數據域中進行的，并沒有將信號成分與相干斑成分進行充分解耦，而且相干斑成分與突變信號成分均被描述為稀疏項并由 L1范數進行建模，也會導致相干斑成分與突變信號成分存在耦合現象，從而導致信號分解結果不準確，并不能真實地反映觀測時間內的地物變化所導致的散射強度突變的情況。此外，DespecKS-NLLRTV在對TV鄰近算子進行求解時采用了Chambolle-Pock算法[23]，其運算效率不高，進一步限制了該方法在實際數據中的應用。

為了對時序SAR圖像進行相干斑抑制的同時，準確地提取出觀測時間段內產生的沿時間維度突變散射點，本文提出了一種基于對數域加性信號分解的方法，首先，通過對對數域SAR強度值進行統計建模，得到時序SAR數據的對數似然，再引入針對不同成分的先驗知識，提出相應的正則項，進而得到各個成分的后驗概率，通過最大后驗概率(Maximum A Posterior,MAP)的估計方法，再得出問題對應的優化模型，最后利用交替方向乘子法(Alternating Direction Method of Multipliers,ADMM)[24]進行求解優化。所提方法在仿真以及實際數據中均驗證了效果，通過對比試驗分析了所提方法在去噪性能以及突變信號成分分離效果上的提升。

2 問題建模與求解方法

本文提出的方法總體框架如圖2所示，為了保持時序SAR圖像中永久散射點(Persistent Scatterer,PS)的散射特性，首先計算振幅離差指數(Amplitude Dispersion Index,ADI)[25]剔除可能的永久散射點，對于每個非永久散射點，再通過Kolmogorov-Smirnov (K-S)檢測[18]將非局部搜索窗口中的同質像素點進行篩選，進而將所有同質點與其排列并構成待分解矩陣，應用所提出的對數域信號分解方法得到相干斑抑制成分以及沿時間維度突變信號成分，最后采用聚合求均值策略得到最后的去噪以及突變信號結果。

圖2 本文所提方法流程圖Fig.2 Flow chart of the proposed method

2.1 問題建模

SAR圖像的強度值服從以下Gamma分布[9,26]：

其中，I ∈R+表示像素強度值，R ∈R+為待恢復的信號強度值，L>0為視數個數，Γ(·)為Gamma函數。在乘性噪聲模型下，強度I表示為

其中，n為均值，為E(n)=1的獨立同分布的隨機變量，進而得到強度的均值與方差為

從式(4)可以看出，噪聲在空間方向上是非平穩的，并與信號真實強度值有關。為了方便處理，通常采用對數函數將乘性噪聲模型進行轉化，得到滿足加性噪聲模型的Fisher-Tippett分布[26]：

其中，g=lg(I),z=lg(R)。根據上述對數域中的強度似然模型，在MAP框架下，待恢復的信號可以通過求解以下模型來得到：

其中，lgpg(g|z)為對數域中的信號對數似然，R(z)為待恢復信號的先驗知識對應的正則化項，λ為正則化項參數。對于時序SAR圖像，空間像素點中的時間序列對數域強度值服從以下對數似然：

其中，t為時間索引，C 為常量，在對變量zt進行優化時，可以將其忽略，并不影響優化結果。

上述MAP問題需要對對數域時間序列強度值的先驗知識進行建模，本文主要考慮以下先驗：

(1) 非局部自相似性：在一定觀測范圍內，存在地物目標尺寸遠大于SAR圖像分辨率或者不同目標屬于同種地物類別等情況，使得區域內的目標散射特性有相似性，進而可以利用非局部搜索方法將一些同質像素點對應的時間序列找出并進行聯合處理。為此，本文利用K-S檢測方法將一定搜索范圍內(例如，21×21窗口范圍)的同質點進行搜索選取，再將候選時間序列構成觀測矩陣。

(2) 平滑性：由于觀測矩陣是由同質像素點對應的時間序列構成的，其在時間以及空間維度上(表示為矩陣的行和列)具有平滑特性，即散射強度變化不大，本文采用TV范數來對此進行建模。

(3) 稀疏性：對于時序信號中存在的沿時間維度突變信號，其本身具有稀疏特性，本文利用 L1范數對其進行刻畫。

根據上述似然以及正則化項，本文提出了如下模型對時序SAR圖像進行相干斑抑制并對觀測數據中的穩定信號及突變信號成分進行分離：

其中，元素gti所構成的矩陣是將K-S檢測方法選取到的同質像素點對應的時序列向量進行按行排列得到的，其時間維度及空間維度的索引分別為t,i。式(8)中X為相干斑抑制后的信號成分，S為沿時間維度突變信號成分，α,β分別為不同正則化項對應的參數。給出矩陣A，其TV范數∥A∥TV定義為

其中，i1,i2分別為矩陣列和行的索引。

2.2 求解方法

為了有效求解所提模型式(8)，本文首先利用變量替換及分解方法[9]將上述模型進行轉化，得到：

其中，D(·):=[Dt(·);Di(·)]表示二維1階導數算子。再將有約束的優化模型轉化為無約束優化模型，得到增廣拉格朗日函數：

其中，〈·〉表示為內積算子，∥·∥F為矩陣F-范數，T1,T2為優化引入的輔助變量，μ為懲罰項對應的參數。進而利用ADMM方法對其進行優化求解，求解過程如下。

2.2.1Z子問題

類似于文獻[26]，上述問題可以通過牛頓法(Newton’s method)進行逐像素多次迭代求解：

2.2.2X子問題

可以通過計算式(14)對于變量X的梯度，并且令梯度為0，得到線性方程：

其中，D?(·)表示D(·)的伴隨算子，利用頻域快速求解方法[27]得到式(14)的解。

2.2.3Y,S 子問題

式(16)與式(17)的解可以由軟閾值(Soft Thresholding)[28]算子得出。

2.2.4 乘子T1,2更新

按照式(18)，式(19)子問題求解規則，在迭代一定次數并滿足收斂條件后，得到的與為當前像素點及其同質點時間序列的相干斑抑制矩陣和突變信號矩陣。

由于上述分解過程是在對數域中進行的，需要將這兩種成分轉換到原始信號域，令為去噪后的SAR圖像像素強度值，根據所提方法的分解模型可得到：

則在原始信號域中，exp()為去噪后的穩定信號成分，exp()(exp()-1)為沿時間維度突變信號成分。

2.3 永久散射點剔除

2.2 節所提方法是針對每個非永久散射點進行的，由于永久散射點幾乎不受相干斑噪聲的影響，這些散射點的時序散射信息需要被保留。本文采用ADI方法對其進行剔除，通過計算SAR時間序列幅度值的均值與標準差之間的比值，即mA/σA，并將大于一定閾值的像素點進行剔除。

對于每個非永久散射點，所提方法均會找到相對應的同質點時序信號，并對所構成的矩陣進行分解，在遍歷所有像素點后，再利用求平均的策略將每個像素點對應的去噪成分以及沿時間維度突變信號成分進行聚合，從而得到最終的相干斑抑制結果以及突變信號分離結果。

3 仿真與實測數據實驗

3.1 仿真數據實驗結果

本文采用的仿真數據與文獻[22]一致，其幅度圖像與其隨時間變化如圖3所示。產生均值為1的高斯隨機變量以及值為0.3的指數相參損失用來模擬SAR圖像中的相干斑效應。共生成24幅圖像作為仿真的時間序列。此外，分別隨機選取了0.5%以及1.0%比例的SAR像素，用服從Rice分布的隨機數將其替換，用來模擬沿時間維度突變信號干擾現象。利用K-S檢測選擇同質點的閾值設置為0.25，正則化項參數α,β分別選取為5.0,10.0。為了驗證所提方法的有效性，本文選擇了以下4種時序SAR圖像去噪方法進行對比：(1) SqueeSAR[18];(2) MSAR-BM3D[19];(3) RABASAR[21];(4) DespecKS-NLLRTV[22]，并采用峰值信噪比(Peak Signal-to-Noise Ratio,PSNR)、平均結構相似性指標(Mean Structure Similarity Index Measure,MSSIM)以及平均互相關(Mean Cross Correlation,MCC)來定量驗證方法在相干斑抑制上的效果。

圖3 本文采用的仿真數據Fig.3 Simulation data used in this paper

圖4(a)—圖4(f)、圖4(g)—圖4(l)分別展示了所有對比方法在沒有沿時間維度突變信號以及在0.5%像素的突變信號干擾下得到的相干斑抑制結果。相應的定量結果如表1所示。在沒有突變信號的情況下，現有的方法如MSAR-BM3D和RABASAR能取得比DespecKS-NLLRTV和本文方法更好的相干斑抑制效果。其中一個原因是MSAR-BM3D和RABASAR均利用了相似圖像塊對噪聲進行抑制，這種方法在例子中給出的同質區域較多的場景中能更好地在空間維度對信號進行平滑處理，然而，在MSSIM指標下，所提出的方法能取得比MSARBM3D更好的效果，主要原因是基于同質像素點的濾波方法能更好地保持場景的空間結構信息。當突變信號的數量開始增加時，DespecKS-NLLRTV和本文方法均能魯棒地對相干斑噪聲進行抑制，其恢復效果很大程度上并不受突變信號的數量影響。由于其他方法并沒有考慮時序數據中存在的沿時間維度突變信號，其圖像恢復效果隨著突變信號的增多而變差。與本文所提出的對數域信號分解方法相比，直接在原始信號域進行信號分解的方法DespecKS-NLLRTV在PSNR等指標上有大約3 dB的差距。由于原始信號域中的信號與相干斑噪聲成分存在耦合，直接對觀測得到的數據進行信號分解并不能消除耦合產生的影響，使得信號的恢復結果并不理想。此外，對于非局部自相似性構建的時序矩陣，DespecKS-NLLRTV引入了低秩先驗進行約束，這會導致恢復得到的時序信號能量有所丟失，也會使得PSNR評價指標降低。與此相比，本文所提出的方法充分考慮到了相干斑噪聲的統計特性，并在對數域中將信號成分與其他成分進行分解，能在噪聲抑制的同時，很好地保持恢復信號的原始結構。由于本文方法并沒有采用低秩約束，并不會導致恢復出的時序信號能量丟失，并且可以大幅提升算法的計算效率。值得注意的是，在沿時間維度突變信號數量較少或不存在的情況下，所提出算法的PSNR值反而較低，而隨著突變信號數量增加，PSNR有所增加。其原因在于，沿時間維度突變信號的數量增加使得觀測信號的本質特性更加符合所提方法的信號分解模型，從而使得分解得到的穩定信號成分更加接近于真值。如果不存在突變信號，仍然會有一部分信號成分被分解成為稀疏信號項，反而會影響相干斑抑制得到的信號成分恢復精度。

表1 仿真數據相干斑抑制結果的定量分析Tab.1 Quantitative analysis of speckle suppression results in simulation data

圖4 仿真數據相干斑抑制結果Fig.4 Speckle suppression results of the simulated data

此外，為了測試提出方法對于不同參數設置下的敏感程度，以PSNR值為例，圖5給出了在不同α以及β值下的時序圖像恢復效果。可以看出，當α>1,β >2時，所提方法的恢復效果受參數值變化敏感度不高，而且當突變信號增多時，給予 L1范數項更大的權重可以使得恢復效果有進一步提升。為了進一步驗證所提方法受圖像數目的影響，分別從原始的24張圖像中選取3,5,10,15,20張圖像構成新的時序數據，用所提方法對其進行相干斑抑制，相應的超參數設置保持不變，即α=5.0,β=10.0，進而計算得到結果的PSNR值如圖6所示。可以看出，隨著圖像數目的增加，所提方法的相干斑抑制結果也隨之提升。如果圖像數目過少，會影響K-S檢測方法對于同質點選取的精度，進而使得所構成的時序矩陣自相似性較差，導致所提方法的去噪效果降低。因此，所提出的方法適用于圖像數目較多(多于20張)的時序SAR數據去噪。此外，盡管實驗中突變信號的數目增加，所提方法并不受到突變信號數目不同而產生效果上的影響，從而進一步驗證了所提方法的魯棒性。

圖5 本文方法的參數敏感性分析Fig.5 Parameter sensitivity analyses of the proposed method

3.2 仿實測數據實驗結果

本文分別選取了上海近海以及浦東機場區域的垂直極化干涉寬幅哨兵1號(Sentinel-1 IW VV)數據，兩組時間序列獲取時間段均為2019-12-02到2021-02-06，每組包含36張圖像。其距離向及方位向分辨率分別約為2.33 m和13.94 m。所提方法的參數設置與3.1節仿真數據一致。

圖7展示了所有對比方法分別在2020-09-15以及2021-02-06獲取的SAR圖像上的相干斑抑制結果。可以看出，近海區域的艦船產生的強散射點構成了時間序列中的沿時間維度突變信號。傳統的方法如SqueeSAR對于艦船區域的去噪效果差，其平均結果導致了此區域產生模糊的亮斑，對于MSARBM3D，可以發現艦船附近區域的去噪結果存在人為引入的噪聲，使得其并不能真實反映出海面散射特性。盡管RABASAR能對大部分同質區域進行很好的相干斑抑制，但對于艦船目標，其散射結果易出現模糊現象，從而損失了真實的船體結構信息。相比于上述幾種方法，采用信號分解思想的DespecKS-NLLRTV和本文方法均能通過沿時間維度突變信號分離的策略降低其對于相干斑抑制的影響。相比于DespecKS-NLLRTV，本文方法更好地利用了SAR圖像的統計特性，減小了直接在原始圖像域中做分解存在的信號與相干斑噪聲的耦合效應。另外，由于DespecKS-NLLRTV中的模型對于突變信號以及相干斑成分均采用L1范數進行建模，會導致兩種成分分離結果互相干擾。相比之下，本文方法能更加準確地對突變信號進行提取。為定量評估不同方法的相干斑抑制結果，我們分別在兩次觀測場景中選取了同質區域A1和A2，并根據所有方法得到的去噪結果計算相應區域內的等效視數(Equivalent Number of Looks,ENL)。從表2可以看出，所提方法可以實現良好的去噪效果，可以對于同質區域內的噪聲進行有效抑制。雖然MSAR-BM3D在這些同質區域內取得了最好的等效視數，但其在散射特性復雜的區域并不能得到很好的去噪結果，尤其在具有沿時間維度突變散射點附近，其去噪結果通常會引入人為噪聲干擾。與上述分析一致，相比于原始信號域信號分解方法DespecKS-NLLRTV，本文提出的信號分解策略能在同質區域內獲得更高的ENL。此外，針對得到的突變信號圖像，本文進一步計算其香農熵，從而定量地反映了圖像的復雜程度。表3列出了兩種方法得到的突變信號圖像所計算的香農熵，可以看出本文方法提取得到的突變信號圖像質量較高，而DespecKS-NLLRTV并不能很好地將相干斑噪聲與突變信號成分進行有效的分離，使得提取得到的突變信號中含有噪聲干擾，從而導致了其圖像熵值較高。

表2 4塊同質區域計算得到的等效視數Tab.2 The equivalent apparent number calculated by four homogeneous regions

表3 沿時間維度突變信號的熵值分析Tab.3 Entropy analysis of the outliers along the temporal dimension

圖7 上海近海區域哨兵1號數據的相干斑抑制結果Fig.7 Speckle suppression results for Sentinel-1 data in the Shanghai offshore region

圖8給出了所有方法在上海浦東機場區域2019-12-02以及2020-09-27獲取到的SAR圖像相干斑抑制結果。場景中由于飛機目標的出現(放大區域)，其強后向散射特性構成了機場區域內的SAR像素強度突變，使其成為時間序列中的沿時間維度突變信號。在這些區域內，SqueeSAR方法中的簡單空間平均模糊掉了強散射目標，從而使得去噪結果中含有虛假亮斑。盡管MSAR-BM3D和RABASAR能對大部分區域進行有效的相干斑抑制，但由于這些方法沒有特別考慮時序圖像中的沿時間維度突變信號，使得去噪結果中包含人為引入的噪點或者導致突變信號的目標結構被過度平滑。針對選取的同質區域A3和A4，所提方法均能在取得良好去噪效果的同時，實現對于突變目標的提取，并在所提取的圖像中完整地保留目標的結構。表3得到的突變信號圖像熵值計算結果進一步說明了兩種分解方法對于異常目標分離效果差異。

圖8 上海浦東機場區域哨兵1號數據的相干斑抑制結果Fig.8 Speckle suppression results for Sentinel-1 data in Shanghai pudong airport region

為了驗證不同方法對于圖像細節的保持特性，分別選取了近海區域(如圖9所示)以及海上艦船區域(如圖10所示)，并利用邊緣檢測技術對所有算法的去噪結果進行圖像結構信息提取。可以看出，本文提出的方法能很好地保持近海區域的地物的結構信息，并沒有因相干斑噪聲抑制導致圖像細節信息丟失。與之對比，RABASAR等方法并沒有很好地保持圖像的結構信息，從而使得邊緣提取結果并不準確。在艦船區域結果中，可以看出所提方法很好地將艦船與背景噪聲進行分離，使得艦船邊緣提取結果并沒有受到背景噪聲的干擾。與之相比，RABASAR等方法邊緣提取結果模糊，不能準確地反映目標的結構特性。

圖9 近海區域圖像邊緣提取結果對比Fig.9 Comparison of image edge extraction results in offshore areas

圖10 海上艦船區域圖像邊緣提取結果對比Fig.10 Comparison of image edge extraction results of ship region

4 結語

本文針對時序SAR圖像中的相干斑噪聲，提出了一種基于對數域加性信號分解技術的相干斑抑制方法，所提方法不僅能很好地對相干斑進行抑制，并且能對時間序列中的穩定信號成分以及沿時間維度突變信號成分進行有效分離，得到的突變信號能展示出觀測時間內由于人造目標出現或地物變化產生的強散射現象，為后續的圖像解譯提供了有效的參考數據。為了驗證方法的有效性，本文選取了現有的4種主流時序SAR圖像相干斑抑制方法進行對比，并且用PSNR等幾種指標定量地分析了不同方法在仿真數據上的去噪效果。此外，本文選取了上海近海以及浦東機場區域的哨兵1號時序SAR數據，并將上述方法所得結果進一步進行對比分析，均能驗證所提方法的有效性。在后續工作中，會進一步驗證本文所提方法在觀測角度大幅變化情況下獲取的SAR時序數據上的相干斑抑制效果。