黃德華

(福建省地質(zhì)測(cè)繪院,福建 福州 350011)

開采地下礦產(chǎn)資源容易破壞巖層結(jié)構(gòu),可能引發(fā)地面沉降等一系列地質(zhì)問題,從而對(duì)沉陷區(qū)的地表建筑、公共設(shè)施、耕地、水環(huán)境等人類生存環(huán)境、經(jīng)濟(jì)財(cái)產(chǎn)安全造成極大的影響[1-4]。為了區(qū)域的可持續(xù)發(fā)展,科學(xué)合理地進(jìn)行開發(fā)、生態(tài)治理,對(duì)礦區(qū)長(zhǎng)時(shí)間序列的動(dòng)態(tài)監(jiān)測(cè)十分必要。

由于礦區(qū)形變周期性長(zhǎng)、范圍廣,而常規(guī)的監(jiān)測(cè)手段具有監(jiān)測(cè)點(diǎn)離散、周期長(zhǎng)、成本高的特點(diǎn),難以進(jìn)行時(shí)空連續(xù)性監(jiān)測(cè)。合成孔徑雷達(dá)差分干涉測(cè)量技術(shù)是通過短周期內(nèi)的影像對(duì)干涉相位的差分計(jì)算,可大范圍觀測(cè)地表形變,其監(jiān)測(cè)成本低,精度可達(dá)厘米至毫米級(jí)[5]。2000年CARNEC C等[6]首次將D-InSAR技術(shù)用于煤礦區(qū)地表形變場(chǎng)監(jiān)測(cè);2005年吳立新等[7]對(duì)中國東部典型工礦區(qū)進(jìn)行了礦區(qū)地表沉陷D-InSAR監(jiān)測(cè)研究。在D-InSAR技術(shù)的基礎(chǔ)上,2001年FERRETTI A等[8]提出了PS-InSAR技術(shù);2002年BERARDINO P等[9]隨后提出了SBAS-InSAR技術(shù),此后更是發(fā)展了眾多時(shí)序監(jiān)測(cè)算法。近年來,多國學(xué)者利用InSAR技術(shù)對(duì)礦區(qū)地面開展了形變監(jiān)測(cè),例如:2015年GRZOVIC M等[10]使用PS-InSAR與SBAS-InSAR方法對(duì)美國伊利諾伊州斯普林菲爾德地下煤礦采區(qū)進(jìn)行了地表形變監(jiān)測(cè),監(jiān)測(cè)結(jié)果顯示地下礦井坍塌導(dǎo)致了該區(qū)域的地表形變;2023年虎小強(qiáng)等[11]使用SBAS-InSAR方法對(duì)新疆阿希礦區(qū)地面沉降時(shí)空變化特征進(jìn)行監(jiān)測(cè),分析了沉降主要驅(qū)動(dòng)因素,等等?，F(xiàn)有研究表明,SBAS-InSAR與PS-InSAR是目前比較有代表性的時(shí)序InSAR算法,可有效監(jiān)測(cè)礦區(qū)地面沉陷。

本文闡述了InSAR技術(shù)原理,運(yùn)用D-InSAR、時(shí)序InSAR技術(shù)對(duì)云岡礦區(qū)開展地面沉降監(jiān)測(cè),其中在2019年1月—2020年12月時(shí)段采用24景哨兵-1A衛(wèi)星影像數(shù)據(jù)進(jìn)行SBAS-InSAR處理獲取形變信息,在2020年1月—2020年9月時(shí)段采用25景Sentinel影像數(shù)據(jù)進(jìn)行PS-InSAR處理獲取形變信息,并選取2020年1月—2020年9月兩種技術(shù)方法不同形變信息進(jìn)行對(duì)比分析,驗(yàn)證兩種時(shí)序InSAR技術(shù)在礦區(qū)地表沉降監(jiān)測(cè)應(yīng)用中的可靠性。

1 研究區(qū)及數(shù)據(jù)概況

1.1 研究區(qū)概況

研究區(qū)位于山西省北部、大同市西部云岡區(qū)內(nèi),面積約137.16 km2,西北部為煤礦開采區(qū),79%以上為地下開采區(qū),其余為露天開采區(qū)。東南部為大同市城市區(qū)域(圖1)。

(a).大同市范圍;(b).研究區(qū)范圍圖1 研究區(qū)概況圖Fig. 1 The overview of the study area

1.2 實(shí)驗(yàn)數(shù)據(jù)

本文采用2019年1月—2020年12月的哨兵-1A衛(wèi)星數(shù)據(jù)及SRTM 30 m分辨率數(shù)字高程數(shù)據(jù)作為數(shù)據(jù)支撐,其中SBAS-InSAR和PS-InSAR采用數(shù)據(jù)情況詳見表1。

表1 哨兵-1A衛(wèi)星影像主要信息表Table 1 Main information of Sentinel-1A satellite image

2 數(shù)據(jù)處理原理及流程

InSAR技術(shù)中相位信息由平地相位、軌道誤差、地形相位、形變相位、大氣相位以及噪聲相位共6部分相位信息組成[12],公式如下:

φdiff=φtopo+φdef+φorb+φatm+φnoise,

式中:φtopo為參考數(shù)字高程模型(Digital Elevationg Model, DEM)殘余地形相位;φdef為形變相位,由線性形變及非線性形變組成;φorb為去平引入的軌道誤差相位,空間上表征為一定的線性特性,可包含于φatm中;φatm為大氣延遲相位;φnoise為噪聲相位。

對(duì)形變前后影像進(jìn)行干涉處理,引入?yún)⒖糄EM模擬地形相位,去平、降噪、減弱大氣延遲后,剔除模擬地形相位,即可獲取形變相位從而生成差分干涉圖,反演一維地表視向形變信息[13]。

相比D-InSAR技術(shù)的短期形變場(chǎng)獲取,多時(shí)相合成孔徑雷達(dá)干涉測(cè)量技術(shù)(Multi-temporal Synthetic Aperture Radar Interferometry, MT-InSAR)增大影像數(shù)量,通過空間域、時(shí)間域?yàn)V波,緩解大氣延遲效應(yīng)等誤差,能對(duì)長(zhǎng)時(shí)間序列下的影像集進(jìn)行綜合分析,得到地表長(zhǎng)期演變趨勢(shì)[14]。永久散射體技術(shù)(PS-InSAR)及小基線集技術(shù)(SBAS-InSAR)是2個(gè)目前比較有代表性的MT-InSAR算法[15]。

2.1 SBAS-InSAR

SBAS-InSAR技術(shù)將所有覆蓋同一地區(qū)的SAR影像組成若干個(gè)子集,子集內(nèi)的影像基線距(包括時(shí)間基線距和空間基線距)小,子集間的影像基線距大,運(yùn)用奇異值分解方法獲取最小形變速率,基于最小形變速率標(biāo)準(zhǔn)獲得小基線干涉圖。SBAS-InSAR方法限制了長(zhǎng)基線導(dǎo)致的幾何去相干,而且使更多的SAR圖像參與到形變計(jì)算,增加了時(shí)間上的采樣。提出SBAS方法是為解決影像數(shù)量不足、基線失相干無法進(jìn)行相位分離的問題,基于一定的空間時(shí)間基線的若干干涉對(duì)集進(jìn)行后續(xù)相位解纏等操作,從而產(chǎn)生解纏相位圖用以分析地表形變情況,以減少時(shí)間和空間去相干因素的影響,可以得到監(jiān)測(cè)區(qū)域的平均形變速率[16]。

SBAS-InSAR主要處理流程(圖2)包括:估算影像對(duì)基線,選擇配準(zhǔn)的主影像,完成時(shí)序影像配準(zhǔn);按照短空間基線組合原則,生成短基線干涉對(duì)數(shù)據(jù)集,進(jìn)行常規(guī)差分干涉步驟,去除平地、地形相位,獲取差分干涉圖集;利用相干系數(shù)圖,選取高相干點(diǎn),進(jìn)行相位解纏;通過相干性設(shè)定閾值進(jìn)行分布式目標(biāo)分析提取;采用奇異值分解法模型求解形變參數(shù)和高程誤差值即殘余地形的最小范數(shù)解;進(jìn)行高通、低通濾波進(jìn)行第二次反演,計(jì)算非線性形變和大氣相位并進(jìn)行分離,從而反演形變總和,將結(jié)果進(jìn)行地理編碼生成時(shí)序形變信息[17]。

圖2 SBAS-InSAR方法基本流程圖Fig. 2 Basic flowchart of SBAS-InSAR

2.2 PS-InSAR

PS-InSAR技術(shù)選取影像集采用長(zhǎng)時(shí)間序列下各個(gè)分辨單元內(nèi)部保持穩(wěn)定、散射特性強(qiáng)烈的永久散射體,利用相位信息可靠的特征,對(duì)這些離散點(diǎn)進(jìn)行時(shí)間序列分析,從而提取目標(biāo)形變量,探測(cè)地表形變信息[18]。

PS-InSAR技術(shù)(圖3)利用K+1幅SAR單視復(fù)數(shù)影像集,結(jié)合外部DEM,計(jì)算得到K幅差分干涉圖,并通過幅度離差系數(shù)法等算法提取穩(wěn)定的N個(gè)PS候選點(diǎn)[19]。

圖3 PS-InSAR方法基本流程圖Fig. 3 Basic flowchart of PS-InSAR

根據(jù)大氣相位時(shí)序隨機(jī)、空間相關(guān)的特點(diǎn),構(gòu)建差分網(wǎng)絡(luò),對(duì)相鄰的兩PS點(diǎn)干涉相位進(jìn)行3次差分,以降低大氣噪聲相位影響。利用二維周期圖估計(jì)高程以及形變參數(shù),根據(jù)求解結(jié)果進(jìn)行相位解纏,得到的相位殘差中包括殘余大氣信號(hào)、非線性形變信號(hào)。由于大氣相位空間域上低頻信號(hào)、時(shí)間域上高頻信號(hào),進(jìn)行空間域低通濾波、時(shí)間域高通濾波,將大氣相位與非線性形變相位分離,將線性形變相位與非線性形變相位相加即為視向形變相位,最終反演形變場(chǎng)[20]。

3 結(jié)果分析

3.1 D-InSAR沉降結(jié)果分析

對(duì)2020年6月2日—2020年7月8日的2景哨兵-1A影像進(jìn)行D-InSAR處理,生成差分干涉圖(圖4)。結(jié)果表明,研究區(qū)西北部榮華皂村、栗莊村、永定莊村、里南溝村等地可見明顯沉降漏斗,最大視向沉降量達(dá)4.6 cm/36 d,沉降區(qū)域普遍形變量為2 ～3 cm/36 d;東南部城市區(qū)域地表穩(wěn)定,無明顯形變。

圖4 D-InSAR監(jiān)測(cè)研究區(qū)沉降結(jié)果圖Fig. 4 Subsident results of D-InSAR monitoring area

3.2 SBAS-InSAR沉降結(jié)果分析

2019年1月—2020年12月的24景哨兵-1A影像進(jìn)行SBAS-InSAR處理,獲得平均形變速率圖(圖5)。研究區(qū)西北部存在12處明顯沉降,分布范圍基本與D-InSAR結(jié)果重合,且由外緣至中心形變量逐漸增大,沉降區(qū)最大平均沉降速率達(dá)180 mm/y,其余區(qū)域平均沉降速率為5～0 cm/y,2019—2020年最大累積沉降量為333 mm,存在沉降中心超出形變測(cè)量梯度的情況;東南部城市區(qū)域地表穩(wěn)定,無明顯形變。

圖5 SBAS-InSAR技術(shù)監(jiān)測(cè)研究區(qū)平均形變速率圖Fig. 5 Average deformation rate of SBAS-InSAR monitoring area

根據(jù)所屬地理位置將沉降漏斗劃分為榮華皂村沉降區(qū)、晉華宮街道沉降區(qū)、忻州窯村沉降區(qū)、曹家窯村沉降區(qū)、里南溝村沉降區(qū)以及永定莊村沉降區(qū)6個(gè)沉降區(qū)域,其中里南溝村沉降區(qū)沉降現(xiàn)象最為嚴(yán)重,沉降區(qū)域面積達(dá)1.27 km2。

選取里南溝北西向(C—D)、北東向(A—B)剖面做沉降漏斗分析,得到其空間累積沉降分布圖(圖6)。A—B剖面(圖6(b))相對(duì)完整,即累積沉降剖面線由北部邊緣-中心-南部邊緣呈淺-深-淺漏斗狀,該剖面沉降中心隨時(shí)間變化基本維持不變,各處沉降量隨時(shí)間變化而增加,總體不斷下沉;C—D剖面(圖6(c))不完整,邊緣至中心沉降逐步加深,由累積沉降剖面可知,2019年沉降中心尚可監(jiān)測(cè),隨著時(shí)間推移,該剖面沉降中心逐步向東轉(zhuǎn)移,預(yù)計(jì)未來沉降中心沉降值將超出探測(cè)梯度。

(a).里南溝村沉降區(qū)剖面位置圖;(b). A—B剖面示意圖;(c).C—D剖面示意圖圖6 2019—2020年研究區(qū)累積沉降剖面示意圖Fig. 6 Schematic diagram of cumulative subsidence profile from 2019 to 2020 in the study area

3.3 PS-InSAR沉降結(jié)果分析

對(duì)2020年1月—9月的25景哨兵-1A影像進(jìn)行PS-InSAR處理,獲得平均形變速率圖(圖7)。研究區(qū)西北部存在12處明顯沉降,分布范圍基本與SBAS-InSAR結(jié)果重合,從形變分布上可以看出,兩者所獲取的形變區(qū)域空間分布特征一致,沉降漏斗位置分布吻合,且形變趨勢(shì)及累積形變量較為一致。

圖7 PS-InSAR技術(shù)監(jiān)測(cè)研究區(qū)平均形變速率圖Fig. 7 Average deformation rate of PS-InSAR monitoring area

3.4 SBAS與PS技術(shù)結(jié)果交叉驗(yàn)證

在3個(gè)明顯沉降區(qū),選取里南溝村沉降區(qū)邊緣至沉降中心的特征點(diǎn)進(jìn)行時(shí)序?qū)Ρ确治?采用2020年的SBAS-InSAR監(jiān)測(cè)結(jié)果,并提取PS-InSAR中對(duì)應(yīng)時(shí)相的監(jiān)測(cè)結(jié)果,將2020年1月設(shè)為起始月份,繪制特征點(diǎn)2020年1月—2020年9月的累積形變曲線(圖8)。

(a). 里南溝村沉降區(qū);(b). A-1點(diǎn)形變量;(c). A-2點(diǎn)形變量;(d). A-3點(diǎn)形變量圖8 里南溝村沉降區(qū)域3個(gè)特征點(diǎn)累積形變圖Fig. 8 Cumulative deformation of three characteristic points in area of Linangou Village

由里南溝村沉降區(qū)域3個(gè)特征點(diǎn)累積形變曲線圖(圖8)可看出,SBAS-InSAR與PS-InSAR監(jiān)測(cè)結(jié)果總體趨勢(shì)一致,表明2020年以來,該區(qū)域持續(xù)沉陷,其中A-1特征點(diǎn)SBAS-InSAR技術(shù)監(jiān)測(cè)累積形變量最大達(dá)35.93 mm,PS-InSAR技術(shù)監(jiān)測(cè)累積形變量最大達(dá)41.99 mm;A-2特征點(diǎn)SBAS-InSAR技術(shù)監(jiān)測(cè)累積形變量最大達(dá)49.88 mm,PS-InSAR技術(shù)監(jiān)測(cè)累積形變量最大達(dá)69.17 mm;A-3特征點(diǎn)SBAS-InSAR技術(shù)監(jiān)測(cè)累積形變量最大達(dá)57.71 mm,PS-InSAR技術(shù)監(jiān)測(cè)累積形變量最大達(dá)70.99 mm。用均方根誤差(RMSE)對(duì)A-1、A-2、A-3點(diǎn)的形變量觀測(cè)值xsbas,t和xps,t進(jìn)行比較:

計(jì)算得A-1點(diǎn)RMSE值為2.89,A-2點(diǎn)RMSE值為7.48,A-3點(diǎn)RMSE值為5.31。

對(duì)兩組形變量觀測(cè)數(shù)據(jù)做線性擬合(圖9),A-1觀測(cè)點(diǎn)的決定系數(shù)R2為0.947 2,A-2觀測(cè)點(diǎn)的決定系數(shù)R2為0.960 1,A-3觀測(cè)點(diǎn)的決定系數(shù)R2為0.962 6。

圖9 里南溝村沉降區(qū)域3個(gè)特征點(diǎn)兩組形變量觀測(cè)結(jié)果擬合曲線圖Fig. 9 Fitted curve of the two observation groups of shape variables at three characteristic points in area of Linangou Village

由上可知,SBAS-InSAR形變監(jiān)測(cè)與PS-InSAR形變監(jiān)測(cè)結(jié)果相近,形變區(qū)域吻合,形變趨勢(shì)基本一致,決定系數(shù)均>0.9,驗(yàn)證了SBAS與PS兩種時(shí)序InSAR技術(shù)在區(qū)域形變監(jiān)測(cè)具有一致性。

4 結(jié)論

大同市云岡區(qū)作為全國最大的產(chǎn)煤區(qū)之一,其地下煤礦開采區(qū)引起地面沉降災(zāi)害。本文應(yīng)用InSAR技術(shù)對(duì)區(qū)域的多景哨兵-1A雷達(dá)影像進(jìn)行處理,獲取了地面沉降特征,詳細(xì)分析了時(shí)空演變規(guī)律并進(jìn)行了初步探討,主要結(jié)論如下:

(1) 沉降主要分布于云岡區(qū)西北部煤礦開采區(qū),沉降面積大,沉降速率達(dá)180 mm/y,最大累積沉降量為333 mm;東南部城市區(qū)域沒有明顯形變跡象。

(2) 采用SBAS-InSAR與PS-InSAR技術(shù)對(duì)其地表沉降情況進(jìn)行長(zhǎng)時(shí)間序列監(jiān)測(cè),監(jiān)測(cè)結(jié)果表明形變空間分布吻合,時(shí)序沉降趨勢(shì)一致,形變速率具有較高的相關(guān)性(R2決定系數(shù)均>0.9),均能有效監(jiān)測(cè)地表沉降變化。