基于觀測時間配準的升降軌InSAR數(shù)據(jù)二維形變場解算

趙建輝,王 博,閔 林,李 寧

(1. 河南大學(xué) 河南省大數(shù)據(jù)分析與處理重點實驗室,河南 開封 475004；2.河南大學(xué) 河南省智能技術(shù)與應(yīng)用工程技術(shù)研究中心,河南 開封 475004；3.河南大學(xué) 計算機與信息工程學(xué)院,河南 開封 475004；4.河南大學(xué) 信息化管理辦公室,河南 開封 475004)

地質(zhì)災(zāi)害給人類的經(jīng)濟生活帶來了巨大的災(zāi)難，究其原因，絕大部分都是由于地球表面的形變引起的，其中不僅有地震形變、地面沉降、火山運動、冰川漂移以及山體滑坡等自然災(zāi)害，還有由于工程開挖、地下水抽取、爆破、棄土等引發(fā)的人為地質(zhì)災(zāi)害。這些不可逆的地表形變已經(jīng)成為影響區(qū)域經(jīng)濟和社會可持續(xù)發(fā)展的重要因素。

干涉合成孔徑雷達(Interferometric Synthetic Aperture Radar，InSAR)測量技術(shù)憑借其全天時、全天候、觀測范圍大和精度高等優(yōu)勢，廣泛應(yīng)用于各種因素引起的地表形變監(jiān)測中[1-6]。然而，傳統(tǒng)的InSAR技術(shù)存在視線方向模糊的問題，僅能獲取地表雷達視線方向的一維形變場，不能完全反映地表的實際形變狀態(tài)。為了獲取更加真實的地表形變狀況，眾多學(xué)者進行了InSAR二維形變場監(jiān)測研究。2011年，Gray等人采用差分InSAR(Differential InSAR，D-InSAR)技術(shù)獲取升降軌視線向一維形變場，根據(jù)SAR衛(wèi)星觀測的空間特征建立形變解算模型，聯(lián)合兩個一維形變場，解算得到垂直和水平東西方向高精度的二維形變場[7]。2012年，Samsonov等人利用多維小基線(MultidimensionalSmall Baseline Subset InSAR，MSBAS-InSAR)技術(shù)處理不同參數(shù)的SAR數(shù)據(jù)，提取了剛果維龍加火山區(qū)地表二維形變速率場[8]。該方法需要用到較多數(shù)量的升降軌數(shù)據(jù)，并不適用于衛(wèi)星升降軌數(shù)據(jù)較少的地區(qū)。2014年，王艷等人基于永久散射干涉測量(Persistent Scatterer InSAR，PS-InSAR)技術(shù)，分別獲取了同一衛(wèi)星一系列相近時段的升軌和降軌地表雷達視線一維形變速率場，利用二維簡化模型得到了蘇州地區(qū)地表垂直和水平東西方向的二維形變速率場[9]。該方法同樣需要一定數(shù)量的升降軌SAR數(shù)據(jù)。2018年，劉媛媛等人利用相干點目標法(Interferometric Point Target Analysis，IPTA)，獲取山西省青徐和交城地區(qū)不同衛(wèi)星數(shù)據(jù)的地表雷達視線向累積形變量，通過參照先驗信息忽略了該地區(qū)變化較小的南北向形變，聯(lián)合已獲取的累積形變量，解算得到該地區(qū)垂直和水平東西方向的地表二維形變場[10]。目前該方面的相關(guān)研究都是通過聯(lián)合升降軌視線向一維形變場解算地表二維形變場，并未考慮不同軌道一維形變場之間時間不一致的問題。單一軌道的一維形變場是衛(wèi)星在該軌道前后兩次對地觀測時間段內(nèi)地表產(chǎn)生的形變，受到衛(wèi)星軌道自身參數(shù)的限制，單一衛(wèi)星或多個衛(wèi)星都很難在相同的時間做到對同一目標的觀測。在二維形變解算過程中，即使采用了時間較近的升降軌SAR數(shù)據(jù)，也未能解決觀測時間不一致對形變場計算結(jié)果的影響問題。

為了解決上述問題，該文提出了一種升降軌觀測時間配準方法。首先，基于時間序列InSAR技術(shù)，分別獲取升軌和降軌方向的一維累積形變量；然后，以其中一軌SAR時序數(shù)據(jù)的采集時刻為參考，采用數(shù)據(jù)插值技術(shù)，對另一軌SAR時序數(shù)據(jù)獲取的時間累積形變量進行內(nèi)插，得到對應(yīng)參考時刻點的形變量，實現(xiàn)升降軌觀測時間配準；最后，聯(lián)合解算時間配準后的InSAR觀測值，獲取精度更高的二維形變場。

1 升降軌SAR信息反演地表二維形變場

常規(guī)InSAR技術(shù)監(jiān)測地表形變，其觀測值代表了兩景數(shù)據(jù)采集時間內(nèi)地表發(fā)生的形變。該形變量是地表真實形變在雷達視線方向的投影大小，并不能完全反映地表的真實形變狀態(tài)。由于現(xiàn)有在軌雷達具有極地飛行的特點導(dǎo)致了其在南北方向形變監(jiān)測的能力較弱，因此通過忽略雷達視線向形變在南北向的形變分量，聯(lián)合升軌和降軌InSAR觀測值，解算得到垂直向和水平東西向的地表二維形變場。

(1)

(2)

圖1 雷達成像幾何示意圖

2 基于立方插值技術(shù)的升降軌觀測值時間配準

在地表形變監(jiān)測中，受重訪周期等因素的影響，很難獲得持續(xù)的時間序列累積形變量。利用數(shù)據(jù)插值技術(shù)對時間序列累積形變量進行處理，估計缺失時刻的形變量，可以得到完整時間序列的地表形變變化，并對地表形變進行分析[11]。線性插值、三次樣條插值及立方插值等都是地表形變監(jiān)測中常用的補充缺失形變信息的方法。相對于其他插值技術(shù)，經(jīng)過三次樣條插值和立方插值處理后，二者結(jié)果都是一個完整的、隨時間平滑變化的累積形變量。然而，當缺失時刻的數(shù)據(jù)分布不均勻時，三次樣方插值容易出現(xiàn)異常結(jié)果，立方插值則可以有效避免上述情況的發(fā)生[12]。

圖2 升降軌InSAR數(shù)據(jù)二維形變場解算處理流程

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

3.1 實驗區(qū)與數(shù)據(jù)源

本文實驗數(shù)據(jù)采用歐空局Sentinel-1A衛(wèi)星提供的IW(Interferometric Wide Swath)成像模式下的12景時間序列SAR影像，其中包括2016年采集的10景升軌SLC(Single Looking Complex)數(shù)據(jù)和2景降軌SLC數(shù)據(jù)，數(shù)據(jù)采集時間如表1所示。實驗區(qū)域為江蘇省東南部的蘇錫常地區(qū)，其整體地勢較低，影像內(nèi)地形以平原為主，升降軌SAR影像覆蓋范圍如圖3所示。

表1 研究所用升降軌SAR影像數(shù)據(jù)采集時間

3.2 實 驗

由于采集的實驗區(qū)域衛(wèi)星數(shù)據(jù)中升軌SAR影像較多而降軌SAR影像較少，實驗中對升降軌數(shù)據(jù)分別使用了不同的InSAR處理技術(shù)，其中使用小基線集InSAR(Small Baseline Subset InSAR，SBAS-InSAR)技術(shù)對2016-02-26—2016-12-10的10景升軌SAR影像進行處理，使用D-InSAR技術(shù)對2016-06-08和2016-07-14的2景降軌SAR影像進行處理，分別得到研究區(qū)域雷達升軌視線向時間序列累積形變量和降軌視線向形變量。由于時間序列累積形變量是基于日期排序的，不利于下一步的插值處理，因此在插值處理時首先將數(shù)據(jù)采集日期轉(zhuǎn)換成了以天為單位的數(shù)據(jù)采集時間間隔，即以起始日期2016-02-26為基準時間0，將之后每一景數(shù)據(jù)的采集日期都轉(zhuǎn)換為其與起始日期之間的間隔天數(shù)，轉(zhuǎn)換后升軌數(shù)據(jù)采集時間間隔如表2所示。對降軌兩景數(shù)據(jù)的采集日期做相同轉(zhuǎn)換，則其對應(yīng)時間間隔分別為104 d和140 d。將降軌數(shù)據(jù)的這兩個時間間隔作為參考，在10景升軌SAR數(shù)據(jù)所得到的升軌視線向時間序列累積形變量中，用立方插值法對升軌視線向累積形變量進行內(nèi)插，得到與降軌觀測時間對應(yīng)的累積形變量，作差處理后得到104～140 d時間段內(nèi)升軌監(jiān)測得到的形變量，完成時間配準。

圖3 升降軌SAR影像覆蓋范圍

表2 升軌SAR數(shù)據(jù)采集時間間隔

在解算地表二維形變場之前，采用地形校正的方法對升降軌結(jié)果形變基準進行統(tǒng)一。升降軌數(shù)據(jù)處理過程中，兩者采用相同的參考位置進行地形校正，避免基準出現(xiàn)偏差。然后，在相同的形變基準下完成地理編碼，得到升降軌視線向形變量。最后，聯(lián)合時間配準后的兩個軌道的一維形變量，解算獲取研究區(qū)域的視線向二維形變場。在二維形變解算過程中，僅對升降軌中具有相同位置的目標點進行處理。圖4所示為時間配準前后升降軌視線向InSAR觀測值的差異對比情況。

圖4 時間配準前后升降軌InSAR觀測值差異對比

圖4(a)是未經(jīng)時間配準的初始升軌InSAR觀測值，代表升軌2016-06-01—07-19共48 d內(nèi)的視線向形變量，圖4(b)是經(jīng)過時間配準后其對應(yīng)的觀測值，代表升軌2016-06-08—07-14共36 d內(nèi)的視線向形變量，配準前后的觀測值時間相差12 d，差異較大。由于研究區(qū)域?qū)儆诘貏葺^為平原地區(qū)，地下水的過度抽取是該地區(qū)產(chǎn)生地表形變的主要原因，地表整體形變是一個隨時間緩慢變化的過程[13]，時間的差異對InSAR觀測值有較大的影響。InSAR觀測時間差異越大，相應(yīng)的InSAR觀測值差異越大。觀察圖4(a)和圖4(b)，A，B，C和D分別為常州市城區(qū)、無錫市江陰地區(qū)、無錫市城區(qū)、蘇州市城區(qū)。時間配準前，圖4(a)中A，B區(qū)域InSAR觀測值較大，地表形變變化較為明顯，C，D區(qū)域InSAR觀測值較小，地表形變變化相對較小。時間配準后，圖4(b)中A，B區(qū)域InSAR觀測值相對于配準前有所減小，C，D區(qū)域InSAR觀測值幾乎觀察不到變化。通過對時間配準前后的升軌InSAR觀測值之間做差處理，可以發(fā)現(xiàn)整個蘇錫常地區(qū)的InSAR觀測值均出現(xiàn)了一定的變化，差異結(jié)果如圖4(c)所示。

為了進一步對比分析時間配準前后的結(jié)果，對圖4(a)、圖4(b)和圖4(c)中A區(qū)域10個相同點位的數(shù)值進行對比，結(jié)果如圖5所示。升軌InSAR觀測值經(jīng)過時間配準后，形變數(shù)值較時間配準前明顯減小。一系列差值結(jié)果的數(shù)值表明時間配準對升軌InSAR觀測值具有明顯的作用。因此，該方法可以為下一步二維形變解算提供更準確的InSAR觀測值。

圖5 形變數(shù)值對比圖

降軌InSAR觀測值如圖4(d)所示，對比觀察時間配準后的升軌InSAR觀測值和降軌InSAR觀測值，兩者形變位置分布相同，但監(jiān)測結(jié)果量值存在一定差異。根據(jù)二維形變解算式(2)，聯(lián)合配準后的升軌InSAR觀測值和降軌InSAR觀測值，得到蘇錫常地區(qū)的垂直向形變量和水平東西向形變量，如圖6(a)和圖6(b)所示。解算結(jié)果表明，在垂直形變方面，蘇錫常地區(qū)2016-06-07主要沉降區(qū)域分布在常州城區(qū)及無錫江陰地區(qū)，蘇州市城區(qū)出現(xiàn)了整體的地表抬升；常州市城區(qū)的垂直沉降值范圍主要分布在-8～-16 mm；無錫市城區(qū)整體地表形變趨勢較為平緩，形變量-5～+5 mm。在水平東西形變方面，蘇州市、無錫市和常州市三地主要城區(qū)整體呈現(xiàn)輕微水平西向形變，而蘇錫常地區(qū)北部(臨近長江區(qū)域)發(fā)生了水平東向形變，靠近山區(qū)、湖邊以及長江的地區(qū)皆發(fā)生了較為明顯的水平東西向形變。本文實驗反演所得的蘇錫常地區(qū)形變信息與相關(guān)部門公布的測量信息基本一致。實驗結(jié)果表明，所提方法可以更加準確地反演地表二維形變場。

圖6 研究區(qū)域地表二維形變場

4 結(jié)束語

本文基于升降軌InSAR數(shù)據(jù)解算地表垂直方向和水平東西方向的二維形變場，利用數(shù)值插值技術(shù)實現(xiàn)了升降軌模式下InSAR觀測值的時間配準，減少了兩者時間不一致對形變場反演精度的影響。實驗結(jié)果表明，聯(lián)合時間配準后的升降軌InSAR觀測值可以更準確地反演地表二維形變場。

升降軌InSAR觀測值不僅存在時間不一致的問題，也存在空間不一致的問題。由于本文實驗區(qū)域主要為平原地形，其地表形變是長時間緩慢變化的過程，因此本文并未考慮不同地形對升降軌InSAR觀測值的影響。如何解決升降軌InSAR觀測值的空間不一致問題，以進一步提高升降軌InSAR數(shù)據(jù)二維形變場的解算精度，是本文下一步的研究方向。