時(shí)序InSAR大范圍地表形變監(jiān)測(cè)的關(guān)鍵問(wèn)題分析

劉媛媛，晏 霞

(1. 東華理工大學(xué) 測(cè)繪工程學(xué)院, 江西 南昌 330013; 2. 長(zhǎng)安大學(xué) 地質(zhì)工程與測(cè)繪學(xué)院, 陜西 西安 710054)

城市工程建設(shè)以及社會(huì)經(jīng)濟(jì)的快速發(fā)展導(dǎo)致人口密集程度越來(lái)越高，由此引發(fā)的地面沉降波及區(qū)域逐年擴(kuò)大，影響范圍高達(dá)數(shù)萬(wàn)平方千米，并且具有成片分布的空間分布特征(何慶成等, 2006)。

大范圍(數(shù)萬(wàn)平方千米)地表形變監(jiān)測(cè)要求監(jiān)測(cè)手段既滿足給定的精度指標(biāo)，又具有大范圍覆蓋和同步測(cè)量的能力。傳統(tǒng)的地面測(cè)量技術(shù)(如GNSS、精密水準(zhǔn)測(cè)量等)布設(shè)的監(jiān)測(cè)網(wǎng)絡(luò)僅能覆蓋數(shù)百平方千米的范圍，而且難以保證監(jiān)測(cè)的實(shí)效性。合成孔徑雷達(dá)干涉測(cè)量(Synthetic Aperture Radar Interferometry，InSAR)及其衍生技術(shù)的出現(xiàn)為解決這一難題提供了強(qiáng)有力的技術(shù)支撐。然而，隨著地面沉降區(qū)域的不斷擴(kuò)大，受到SAR影像幅寬范圍(如Envisat ASAR標(biāo)準(zhǔn)幅寬為100 km×100 km)的限制，單幅SAR影像無(wú)法滿足整個(gè)區(qū)域或者整個(gè)盆地尺度的地表形變監(jiān)測(cè)需求，給大范圍的地質(zhì)災(zāi)害形變監(jiān)測(cè)帶來(lái)不便。另外，盡管時(shí)序InSAR技術(shù)(如永久散射體技術(shù)、小基線集技術(shù)等)已廣泛應(yīng)用于地表形變監(jiān)測(cè)中(朱煜峰等, 2017;董少春等, 2019; 王升等，2019;Liu et al., 2018b, 2019)，卻只聚焦于小范圍、局部地區(qū)的形變分析，同時(shí)進(jìn)行數(shù)千乃至數(shù)萬(wàn)平方千米的形變監(jiān)測(cè)研究尚不多見(李廣宇等, 2018; 張永紅等, 2016; Chaussard et al., 2014; Peng et al., 2020; Yang et al., 2014)。

要實(shí)現(xiàn)高精度的大范圍地表形變監(jiān)測(cè)，采用的SAR影像必須同時(shí)兼顧在方位向能夠連續(xù)獲取長(zhǎng)條帶SAR數(shù)據(jù)和在距離向覆蓋范圍廣(葛大慶, 2013)。然而，隨著數(shù)據(jù)處理范圍的不斷擴(kuò)大，影響形變信號(hào)解算精度的干擾因素(如基線誤差、大氣附加相位、相位解纏誤差等)也隨之增加，如何有效降低這些不利因素的影響是時(shí)序InSAR大范圍地表形變監(jiān)測(cè)面臨的主要問(wèn)題之一。同一形變信號(hào)在不同入射角下的InSAR監(jiān)測(cè)結(jié)果在平面位置上并不能嚴(yán)格的一一對(duì)應(yīng)，而且僅能獲取相對(duì)給定參考點(diǎn)的形變量，導(dǎo)致相鄰軌道的形變解算結(jié)果間存在明顯的參考基準(zhǔn)偏差，如何實(shí)現(xiàn)相鄰軌道InSAR形變監(jiān)測(cè)結(jié)果的坐標(biāo)系與參考基準(zhǔn)的統(tǒng)一是時(shí)序InSAR大范圍地表形變監(jiān)測(cè)亟待解決的又一主要問(wèn)題。

針對(duì)上述問(wèn)題，筆者開展時(shí)序InSAR大范圍地表形變監(jiān)測(cè)的關(guān)鍵技術(shù)分析，并將此應(yīng)用于山西省太原盆地，獲取盆地尺度的高精度地表形變監(jiān)測(cè)結(jié)果。

1 長(zhǎng)條帶SAR數(shù)據(jù)的時(shí)序InSAR監(jiān)測(cè)

長(zhǎng)條帶SAR數(shù)據(jù)的時(shí)序InSAR監(jiān)測(cè)即對(duì)連續(xù)掃描獲取的SAR影像序列進(jìn)行分析，以獲取長(zhǎng)條帶模式下的大范圍地表形變場(chǎng)。

1.1 長(zhǎng)條帶SAR數(shù)據(jù)的分塊處理法

考慮到計(jì)算效率、誤差傳遞等因素對(duì)大范圍形變監(jiān)測(cè)的影響，尤其是軌道誤差、大氣附加相位等長(zhǎng)波長(zhǎng)誤差的影響，筆者采用分塊處理方式。其基本思路是首先將長(zhǎng)條帶SAR影像在方位向與距離向分別進(jìn)行分塊，分割成若干個(gè)具有一定重疊區(qū)域的數(shù)據(jù)子塊，然后在各子塊內(nèi)利用時(shí)序InSAR技術(shù)進(jìn)行參數(shù)估算，最后對(duì)估算出的線性形變速率、形變時(shí)間序列等進(jìn)行拼接處理從而獲取大范圍的地表形變場(chǎng)。

1.2 軌道趨勢(shì)性誤差改正

由于SAR影像獲取時(shí)衛(wèi)星定軌精度可能不夠精確，導(dǎo)致干涉組合的空間基線估值產(chǎn)生偏差，因此在長(zhǎng)條帶時(shí)序InSAR形變監(jiān)測(cè)中，軌道殘余相位成為影響形變監(jiān)測(cè)精度不可忽略的因素之一。當(dāng)監(jiān)測(cè)區(qū)域范圍逐漸擴(kuò)大時(shí)，軌道殘余相位隨之逐漸累積，最終在形變圖中表現(xiàn)為一種趨勢(shì)性的影響。鑒于此，筆者通過(guò)建立多項(xiàng)式模型擬合相位解纏的趨勢(shì)面，以消除其對(duì)形變解算的影響。數(shù)學(xué)模型表示如下：

φorb=a0+a1x+a2y+a3x2+a4y2+a5xy+a6h

(1)

式中，φorb表示軌道誤差相位，(x,y)表示參與多項(xiàng)式模型構(gòu)建的像元行列坐標(biāo)，h表示參與多項(xiàng)式擬合運(yùn)算的像元高程值，a0,a1,a2,a3,a4,a5,a6為模型擬合系數(shù)(即待求參數(shù))。

1.3 連續(xù)刈幅的參考基準(zhǔn)統(tǒng)一

在各數(shù)據(jù)子塊中利用時(shí)序InSAR技術(shù)估算的線性形變速率、形變時(shí)間序列等參數(shù)具有各自獨(dú)立的參考基準(zhǔn)，且各參考點(diǎn)的活動(dòng)性通常未知。由于各數(shù)據(jù)子塊具有相同的軌道編號(hào)且已進(jìn)行精確配準(zhǔn)，因此連續(xù)刈幅參考基準(zhǔn)統(tǒng)一的實(shí)質(zhì)就是補(bǔ)償不同數(shù)據(jù)子塊間的基準(zhǔn)偏差(可看作一個(gè)系統(tǒng)偏差)，通過(guò)對(duì)重疊區(qū)域內(nèi)高相干點(diǎn)目標(biāo)的形變參數(shù)進(jìn)行統(tǒng)計(jì)分析實(shí)現(xiàn)。

以線性形變速率為例，連續(xù)刈幅參考基準(zhǔn)統(tǒng)一的基本思路是首先選定一個(gè)子塊作為參考，提取出所有相干點(diǎn)目標(biāo)的形變速率，對(duì)其余數(shù)據(jù)子塊的形變監(jiān)測(cè)結(jié)果進(jìn)行空間內(nèi)插，生成連續(xù)的地表形變場(chǎng)；然后根據(jù)參考子塊上相干點(diǎn)目標(biāo)的空間位置提取其余數(shù)據(jù)子塊對(duì)應(yīng)的形變監(jiān)測(cè)結(jié)果；最后按照公式(2)計(jì)算重疊區(qū)域內(nèi)同名高相干點(diǎn)線性形變速率的整體偏差，從而完成連續(xù)刈幅線性形變速率參考基準(zhǔn)的統(tǒng)一。

(2)

2 相鄰軌道InSAR監(jiān)測(cè)成果的拼接

2.1 坐標(biāo)系的統(tǒng)一

為實(shí)現(xiàn)相鄰軌道長(zhǎng)條帶時(shí)序InSAR監(jiān)測(cè)成果的有效拼接，首要工作是在空間上統(tǒng)一相鄰軌道InSAR監(jiān)測(cè)成果坐標(biāo)系(空間參考基準(zhǔn))?；舅悸肥菍⒉煌壍老碌腎nSAR觀測(cè)值均轉(zhuǎn)換至地理坐標(biāo)系，然后以某一軌道對(duì)應(yīng)的地理坐標(biāo)系作為參考基準(zhǔn)(主軌道)，將其他地理坐標(biāo)系下(輔軌道)的InSAR觀測(cè)值通過(guò)坐標(biāo)轉(zhuǎn)換使其位于統(tǒng)一的坐標(biāo)系。

2.2 參考基準(zhǔn)的統(tǒng)一

與長(zhǎng)條帶時(shí)序InSAR監(jiān)測(cè)中的參考基準(zhǔn)統(tǒng)一類似，相鄰軌道InSAR監(jiān)測(cè)成果參考基準(zhǔn)統(tǒng)一的核心問(wèn)題同樣在于如何有效提取相鄰軌道重疊區(qū)域內(nèi)高相干點(diǎn)目標(biāo)的InSAR形變參數(shù)。在此基礎(chǔ)上估算相鄰軌道監(jiān)測(cè)結(jié)果間的整體參考基準(zhǔn)偏差。由于不同軌道SAR數(shù)據(jù)對(duì)應(yīng)的入射角各異，即使是地面上的同一地物在不同軌道SAR影像上的空間位置也會(huì)略有差異。因此，根據(jù)主軌道上高相干點(diǎn)目標(biāo)的空間位置，提取輔軌道上對(duì)應(yīng)的高相干點(diǎn)目標(biāo)的形變參數(shù)，利用式(2)估計(jì)參考基準(zhǔn)整體偏差，進(jìn)而統(tǒng)一相鄰軌道InSAR監(jiān)測(cè)成果參考基準(zhǔn)?？紤]到相鄰軌道各高相干點(diǎn)目標(biāo)與形變解算參考點(diǎn)的距離對(duì)誤差傳播的影響，故采用反距離加權(quán)平均法公式(3)估算重疊區(qū)域高相干點(diǎn)目標(biāo)的形變值。

(3)

3 應(yīng)用分析

為檢驗(yàn)上述時(shí)序InSAR方法的可行性，本次選取山西省太原盆地作為實(shí)驗(yàn)區(qū)，進(jìn)行時(shí)序InSAR大范圍地表形變監(jiān)測(cè)試驗(yàn)與分析。

3.1 研究區(qū)概況

太原盆地位于山西省中部，東臨太行山，西傍呂梁山，覆蓋面積約6 195 km2，涉及太原市、晉中市、孝義市、介休市、清徐縣、交城縣等。該盆地是目前我國(guó)地面沉降等地質(zhì)災(zāi)害最為嚴(yán)重的地區(qū)之一。太原盆地作為一個(gè)典型的內(nèi)陸新生代斷陷盆地，其地貌格局完全受控于活動(dòng)構(gòu)造，主要受到交城斷裂帶和范村-洪山斷裂的控制(劉媛媛, 2018)。

3.2 數(shù)據(jù)源及數(shù)據(jù)處理

3.2.1 SAR數(shù)據(jù)

本研究采用的SAR影像主要包括2008年8月至2010年9月的20景Envisat ASAR數(shù)據(jù)和2007年2月至2011年3月的55景ALOS PALSAR-1數(shù)據(jù)，數(shù)據(jù)覆蓋范圍如圖1所示。其中，Envisat ASAR數(shù)據(jù)包含2個(gè)連續(xù)的圖幅，覆蓋范圍約為150 km×100 km；ALOS PALSAR-1數(shù)據(jù)包含3個(gè)相鄰的軌道(Track455, Track456, Track457)，各軌道包含3個(gè)連續(xù)的圖幅，數(shù)據(jù)覆蓋范圍約為210 km×70 km，相鄰軌道的重疊度約為30%，地面距離約為21 km。

圖1 太原盆地的SAR數(shù)據(jù)覆蓋范圍Fig.1 Spatial coverage of the SAR datasets in Taiyuan basin

3.2.2 時(shí)序InSAR處理

為實(shí)現(xiàn)太原盆地高精度的地表形變監(jiān)測(cè)，采用多級(jí)配準(zhǔn)策略實(shí)現(xiàn)主、從影像間的精確配準(zhǔn)與重采樣，并采用小基線模式組合干涉對(duì)，Envisat ASAR和ALOS PALSAR-1數(shù)據(jù)的時(shí)間閾值分別設(shè)為250天和1 000天，空間基線閾值分別設(shè)為500 m和1 500 m。利用歐空局提供的DORIS精密軌道數(shù)據(jù)進(jìn)行基線精化，并利用研究區(qū)30 m分辨率的SRTM DEM數(shù)據(jù)模擬與去除平地和地形相位的影響。

由于太原盆地內(nèi)既存在城鎮(zhèn)等高相干點(diǎn)目標(biāo)密集的區(qū)域，又包含大片的農(nóng)田等低相干區(qū)域，因此本次主要采用干涉點(diǎn)目標(biāo)分析技術(shù)(劉媛媛, 2018; Werner et al., 2003)進(jìn)行時(shí)序InSAR處理，獲取各高相干點(diǎn)目標(biāo)的形變時(shí)間序列和年平均形變速率圖。

3.3 結(jié)果分析

3.3.1 相鄰軌道InSAR監(jiān)測(cè)成果的拼接

為獲取太原盆地完整的地表形變場(chǎng)，需要對(duì)3個(gè)相鄰軌道的ALOS PALSAR-1形變監(jiān)測(cè)結(jié)果進(jìn)行拼接。鑒于不同軌道SAR數(shù)據(jù)獲取時(shí)刻存在差異，因此僅能對(duì)相鄰軌道獲取的年平均形變速率進(jìn)行有效拼接。此外，由于相鄰軌道SAR數(shù)據(jù)入射角略有不同，假定研究區(qū)地表形變僅發(fā)生在垂直方向，利用各軌道的入射角將視線向形變速率投影至垂直向，完成相鄰軌道的成果拼接。在本實(shí)驗(yàn)中，以Track456作為參考軌道，首先將Track455和Track457的形變速率投影轉(zhuǎn)換至該坐標(biāo)系，然后統(tǒng)計(jì)相鄰軌道重疊區(qū)域的同名高相干點(diǎn)目標(biāo)，利用公式(3)進(jìn)行參考基準(zhǔn)偏差校正，得到具有統(tǒng)一參考基準(zhǔn)的太原盆地尺度的年平均形變速率(圖2)。從圖2可以看出，相鄰軌道InSAR形變監(jiān)測(cè)結(jié)果在重疊區(qū)域具有較好的一致性，表明筆者提出的時(shí)序InSAR大范圍形變監(jiān)測(cè)方法具有一定的可行性。

圖2 ALOS PALSAR-1數(shù)據(jù)的年平均形變速率(2007—2011年)Fig.2 Annual deformation rate derived from ALOS PALSAR-1 from 2007 to 2011

3.3.2 形變結(jié)果分析

圖3為Envisat ASAR數(shù)據(jù)獲取的2008年8月至2010年9月的太原盆地年平均形變速率。圖2和圖3均顯示，2007—2011年期間整個(gè)盆地普遍存在不同程度的地面沉降災(zāi)害，沉降區(qū)具有明顯的差異特征，主要分布在太原盆地的重點(diǎn)城市區(qū)域，如太原市中南部地區(qū)、清徐縣、交城縣、祁縣、太谷縣以及介休市等區(qū)域，其年平均沉降速率一般超過(guò)90 mm/a，最大沉降中心位于清徐縣，沉降速率高達(dá)240 mm/a。筆者分析發(fā)現(xiàn)形變區(qū)與太原盆地地下水源地的空間分布具有一定的關(guān)聯(lián)性，沉降區(qū)基本與地下水開采嚴(yán)重的區(qū)域?qū)?yīng)。Liu等(2018a, 2018b)詳細(xì)分析了盆地典型區(qū)域(太原市、清徐-交城區(qū)域)的地面沉降與地下水之間的相關(guān)性。此外，從空間分布來(lái)看，太原盆地的地面沉降受到斷裂帶的控制，斷裂帶兩側(cè)具有明顯的差異沉降特征。

圖3 Envisat ASAR數(shù)據(jù)的年平均形變速率(2008—2010年)Fig.3 Annual deformation rate derived from Envisat ASAR from 2008 to 2010

3.3.3 精度評(píng)定

本次首先對(duì)ALOS PALSAR-1相鄰軌道重疊區(qū)域的形變觀測(cè)值進(jìn)行求差統(tǒng)計(jì)，檢驗(yàn)獨(dú)立觀測(cè)值的精度(內(nèi)符合精度評(píng)定)。圖4為Track455與Track456、Track457與Track456重疊區(qū)內(nèi)線性形變速率差值的統(tǒng)計(jì)直方圖，其基本符合正態(tài)分布，對(duì)應(yīng)的標(biāo)準(zhǔn)差分別為3.7 mm/a和4.2 mm/a，表明本次提出的時(shí)序InSAR大范圍地表形變監(jiān)測(cè)方法是可靠的。另外，筆者還對(duì)ALOS PALSAR-1與Envisat ASAR數(shù)據(jù)重疊區(qū)域的形變速率進(jìn)行了求差統(tǒng)計(jì)(圖5)，對(duì)應(yīng)的差值標(biāo)準(zhǔn)差為5.6 mm/a。經(jīng)分析發(fā)現(xiàn)，差值較大的區(qū)域恰好與形變速率較大的區(qū)域?qū)?yīng)(例如清徐-交城區(qū)域)，可能是受水平位移影響所致。

圖4 ALOS PALSAR-1相鄰軌道形變速率差值統(tǒng)計(jì)直方圖Fig.4 Histograms of deformation rate difference between the results from adjacent track ALOS PALSAR-1

圖5 ALOS PALSAR-1與Envisat ASAR形變速率差值(a)及對(duì)應(yīng)的統(tǒng)計(jì)直方圖(b)Fig.5 Deformation rate difference between the results from ALOS PALSAR-1 and Envisat ASAR (a) and the corresponding histogram (b)

此外，筆者還利用中國(guó)地震局第二監(jiān)測(cè)中心提供的16個(gè)GPS監(jiān)測(cè)點(diǎn)位(圖3所示)觀測(cè)值評(píng)價(jià)InSAR形變監(jiān)測(cè)結(jié)果的可靠性(外符合精度評(píng)定)?？紤]到InSAR監(jiān)測(cè)的是各高相干點(diǎn)目標(biāo)的形變值，而GPS技術(shù)反映的僅是特定監(jiān)測(cè)點(diǎn)的變形，二者在空間位置上可能存在不匹配現(xiàn)象。因此，本次以GPS監(jiān)測(cè)點(diǎn)為中心提取100 m范圍內(nèi)的高相干點(diǎn)目標(biāo)，采用反距離加權(quán)平均方法獲取各GPS監(jiān)測(cè)點(diǎn)對(duì)應(yīng)位置的InSAR形變值。InSAR與GPS間的相關(guān)性分析(圖6)表明，InSAR形變監(jiān)測(cè)結(jié)果與GPS監(jiān)測(cè)數(shù)據(jù)間具有很高的一致性，相關(guān)系數(shù)分別為0.992 8(ALOS PALSAR-1)和0.988 2(Envisat ASAR)，與GPS間差值的均方根誤差(RMSE)分別為2.8 mm/a(ALOS PALSAR-1)和2.5 mm/a(Envisat ASAR)，表明本次監(jiān)測(cè)結(jié)果具有較高的精度。

圖6 InSAR與GPS之間相關(guān)性分析Fig.6 Correlation analysis between InSAR result and GPS observations

4 結(jié)語(yǔ)

筆者針對(duì)時(shí)序InSAR大范圍地表形變監(jiān)測(cè)的關(guān)鍵問(wèn)題開展了有針對(duì)性的研究，包括長(zhǎng)條帶SAR數(shù)據(jù)分塊處理法以及相鄰軌道InSAR監(jiān)測(cè)成果的拼接，并以山西省太原盆地為研究對(duì)象開展了應(yīng)用試驗(yàn)分析，獲取了山西省太原盆地尺度的高精度地表形變監(jiān)測(cè)結(jié)果。結(jié)果顯示，ALOS PALSAR-1和Envisat ASAR數(shù)據(jù)與GPS監(jiān)測(cè)數(shù)據(jù)的相關(guān)系數(shù)分別為0.992 8和0.988 2，與GPS間差值的均方根誤差分別為2.8 mm/a和2.5 mm/a，表明本次提出的時(shí)序InSAR技術(shù)有效地解決了大范圍地表形變監(jiān)測(cè)同步觀測(cè)的難題，實(shí)現(xiàn)了地表形變大范圍、高精度的監(jiān)測(cè)，有效地揭示了區(qū)域性地表形變的空間分布特征，展現(xiàn)了常規(guī)地面監(jiān)測(cè)手段無(wú)法比擬的技術(shù)優(yōu)勢(shì)，對(duì)于宏觀揭示盆地尺度的地表形變的空間變化特征具有重要意義。