利用漸進式SBAS-InSAR監測分析北京平原地區地表形變

2023-02-28 05:24:46劉媛媛夏元平趙振宇

大地測量與地球動力學 2023年3期

劉媛媛晏霞夏元平,3 趙振宇

1 東華理工大學江西省數字國土重點實驗室，南昌市廣蘭大道418號，330013 2 東華理工大學測繪工程學院，南昌市廣蘭大道418號，330013 3 自然資源部環鄱陽湖區域礦山環境監測與治理重點實驗室，南昌市廣蘭大道418號，330013

傳統SBAS-InSAR技術通過同時處理覆蓋研究區一定時間段內的所有SAR影像來獲取該時段內的地表形變結果[1]。若新增SAR影像，則需將所有影像重新處理以獲取新增SAR影像時刻的地表形變信息，難以滿足高效、持續、動態的監測要求。為解決這一難題，許多學者提出各種解決辦法[2-5]，取得了較好效果。本文提出漸進式SBAS-InSAR技術，即在傳統SBAS-InSAR技術獲取的解算結果基礎上，將較新的少部分已有SAR影像結合新增SAR影像進行增量解算，以得到所有影像獲取時刻的地表形變結果，并在絕大部分已有數據形變參數和已有數據參數矩陣不變的情況下，實現高效、持續、精確的地表形變監測。以2019-06～2020-04覆蓋北京平原地區的26景Sentinel-1A SAR影像為例，利用漸進式SBAS-InSAR進行實驗，并從精度和計算效率2個方面與SBAS-InSAR地表形變監測結果進行對比。

1 漸進式SBAS-InSAR原理

如圖1所示，漸進式SBAS-InSAR分為2個階段：1)利用SBAS-InSAR求解已有SAR影像數據集的地表形變時間序列(階段1)；2)以階段1獲取的已有SAR影像數據集的形變參數為基礎，結合已有SAR影像數據集中較新部分SAR影像與新增SAR影像進行高相干點選取、相位解纏、大氣延遲、軌道誤差和高程誤差糾正，獲取新增干涉對的解纏相位，并以此作為觀測值進行序貫平差解算，得到更新后全部影像時刻的地表形變時間序列(階段2)。

圖1 漸進式SBAS-InSAR流程Fig.1 Flow chart of progressive SBAS-InSAR

首先，采用SBAS-InSAR對已有SAR影像數據集進行處理，獲取對應地表形變時間序列結果；在此基礎上，將已有SAR影像數據集中較新的部分(tn1時刻以后，共包含n景)SAR影像與新增SAR影像數據組合進行差分干涉處理，若研究區內新增N1景SAR影像(對應時間段為tN+1-tN+N1)，在相同時空基線閾值下形成M1個差分干涉對，以去除大氣延遲、軌道誤差后的解纏相位作為觀測值，則可構建如下觀測方程：

(1)

(2)

式中，A、B1、B2、C分別為系數矩陣，由0或相鄰兩景SAR影像的時間間隔構成(如12、24等)；X1、X2、X3為未知參數，由相鄰兩景SAR數據間的形變相位速率構成。

若將已有SAR數據和新增SAR數據的解算認為是2次平差過程，則從形式上看，式(1)和式(2)與測量平差中的序貫平差方法一致，故可將SBAS-InSAR求解已有SAR影像集的形變相位速率視為序貫平差的階段1，將已有SAR影像集中較新的部分SAR影像集與新增SAR影像的形變相位速率視為序貫平差的階段2，即可利用序貫平差思想更新所有SAR影像獲取時刻的形變相位速率。

(3)

式中，δX1、δX2分別為X1、X2的改正數。

假設各觀測時刻獲取的SAR影像相互獨立，且權矩陣為單位矩陣，利用序貫平差思想可得[6]：

(4)

式中，Q為SBAS-InSAR求解中存儲的協方差矩陣。

綜合式(3)和式(4)即可得到新增N1景SAR影像后各相鄰SAR影像獲取時刻的形變相位速率，在此基礎上，通過數值積分可得到更新后的形變時間序列。

2 實驗驗證

采用模擬數據和真實數據，從解算精度和計算效率2個方面驗證漸進式SBAS-InSAR的可行性與有效性。

2.1 模擬數據實驗

以1個像元為例模擬3種形變模型，即線性模型、周期模型和指數模型。模擬時選用覆蓋北京平原地區的26景SAR影像參數，假設前16景(2019-06-03～12-24)SAR影像為存檔SAR影像，第17～26景(2020-01-05～04-22)為新增SAR影像，設置時空基線閾值為60 d與120 m，共形成112個干涉對，其時空基線組成如圖2所示。此外，對3種形變模型分別添加均值為0 mm、方差為1.5 mm的高斯噪聲，并進行100次蒙特卡洛隨機模擬。分別利用SBAS-InSAR與漸進式SBAS-InSAR求解參數，得到地表形變時間序列，圖3為添加高斯噪聲與未添加噪聲情況下3種形變模型對應的模擬真值、SBAS-InSAR和漸進式SBAS-InSAR獲取的形變時間序列。結果顯示，2種技術獲取的未知參數的標準差小于1 mm，表明漸進式SBAS-InSAR與SBAS-InSAR獲得的參數精度相當。

圖2 已有SAR影像和已有SAR影像中較新部分SAR影像與新增SAR影像的時空基線分布Fig.2 Spatiotemporal baseline distribution of existing SAR images and newer SAR images in existing SAR images and new SAR images

圖3 3種形變模型下SBAS-InSAR和漸進式SBAS-InSAR的形變時間序列Fig.3 Deformation time series derived from SBAS-InSAR and progressive SBAS-InSAR with three different deformation models

2.2 真實數據實驗

2.2.1 研究區和數據介紹

選取北京通州區、朝陽區等地面沉降嚴重的地區進行研究，覆蓋范圍約38 km×42 km(圖4)。研究區位于第四系凹陷區，受地下水水位降低及地下水不斷超采的影響存在明顯的地表形變[7]。利用SBAS-InSAR和漸進式SBAS-InSAR處理與分析2019-06～2020-04共26景C波段的升軌Sentinel-1A SAR影像。為去除地形和軌道誤差對干涉相位的影響，同時使用歐空局提供的30 m分辨率SRTM1-DEM高程數據和精密軌道數據(定位精度優于5 cm)進行校正。

圖4 研究區位置Fig.4 Location of the study area

2.2.2 數據處理

假設已有SAR影像集包含16景(2019-06-03～12-24)SAR影像，設置時空基線閾值為60 d與120 m，得到54對干涉對(圖2黑色實線)。影像配準和重采樣后，結合外部DEM及精密軌道數據去除地形及軌道誤差，對高程誤差進行建模和估算以將其扣除，采用時空域的高低通濾波去除大氣延遲效應，最后利用SBAS-InSAR求解出相應時刻的形變相位值。在此基礎上，逐漸新增1～10景(2020-01-05～04-22)SAR影像，并將已有SAR影像中較新部分SAR影像納入新增SAR影像中(圖2紅色虛線)，采用漸進式SBAS-InSAR進行處理，獲取全部影像的形變時間序列。

為分析漸進式SBAS-InSAR的精度與時間效率，同時利用SBAS-InSAR對所有SAR影像進行解算，獲取其地表形變時間序列和所需時間。

3 實驗結果與分析

3.1 監測結果分析

3.1.1 形變速率結果分析

圖5為新增10景SAR影像時SBAS-InSAR與漸進式SBAS-InSAR分別反演得到的2019-06～2020-04北京平原地區沿雷達視線向的地表形變速率場。可以看出，研究區形變空間分布不均勻，主要集中在朝陽區與通州區交界地帶，包括金盞鄉、樓梓莊鄉、徐辛莊鎮、黑莊戶鄉、臺湖鄉、梨園、平西府鎮、燕丹鄉和李橋鎮等地，最大形變速率約為-110 mm/a。結合圖5(a)和5(b)可知，無論是形變量級還是形變分布區域，2種技術反演得到的形變結果均具有較高的一致性。

圖5 年平均形變速率Fig.5 Annual average deformation rate

為更直觀地顯示本文方法的可靠性，對漸進式SBAS-InSAR與SBAS-InSAR獲取的年平均形變速率作差(圖6)。可以看出，二者差值的平均值和標準差分別為0.01 mm/a和0.1 mm/a，且2種技術獲取的形變速率之差的絕對值基本小于0.5 mm/a，最大不超過1.5 mm/a，表明2種技術獲取的形變結果具有較高的一致性。

圖6 漸進式SBAS-InSAR與SBAS-InSAR年平均形變速率對比Fig.6 Difference of annual average deformation rate between progressive SBAS-InSAR and SBAS-InSAR

3.1.2 形變時間序列結果分析

圖7為利用漸進式SBAS-InSAR獲取的形變時間序列。可以看出：1)2019-06～2020-04研究區主要形變位于朝陽區與通州區的交界處，最大累積形變量約為-106 mm。2)研究區西部及東北部表現為輕微的地面抬升現象。3)隨著監測時間的累積，金盞鄉、樓梓莊鄉、徐辛莊鎮、黑莊戶鄉、臺湖鄉和梨園形變區的累積形變量不斷增加，形變范圍亦不斷擴大；平西府鎮、燕丹鄉和李橋鎮形變區的累積形變量有所增加，形變范圍基本不變。

圖8為漸進式SBAS-InSAR與SBAS-InSAR獲取的累積形變量之間的差異直方圖。可以看出，形變差異絕對值基本小于1 mm，最大形變差異絕對值不超過3 mm，表明漸進式SBAS-InSAR獲取的累積時間序列與SBAS-InSAR的監測結果具有較好的一致性。

圖8 SBAS-InSAR與漸進式SBAS-InSAR時間序列累積形變量之差直方圖Fig.8 Histograms of the difference of cumulative deformation time series between SBAS-InSAR and progressive SBAS-InSAR

此外，選取3個特征點(A點位于快速形變區、B點位于較快速形變區、C點位于地表穩定區，位置見圖5)對2種技術獲取的時間序列進行對比分析，結果如圖9所示。可以看出，3個特征點的累積形變量差值均小于2 mm，可忽略不計，表明漸進式SBAS-InSAR與SBAS-InSAR獲取的累積形變結果高度吻合。

圖9 3個特征點的形變時間序列Fig.9 Deformation time series of three feature points

3.2 時間效率分析

為便于對比分析，采用相同的準則選取漸進式SBAS-InSAR與SBAS-InSAR的高相干點目標，故2種技術處理的高相干點數目相同。以任意一個高相干點為例，分析漸進式SBAS-InSAR與SBAS-InSAR求解所需時間。對于單個高相干點，在相同時空基線閾值下，當新增N1景SAR影像后，SBAS-InSAR需要對M+M1(M為已有數據形成的干涉對數量，M1為新增SAR影像和已有SAR影像中較新部分形成的干涉對數量)幅干涉對去除地形、糾正大氣效應和軌道誤差等，而漸進式SBAS-InSAR只需對M1幅干涉圖進行處理。從理論上講，當已有SAR影像很多而新增SAR影像不多時，M1?M，因此漸進式SBAS-InSAR比SBAS-InSAR具有更高的效率。

求解形變相位值時，已知形變相位值的求解主要涉及線性方程組(式(2)和式(3))的最小二乘法計算，計算復雜度與待求解未知數的個數有關，未知數越多，計算越復雜。SBAS-InSAR需要求解N+N1景SAR影像對應的地表形變速率增量，而漸進式SBAS-InSAR則只針對新增的N1景SAR影像對應的形變速率增量進行解算，并更新已有N景SAR影像對應的地表形變速率增量，計算復雜度明顯降低，計算所需存儲空間更小，時效性更高。

圖10為在16景SAR影像的基礎上逐漸新增1～10景SAR影像分別進行漸進式SBAS-InSAR與SBAS-InSAR處理的時間成本。數據處理的計算機配置為Intel Core i7(2.60 GHz)，內存32 G。從圖10可以看出，執行配準、大氣和軌道誤差校正等步驟時，SBAS-InSAR的處理時間約為漸進式SBAS-InSAR的2～5倍。在形變參數反演解算階段，漸進式SBAS-InSAR一方面要引入和使用SBAS-InSAR對已有SAR影像求解的協方差矩陣，另外還要對已解算的結果進行修正與更新，故無法像前一階段那樣明顯地縮短運行時間。但隨著SAR影像逐漸增加，漸進式SBAS-InSAR的時間效率相較于SBAS-InSAR仍具有明顯的改進，在增加第10景SAR影像時，其形變解算時間提高約66%。綜上，在處理隨時間動態增加的SAR影像時，漸進式SBAS-InSAR明顯具有更高的時間效率。