利用SBAS-InSAR技術(shù)監(jiān)測菏澤市地面沉降

張亞鳳 劉國林 牛 沖 陳 洋 周一鳴

1 山東科技大學(xué)測繪與空間信息學(xué)院，青島市前灣港路579號，266590

2 山東省地質(zhì)測繪院，濟(jì)南市二環(huán)東路11101號，250014

研究地面沉降空間分布形態(tài)和產(chǎn)生因素對社會穩(wěn)定及災(zāi)害防控具有重要的實用價值和科學(xué)意義。合成孔徑雷達(dá)差分干涉測量(D-InSAR)技術(shù)是一種新型的空間對地觀測技術(shù)，具有全天候、全天時、分辨率高、監(jiān)測范圍廣等特點[1-2]。但這種技術(shù)易受時空失相干、大氣延遲和植被覆蓋等影響，在地形復(fù)雜地區(qū)監(jiān)測能力受限[3]。短基線集合成孔徑雷達(dá)干涉測量(SBAS-InSAR)技術(shù)能有效解決時空失相干的影響，實現(xiàn)長時間連續(xù)監(jiān)測[4]。目前已有大量學(xué)者采用SBAS-InSAR技術(shù)監(jiān)測地面沉降，結(jié)果表明，該方法具有較高的準(zhǔn)確性和可靠性[5-7]。

本文利用SBAS-InSAR技術(shù)對菏澤市2017-05-20～2021-05-23的Sentinel-1A SAR影像進(jìn)行處理，獲取該時段的年平均沉降速率和累積沉降量，并結(jié)合實際煤礦開采情況對研究區(qū)域累積沉降的時空演化和時序沉降場進(jìn)行精細(xì)化分析。

1 研究區(qū)域概況與數(shù)據(jù)來源

1.1 研究區(qū)域概況

菏澤市地處山東省西南部，共有2個市轄區(qū)、7個縣、2個開發(fā)區(qū)，面積約12 239 km2，地理范圍為34°39′～35°52′N、114°45′～116°25′E。圖1為研究區(qū)域的具體位置。

1.2 數(shù)據(jù)來源

Sentinel-1A SAR影像數(shù)據(jù)由歐空局提供，搭載C波段的合成孔徑雷達(dá)，對地觀測模式為干涉寬幅模式(IW)，幅寬為250 km，空間分辨率為5 m×20 m，極化方式為VV，重訪周期為12 d，具有高重訪頻率、高覆蓋能力以及極好的時效性和可靠性。本文選取覆蓋菏澤市2017-05-20～2021-05-23期間65景Sentinel-1A SAR影像，同時引入SRTM3 DEM數(shù)據(jù)(分辨率為90 m)進(jìn)行計算，以消除地形相位的影響。

2 數(shù)據(jù)處理方法

SBAS-InSAR技術(shù)是一種用于研究低分辨率、大尺度形變的時間序列分析方法。該方法首先通過設(shè)置時空基線閾值將影像自由組合生成一系列干涉圖，利用最小二乘法求解地表形變信息；然后通過矩陣的奇異值分解(SVD)法進(jìn)行時間序列分析，求得最小范數(shù)最小二乘法解，去除殘余地形，解決時間不連續(xù)的問題；最后將多個短基線集數(shù)據(jù)聯(lián)合起來進(jìn)行反演，得到觀測時段內(nèi)的地表形變信息和平均形變速率。具體計算方法見文獻(xiàn)[8]，本文不再贅述。

根據(jù)短基線原則，SBAS-InSAR技術(shù)將收集到的影像數(shù)據(jù)組合成若干個短基線數(shù)據(jù)集，能有效提高時空的相關(guān)性，提取到更可靠的時序沉降結(jié)果。其主要數(shù)據(jù)處理步驟如下：

1)將原始數(shù)據(jù)轉(zhuǎn)換成SLC格式數(shù)據(jù)，并用實驗區(qū)矢量范圍裁剪影像，選取2017-07-19為超級主影像，其他影像依次與主影像配準(zhǔn)。

2)設(shè)置時空基線閾值和多視比例值(距離向為4，方位向為1)，根據(jù)短基線原則生成143對干涉像對。利用外部DEM模擬地形相位去除干涉圖中的地形相位，得到差分干涉圖。

3)選用Goldstein方法進(jìn)行濾波，采用基于Delaunay三角網(wǎng)的最小費用流法MCF對濾波增強(qiáng)后的差分干涉圖進(jìn)行相位解纏去除噪聲，減少由軌道、大氣和植被等因素引起的誤差。

4)在地形平坦、沒有相位躍變的區(qū)域選取多個控制點進(jìn)行重去平，消除相位偏移量。采用SVD方法估算平均形變速率和殘余地形相位。通過大氣濾波估算和去除大氣延遲相位，從而得到最終地表形變結(jié)果。將生成的結(jié)果編碼至WGS-84坐標(biāo)系下，得到LOS上的平均形變速率和累積形變量。

數(shù)據(jù)處理流程如圖2所示。

3 SBAS-InSAR沉降結(jié)果與時序分析

3.1 SBAS-InSAR沉降結(jié)果

采用SBAS-InSAR技術(shù)處理65景Sentinel-1A數(shù)據(jù)，獲取研究區(qū)域2017-05-20～2021-05-23期間的地面沉降信息(圖3)。由圖可知：

1)在監(jiān)測時段內(nèi)，研究區(qū)域存在一處明顯的沉降盆地，沉降中心位于鄆城縣及周邊區(qū)域，最大沉降速率為-311 mm/a(本文中負(fù)號僅表示方向)，最大累積沉降量達(dá)-1 269 mm。結(jié)合相關(guān)資料分析可知，鄆城縣周圍煤礦眾多，煤礦工作面的持續(xù)開采導(dǎo)致地面沉降加劇，出現(xiàn)沉降盆地。

2)沉降盆地周圍的沉降速率為-40～-20 mm/a，沉降量為-150～-80 mm；距沉降中心較遠(yuǎn)的區(qū)域，沉降速率為-20～0 mm/a，沉降量為-80～0 mm。除鄆城縣這一主要沉降盆地外，巨野縣西北部也有較為嚴(yán)重的沉降，最大年平均沉降速率為-60 mm/a，最大累積沉降量為-250 mm；鄄城縣、牡丹區(qū)和曹縣部分地區(qū)有輕微沉降，年平均沉降速率約為-20 mm/a，累積沉降量為-80～0 mm。

3.2 累積沉降分析

為驗證SBAS-InSAR監(jiān)測結(jié)果的可靠性，選取煤礦分布集中、沉降最嚴(yán)重的鄆城區(qū)域進(jìn)行詳細(xì)分析。以2017-05-20的影像為基準(zhǔn)獲取鄆城沉降的時空分布累積沉降圖，時間間隔為0.5 a(圖4)。其中，A、B、C、D為鄆城區(qū)域4個主要煤礦區(qū)，監(jiān)測時段內(nèi)主要開采工作面如圖5所示(D區(qū)因數(shù)據(jù)資料不足，未給出工作面)。

根據(jù)煤礦開采進(jìn)度，將監(jiān)測過程分為4個階段進(jìn)行分析。

1)階段1：時間為2017-05-20～2018-05-15。由圖4(a)、4(b)可見，煤礦A和煤礦C交界處存在明顯的沉降盆地，沉降中心最大累積沉降量達(dá)-150 mm，周邊沉降量普遍為-100～-20 mm。該時段內(nèi)位于王樓村東南處煤礦A的工作面1開始回采，截至2018-05該工作面已回采1 530 m。SBAS-InSAR監(jiān)測到該時段內(nèi)此工作面附近沉降最為嚴(yán)重，最大累積沉降量達(dá)-300 mm。鄆城其他區(qū)域也存在沉降，主要分散在各個煤礦的中心及周邊城鎮(zhèn)。位于陳河口村附近煤礦B的工作面1和2于2018年初開始開采，監(jiān)測到附近最大累積沉降量為-150 mm。煤礦C前期已完成多個工作面的開采，部分地面沉降主要位于車樓村附近[9]。2018年初煤礦C的工作面1和2開始開采，導(dǎo)致該地區(qū)煤礦沉降最為嚴(yán)重，最大累積沉降量達(dá)-600 mm。煤礦D處前期開采導(dǎo)致地面緩慢沉降，累積沉降量和沉降范圍較小，沉降量普遍為-60～0 mm。

2)階段2：時間為2018-05-15～2019-05-10。由圖4(c)、4(d)可見，隨著時間推移，沉降范圍逐漸增大，向北延伸至王樓村附近，向東延伸至煤礦B西側(cè)，向西延伸至后黃崗村附近，向南延伸至煤礦C車樓村附近。沉降中心最大累積沉降量為-300 mm，周邊累積沉降量集中在-250～-60 mm。2019-05煤礦A的工作面1開采結(jié)束，回采推進(jìn)至2 730 m，最大累積沉降量為-600 mm。煤礦B的工作面1和2以及煤礦C的工作面1和2仍在開采中，且沉降范圍逐步擴(kuò)大，沉降增加，最大累積沉降量均達(dá)-800 mm。2018年底煤礦C的工作面2結(jié)束開采，2019年初工作面3接替開采。煤礦D的最大累積沉降量為-150 mm。

3)階段3：時間為2019-05-10～2020-05-04。由圖4(e)、4(f)可知，沉降范圍持續(xù)擴(kuò)張，南北向最為明顯，王樓村和車樓村處沉降連接。沉降中心最大累積沉降量達(dá)-400 mm，周邊沉降量基本為-300～80 mm。煤礦A的工作面2開采時段為2019-06～2020-04，共開采1 010 m，其附近最大累積沉降量為-800 mm。由于煤礦B的工作面1、2、3和煤礦C的工作面3、4、5持續(xù)開采(其中,2019年底煤礦C的工作面3結(jié)束開采，2020年初工作面4、5開始開采)，造成沉降范圍繼續(xù)擴(kuò)大，部分地區(qū)最大累積沉降量超過-800 mm。煤礦D的沉降量同樣持續(xù)增加，最大累積沉降量為-600 mm。

4)階段4：時間為2020-05-04～2021-05-23。由圖4(g)、4(h)可知，沉降范圍往南北方向擴(kuò)張更加明顯，尤其在郭屯鎮(zhèn)和王樓村附近。沉降中心最大累積沉降量達(dá)-600 mm，周邊累積沉降量基本為-400～-100 mm。煤礦A的工作面3開采時段為2020-05～2021-02，共計開采650 m。之后煤礦B的工作面1、3、4以及煤礦C的工作面4、5、6接替開采。煤礦A、B、C部分區(qū)域最大累積沉降量達(dá)-1 269 mm，煤礦D最大累積沉降量為-800 mm。

綜上所述，隨著時間的推移，地面沉降愈加明顯，無論是沉降量還是沉降范圍都呈現(xiàn)逐漸增大的趨勢。沉降主要集中在各個煤礦開采的工作面附近，最大累積沉降量可達(dá)-1 269 mm。說明鄆城區(qū)域的沉降主要由煤礦開采引起，隨著工作面的持續(xù)開采，沉降量逐漸增大。

為更深入地了解沉降區(qū)域的時序變化情況，選取不同煤礦工作面附近沉降較為嚴(yán)重的6個特征點(圖4(h))，以2017-05-20的影像為基準(zhǔn)，分析其時序累積沉降量(圖6)。由圖可見，隨著時間的推移，P1～P6點累積沉降量逐漸增大，沉降趨勢近似呈線性。其中，P1～P4位于沉降中心，監(jiān)測初期沉降較小，2018年初隨著煤礦開采量的增加，累積沉降量也隨之增大，沉降曲線相對于其他特征點波動較多。截至2021-05-23，P1～P4的累積沉降量分別達(dá)到-939 mm、-453 mm、-820 mm和-865 mm。P5、P6位于沉降邊緣地區(qū)，沉降曲線變化平緩，且沉降量較小，累積沉降量分別為-412 mm和-326 mm。綜上可知，P1點的累積沉降量最大，P6點的累積沉降量最小，符合實際地質(zhì)情況。

3.3 精度驗證

收集鄆城水準(zhǔn)測量數(shù)據(jù)，驗證SBAS-InSAR監(jiān)測沉降的精度。水準(zhǔn)觀測時段為2017-05-13～2018-06-23，水準(zhǔn)點布設(shè)情況如圖7所示。選取2017-05-20～2018-06-08的SBAS-InSAR觀測結(jié)果與水準(zhǔn)測量結(jié)果進(jìn)行對比，并利用線性插值法獲取與SAR影像日期一致的水準(zhǔn)測量結(jié)果，解決水準(zhǔn)與SAR影像數(shù)據(jù)獲取時間不一致的問題。

考慮到SBAS-InSAR的結(jié)果不連續(xù)，分別在鄆城區(qū)域內(nèi)選取沉降量較小的S1～S4點及接近沉降中心處沉降較嚴(yán)重的S5、S6點，對比分析不同沉降程度的SBAS-InSAR和水準(zhǔn)測量的累積沉降結(jié)果(圖8)。

由圖8(a)～(d)可知，水準(zhǔn)點S1～S4位于沉降盆地邊緣地帶，累積沉降量較小。由于2018-04～06期間植被生長茂密，導(dǎo)致SAR影像干涉質(zhì)量差，影響沉降監(jiān)測的精度。截至2018-06-08，4個水準(zhǔn)點的水準(zhǔn)測量結(jié)果分別為-35 mm、-273 mm、-200 mm和-351 mm，SBAS-InSAR監(jiān)測結(jié)果分別為-45 mm、-113 mm、-55 mm和-82 mm，二者累積沉降量絕對差值分別為10 mm、160 mm、145 mm和269 mm?？梢钥闯?，SBAS-InSAR監(jiān)測到的沉降趨勢和沉降量與水準(zhǔn)測量結(jié)果基本相符。

由圖8(e)、8(f)可知，水準(zhǔn)點S5、S6接近沉降盆地中心區(qū)域，累積沉降量較大。2017-05-20～2018-02-08期間S5的SBAS-InSAR監(jiān)測結(jié)果與水準(zhǔn)測量結(jié)果相差不大，且變化趨勢基本一致；2017-05-20～2017-09-17期間S6的SBAS-InSAR監(jiān)測結(jié)果與水準(zhǔn)測量的沉降曲線保持著較好的一致性。隨著時間的推移，水準(zhǔn)測量結(jié)果大幅度增加，而SBAS-InSAR監(jiān)測結(jié)果增加緩慢，截至2018-06-08，S5、S6的水準(zhǔn)測量結(jié)果分別為-900 mm和-2 325 mm，SBAS-InSAR監(jiān)測結(jié)果分別為-122 mm和-185 mm。由于受到植被、大氣等誤差因素的影響，二者結(jié)果相差較大，絕對差值分別為778 mm和2 140 mm。

根據(jù)6個點的SBAS-InSAR監(jiān)測結(jié)果與水準(zhǔn)測量結(jié)果可知，SBAS-InSAR可以監(jiān)測到與水準(zhǔn)測量相符的沉降趨勢，但在沉降量上與水準(zhǔn)測量結(jié)果存在一定的差異。在沉降量較小的區(qū)域，SBAS-InSAR與水準(zhǔn)測量結(jié)果較為吻合；在沉降量較大的區(qū)域，由于InSAR技術(shù)本身的局限性，當(dāng)單個像元內(nèi)的形變超過影像的半個波長時，就無法探測出正確形變，因此InSAR在沉降中心的監(jiān)測能力不足。由于SBAS-InSAR(面單元)和水準(zhǔn)監(jiān)測(點單元)的數(shù)據(jù)類型不同，且監(jiān)測時段不完全一致，SBAS-InSAR監(jiān)測到的沉降量與水準(zhǔn)實測數(shù)據(jù)相差較大，說明SBAS-InSAR技術(shù)在沉降較為嚴(yán)重區(qū)域的精度還有待提升。

4 降雨量對SBAS-InSAR監(jiān)測精度的影響

菏澤地屬溫帶季風(fēng)型大陸性氣候，四季分明，春秋多風(fēng)、夏熱多雨，多年平均降雨量約為662.7 mm，年內(nèi)降水分布不均，主要集中在6～8月。為探究降雨量對SBAS-InSAR監(jiān)測精度的影響，收集菏澤市2017-05～2021-05的月均降雨量，并計算相應(yīng)時段SAR干涉對的月均相干系數(shù)，分析降雨量和相干系數(shù)的關(guān)系，結(jié)果如圖9所示。由圖9可見，菏澤市全年降雨主要集中在6～8月，其中，8月月均降雨量最大，為199 mm；1～5月和9～12月降雨量較少，其中，12月月均降雨量最小(9 mm)。在降雨量較大的6～8月，充沛的雨量和茂盛的植被導(dǎo)致SAR影像質(zhì)量較差，相干性較低，基本為0.37～0.40，其中，8月月均相干系數(shù)最小(0.37)；在降雨量較小的1～5月和9～12月，植被稀疏，SAR影像質(zhì)量較好，相干性較高，基本為0.40～0.63，其中1月月均相干系數(shù)最大(0.63)。由此可知，SAR干涉對的相干性與降雨量有直接關(guān)系，二者呈負(fù)相關(guān)，即降雨量越小，SAR干涉對的相干性越高，監(jiān)測精度就越高。

5 結(jié) 語

1)在監(jiān)測時段內(nèi)，鄆城地區(qū)沉降較為嚴(yán)重，主要是由煤礦開采造成的，最大年平均沉降速率為-311 mm/a，最大累積沉降量達(dá)-1 269 mm；巨野縣西北部地區(qū)也有較為嚴(yán)重的沉降，最大年平均沉降速率為-60 mm/a，最大累積沉降量為-250 mm；鄄城縣、牡丹區(qū)和曹縣沉降量較小，年平均沉降速率基本為-20 mm/a，累積沉降量基本為-80～0 mm。

2)煤礦分布集中、沉降最嚴(yán)重的鄆城地區(qū)累積沉降量最大的區(qū)域主要在王樓村、車樓村和陳河口村附近，與實際地質(zhì)情況相符，最大累積沉降量可達(dá)-1 269 mm。

3)SBAS-InSAR技術(shù)可以準(zhǔn)確探測到沉降位置，反映大范圍的沉降情況和沉降變化趨勢。在沉降量較小的邊緣地區(qū)，其監(jiān)測結(jié)果與水準(zhǔn)測量結(jié)果基本相符，監(jiān)測精度較高；在靠近沉降中心區(qū)域，受到SAR影像特征和技術(shù)自身等限制，SBAS-InSAR監(jiān)測結(jié)果與水準(zhǔn)測量結(jié)果存在較大差異，監(jiān)測精度仍有待提升。

4)根據(jù)月均降雨量和SAR影像的月均相干系數(shù)可知，SAR干涉對的相干性與降雨量有直接關(guān)系，二者呈負(fù)相關(guān)，即降雨量越小，SAR干涉對的相干性越高，監(jiān)測精度就越高。