基于時間序列InSAR的巴東地表形變分析

婁連惠，劉 強(qiáng)，譚玉敏

(1.中國地質(zhì)環(huán)境監(jiān)測院，北京 100081； 2.北京航空航天大學(xué) 交通科學(xué)與工程學(xué)院，北京 102206)

1 研究背景

長江三峽地區(qū)位于川鄂交接的山地峽谷地區(qū)，山多坡陡，人口分布不均勻，主要集中在長江沿岸河谷地區(qū)。由于過度墾牧，生態(tài)惡化，一旦遇上暴雨天氣或地震災(zāi)害極易發(fā)生滑坡、泥石流或巖崩災(zāi)害[1]。隨著三峽庫區(qū)沿岸經(jīng)濟(jì)日益發(fā)展，地質(zhì)災(zāi)害帶來的損失也變得越來越嚴(yán)重，因此對三峽庫區(qū)沿岸地表形變情況進(jìn)行監(jiān)測和分析對三峽水庫的安全運行和區(qū)域地質(zhì)災(zāi)害的監(jiān)測預(yù)防有著重要的作用。巴東屬于三峽庫區(qū)建成后的遷移建設(shè)城區(qū)，地質(zhì)環(huán)境條件非常復(fù)雜，屬于三峽庫區(qū)地質(zhì)災(zāi)害易發(fā)區(qū)[2]。

20世紀(jì)50年代末出現(xiàn)的干涉合成孔徑雷達(dá)(Interferometric Synthetic Aperture Radar，InSAR)技術(shù)隨著相關(guān)研究的不斷發(fā)展，已被證明在地面沉降測量方面可以達(dá)到與GPS相當(dāng)?shù)木?，同時兼具覆蓋范圍廣、可實現(xiàn)連續(xù)觀測以及成本低廉等優(yōu)勢，主要包括差分InSAR(Differential InSAR，D-InSAR)和近年發(fā)展起來的時間序列InSAR技術(shù)[3]。時間序列InSAR技術(shù)主要包括永久散射體(Permanent Scatterers,PS)和小基線子集(Small Baseline Subsets，SBAS)等。D-InSAR技術(shù)因受大氣效應(yīng)和失相關(guān)影響顯著，在干涉和差分干涉測量精度方面有一定的局限性[4]；PS技術(shù)可以消除大部分大氣延遲，可以監(jiān)測在實際地面上變形較小并且變化較平穩(wěn)的區(qū)域，但要想得到可靠的 PS 點和結(jié)果，至少需要 20個時相的獲取時間規(guī)律、連續(xù)的數(shù)據(jù)，適用于干涉條件和輻射比較穩(wěn)定的區(qū)域[5]；SBAS技術(shù)利用比較少的數(shù)據(jù)可組合出較多的干涉對，且獲取的主影像干涉對的時間和空間基線小，可以很好地削弱或消除時空失相關(guān)對結(jié)果的影響[6]。

時間序列InSAR技術(shù)已經(jīng)廣泛應(yīng)用于多個領(lǐng)域，如滑坡檢測[7-9]、公路形變檢測[10]、地震形變監(jiān)測[11-12]以及城市形變監(jiān)測[13]等。隨著InSAR技術(shù)發(fā)展及以Sentinel-1A為代表的免費可用InSAR數(shù)據(jù)源的快速增加，進(jìn)行區(qū)域長時間序列地表形變監(jiān)測成為了可能。本文選取地質(zhì)條件復(fù)雜的巴東地區(qū)長江沿岸為研究區(qū)域，利用2017—2019年間Sentinel-1A數(shù)據(jù)分析巴東地區(qū)地表形變情況。

2 研究區(qū)域

巴東處于長江三峽巫峽與西陵峽過渡地帶，是長江三峽沿岸地質(zhì)條件最為復(fù)雜的地區(qū)之一，區(qū)內(nèi)地質(zhì)構(gòu)造復(fù)雜,褶皺、斷層密集，整個巴東區(qū)域地勢陡峭，相對高程最高可達(dá)600 m。如此復(fù)雜的地形條件導(dǎo)致巴東地區(qū)地質(zhì)災(zāi)害(滑坡、泥石流等)頻發(fā)。研究區(qū)域如圖1所示(底圖為谷歌地圖)。

圖1 研究區(qū)范圍

3 數(shù)據(jù)和方法

3.1 數(shù)據(jù)源

本文共選取了53景Sentinel-1A數(shù)據(jù)，所有數(shù)據(jù)都是從升軌模式獲得的，極化方式為VV，時間跨度為2017年4月17日—2019年3月26日。C波段雖然波長較短，對樹葉穿透性差，但是對于城市等高相干性區(qū)域的形變有很高的敏感性。將2018年3月7日一景數(shù)據(jù)作為超級主影像，設(shè)置時間基線閾值為120 d，空間基線閾值為垂直基線的1.2%，共組成了231個干涉對。此外，本文使用了30 m分辨率的SRTM DEM用來去除地形相位信息。

本文也使用了與雷達(dá)數(shù)據(jù)獲取時間相對應(yīng)的巴東地區(qū)長江水位信息，數(shù)據(jù)來源于中國長江三峽集團(tuán)有限公司官網(wǎng)，所有水位信息獲取的時間均為當(dāng)天的14時。

3.2 處理方法

本文基于SBAS技術(shù)提取巴東地區(qū)地表形變信息。SBAS技術(shù)通過采用多對短時間和空間基線下形成的高相干性的干涉像對進(jìn)行解算，得到實際地表形變相位后，通過提取高相干性點，采用奇異值分解求得地表形變速率的最小二乘解，進(jìn)而得到整個時間序列的形變速率。具體處理步驟如下：

(1)生成連接圖。根據(jù)時間和空間基線情況進(jìn)行干涉像對配對，配對的像對將會進(jìn)行干涉處理。

(2)干涉工作流處理。該步驟對所有步驟(1)中生成的干涉像對進(jìn)行干涉處理，包括生成相干性系數(shù)圖、去平、濾波和相位解纏。

(3)軌道精煉和重去平。該步驟利用手動選取的地面相干性較高且相位信息較為連續(xù)的控制點，結(jié)合從SRTM DEM中提取的高程信息估算和去除殘余相位以及解纏后依然存在的相位坡道。

(4)形變速率反演。該步驟是SBAS技術(shù)的核心步驟之一，主要目的是估算形變速率和殘余相位并進(jìn)行大氣濾波，通過時間序列形變情況估算和去除大氣相位的影響，進(jìn)而得到時間序列形變圖。

(5)地理編碼。對所獲得的雷達(dá)坐標(biāo)系下的速率圖進(jìn)行地理編碼。具體流程如圖2所示。

圖2 SBAS處理流程

SBAS流程中，相位解纏是SBAS技術(shù)中最重要的步驟之一，解纏結(jié)果直接決定了形變結(jié)果的精度。目前存在多種解纏算法，如路徑跟蹤解纏算法[14-15]、最小二乘法[16]、基于網(wǎng)絡(luò)規(guī)劃的最小費用流法[17]等。最小費用流法通過先求解纏繞相位與解纏相位的差的最小值，進(jìn)而將求解轉(zhuǎn)換為最小費用流的問題，該算法在計算過程中會根據(jù)流的大小和方向?qū)Ω鱾€相位矩陣進(jìn)行積分，得到高質(zhì)量的解纏結(jié)果，再向低質(zhì)量區(qū)域進(jìn)行積分，得到所有的結(jié)果。這樣能夠很有效地避免相位誤差從低相干性區(qū)域傳到高相干性區(qū)域，對于巴東地區(qū)，地表情況復(fù)雜，城區(qū)以外的范圍植被覆蓋茂密，形成了大范圍的低相干性區(qū)域，采用最小費用流算法進(jìn)行解纏可以得到高精確度的城區(qū)形變結(jié)果。

根據(jù)上述步驟，解算得到了巴東地區(qū)在2017年4月17日—2019年3月26日之間的年平均形變速率，如圖3所示。

圖3 巴東地區(qū)年平均形變速率

4 結(jié)果分析

由圖3可以看出，巴東地區(qū)整體形變速率較小，主要的城市區(qū)域形變速率均維持在-5～5 mm/a之間。其中形變速率正值代表地表形變朝向衛(wèi)星方向，負(fù)值代表地表形變遠(yuǎn)離衛(wèi)星方向。由于水面的相干性較低，圖3中對于水體部分沒有相應(yīng)的形變結(jié)果，也從側(cè)面反映出解纏結(jié)果的準(zhǔn)確性。其次由于巴東地區(qū)地表情況復(fù)雜，有些區(qū)域不僅植被覆蓋茂密也存在較大的坡度，這些都會導(dǎo)致失相關(guān)現(xiàn)象的產(chǎn)生[18]，表現(xiàn)在圖3中最右側(cè)山體解纏結(jié)果的缺失以及部分森林區(qū)域像素點缺失，但城區(qū)等整體形變結(jié)果良好。

三峽大壩建成以后，多位學(xué)者從不同角度做了水庫和地表形變間的相關(guān)性研究：王林松等[19]基于DEM數(shù)據(jù)建立了三峽庫區(qū)蓄水負(fù)荷和地表形變的響應(yīng)模型；王偉等[20]通過建立水體重力與地表形變的關(guān)系預(yù)測了三峽沿岸的地表形變情況。但是三峽庫區(qū)地質(zhì)情況極其復(fù)雜，本文在InSAR地表形變分析基礎(chǔ)上，對研究區(qū)按與岸線距離遠(yuǎn)近進(jìn)行了分區(qū)分析，每隔185 m設(shè)置緩沖區(qū)，分區(qū)劃分以及點群獲取情況如圖4所示。

圖4 點群分割示意圖

為了保證分析的整體精度，本文以0.80的高相干性系數(shù)閾值對獲取的形變圖進(jìn)行點群提取，提取的點群如圖4所示，按照距離江邊遠(yuǎn)近順序?qū)Ⅻc群分割成了點群5—點群1。由圖4可以看出，提取的點群分布范圍覆蓋了研究區(qū)的大部分區(qū)域，且整體獲得的點群數(shù)量充足，因此點群情況足以反映研究區(qū)域內(nèi)的整體形變情況。在一些樹木較多或者存在較大坡度的區(qū)域，數(shù)據(jù)相干性較低，沒有相應(yīng)的點群出現(xiàn)，這也保證了用以分析的數(shù)據(jù)的精確性。具體點群數(shù)量如表1所示。

表1 不同緩沖區(qū)內(nèi)點群數(shù)量

為了統(tǒng)計不同緩沖區(qū)內(nèi)的整體形變趨勢，將緩沖區(qū)1—緩沖區(qū)5內(nèi)的點分別進(jìn)行統(tǒng)計，其中每一個點都包含2017年4月17日—2019年3月26日之間的時間序列形變信息，時間序列形變情況的時間分辨率由相同時間內(nèi)所輸入的Sentinel-1A數(shù)據(jù)的數(shù)據(jù)量決定。通過求得每一個時期衛(wèi)星視線向的形變的平均值反映該緩沖內(nèi)的整體形變趨勢。按照不同緩沖區(qū)內(nèi)形變趨勢得出的時間序列形變折線圖如圖5所示。

圖5 不同緩沖區(qū)內(nèi)時間序列形變折線圖

從圖5可知，各個緩沖區(qū)的形變均在-12～4 mm之間，形變較小，形變的整體趨勢比較明顯，隨著點群1—點群4與江邊的距離變遠(yuǎn)，緩沖區(qū)1—緩沖區(qū)4內(nèi)的形變趨勢呈現(xiàn)出逐步減小的態(tài)勢，并且緩沖區(qū)4和緩沖區(qū)5內(nèi)的形變趨勢逐步失去規(guī)律性。

為了更清晰地表達(dá)距離江邊較近的緩沖區(qū)內(nèi)形變趨勢和水位變化在時間上的相互關(guān)系，將緩沖區(qū)1—緩沖區(qū)3內(nèi)的形變與水位做了圖6所示的分析。圖6中最上方的折線表示水位高度線，隨著三峽大壩的蓄放水，水位規(guī)律呈現(xiàn)出周期性，由圖6可以看出，水位變化周期和緩沖區(qū)1—緩沖區(qū)3內(nèi)的形變趨勢在時間上呈現(xiàn)出一定的相關(guān)性。

圖6 緩沖區(qū)1—緩沖區(qū)3內(nèi)的形變與水位復(fù)合圖

5 結(jié) 語

本文利用免費數(shù)據(jù)Sentinel-1A完成了巴東地區(qū)地表形變分析，結(jié)果顯示巴東地區(qū)整體形變速率較小，趨于穩(wěn)定狀態(tài)。由于該區(qū)域所得Sentinel-1A數(shù)據(jù)均為升軌數(shù)據(jù)，無法結(jié)合同期降軌數(shù)據(jù)進(jìn)行分析，因此文中所得到的形變速率以及形變方向均為衛(wèi)星視線向。

分析結(jié)果顯示在距離岸邊一定范圍內(nèi)的區(qū)域，其地形變速率與水位變化規(guī)律有一定的一致性，但考慮到降水等其他重要因素對地表形變的影響在時間周期上與水位變化重合，而本文沒有將同期降雨數(shù)據(jù)一并分析，所以本文結(jié)論中按江邊距離分區(qū)進(jìn)行的水位變化關(guān)聯(lián)性分析僅是從區(qū)域統(tǒng)計角度提供相對定量的地形變分析結(jié)果，可能不夠充分。