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

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

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

1 漸進式SBAS-InSAR原理

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

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

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

(1)

(2)

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

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

(3)

式中，δX1、δX2分別為X1、X2的改正數(shù)。

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

(4)

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

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

2 實驗驗證

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

2.1 模擬數(shù)據(jù)實驗

以1個像元為例模擬3種形變模型，即線性模型、周期模型和指數(shù)模型。模擬時選用覆蓋北京平原地區(qū)的26景SAR影像參數(shù)，假設前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求解參數(shù)，得到地表形變時間序列，圖3為添加高斯噪聲與未添加噪聲情況下3種形變模型對應的模擬真值、SBAS-InSAR和漸進式SBAS-InSAR獲取的形變時間序列。結(jié)果顯示，2種技術獲取的未知參數(shù)的標準差小于1 mm，表明漸進式SBAS-InSAR與SBAS-InSAR獲得的參數(shù)精度相當。

圖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 真實數(shù)據(jù)實驗

2.2.1 研究區(qū)和數(shù)據(jù)介紹

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

圖4 研究區(qū)位置Fig.4 Location of the study area

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

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

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

3 實驗結(jié)果與分析

3.1 監(jiān)測結(jié)果分析

3.1.1 形變速率結(jié)果分析

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

圖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種技術獲取的形變結(jié)果具有較高的一致性。

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

3.1.2 形變時間序列結(jié)果分析

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

圖8為漸進式SBAS-InSAR與SBAS-InSAR獲取的累積形變量之間的差異直方圖?？梢钥闯?，形變差異絕對值基本小于1 mm，最大形變差異絕對值不超過3 mm，表明漸進式SBAS-InSAR獲取的累積時間序列與SBAS-InSAR的監(jiān)測結(jié)果具有較好的一致性。

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

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

3.2 時間效率分析

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

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

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

圖10 漸進式SBAS-InSAR與SBAS-InSAR所需時間對比Fig.10 Comparison of data processing time between progressive SBAS-InSAR and SBAS-InSAR

4 結(jié) 語

本文分別采用模擬數(shù)據(jù)和真實數(shù)據(jù)對漸進式SBAS-InSAR和SBAS-InSAR進行對比分析，得出如下結(jié)論：

1)在添加高斯噪聲與未添加噪聲的情況下，3種不同形變模型的計算結(jié)果顯示，2種技術獲取的未知參數(shù)的標準差均小于1 mm，表明漸進式SBAS-InSAR與SBAS-InSAR解算精度相近。

2)利用2019-06～2020-04的26景Sentinel-1A數(shù)據(jù)對北京平原地區(qū)的地表形變進行監(jiān)測與分析，發(fā)現(xiàn)2種技術反演得到的地表形變速率在量級和形變區(qū)域分布上均具有較高的一致性，差值的絕對值基本小于0.5 mm/a，最大不超過1.5 mm/a。

3)針對隨時間動態(tài)增加的SAR影像，漸進式SBAS-InSAR獲取形變結(jié)果所需時間僅是SBAS-InSAR的1/2～1/5，可大幅提高地表形變的計算效率。