王 羽 王 磊 李靖宇 江克貴 李世保 楊 科 劉洪江

(1.安徽理工大學空間信息與測繪工程學院,安徽 淮南 232001;2.安徽理工大學礦山采動災害空天地協(xié)同監(jiān)測與預警安徽普通高校重點實驗室,安徽 淮南 232001;3.安徽理工大學礦區(qū)環(huán)境與災害協(xié)同監(jiān)測煤炭行業(yè)工程研究中心,安徽 淮南 232001;4.中國礦業(yè)大學(北京)地球科學與測繪工程學院,北京 100083;5.安徽理工大學礦業(yè)工程學院,安徽 淮南 232001)

地下礦產(chǎn)資源開采會使采空區(qū)上方巖層受到重力影響而產(chǎn)生形變,易導致位于采空區(qū)上方的建(筑)物、民房、耕地以及交通設施等遭到破壞,還會引發(fā)地質災害,影響了區(qū)域居民生產(chǎn)生活,并對礦區(qū)安全生產(chǎn)造成了很大威脅[1-5]。通常情況下,礦區(qū)的地表形變主要表現(xiàn)為地面沉降。傳統(tǒng)的礦區(qū)地表形變觀測技術主要包括地面觀測技術和三維激光掃描技術。其中,地面觀測技術主要在地面布設觀測站,采用全站儀、水準儀、GPS-RTK[6-7]等測量儀器獲取觀測站各點的形變信息,再根據(jù)監(jiān)測點上的形變來反演整個礦區(qū)的地表形變。該類方法雖然觀測精度很高,但局限性較明顯,即無法動態(tài)地反映整個區(qū)域的下沉狀態(tài),并且人力物力成本較高,不利于對一個區(qū)域進行長時間的追蹤觀測。三維激光掃描技術可以獲得實時動態(tài)的地表形變信息,且觀測精度可達到毫米級別,但是三維激光點云數(shù)據(jù)處理耗時費力、監(jiān)測效率低、觀測成本較高,使得該技術難以大范圍推廣應用[8]。受到多方面的條件限制,傳統(tǒng)的礦區(qū)地表形變觀測技術已經(jīng)無法適應現(xiàn)代礦區(qū)災害預警體系中對于地表形變?nèi)媲腋咝ПO(jiān)測的要求。

20世紀60年代以來,合成孔徑雷達技術(Synthetic Aperture Radar,SAR)的發(fā)展使得地表形變監(jiān)測的研究迎來重大技術性突破,由于雷達波束的特殊性,使得該技術可不受氣候條件和晝夜的影響,甚至可以透過部分地表或者植被獲取地面信息[9-10]。該類特性使其在地震、冰川和地表形變監(jiān)測等領域具有非常廣闊的應用場景,尤其是近些年來迅速發(fā)展的差分干涉測量技術(Differential Interferometry Synthetic Aperture Radar,D-InSAR)能夠以較低成本獲取大區(qū)域較高精度的地表形變信息[11]。國外學者較早應用該技術對地表形變進行研究,并取得了較好的成果。1989年,GABRIEL等[12]應用D-InSAR技術成功測量了加利福尼亞河谷的地表形變,并且監(jiān)測精度可達到厘米級,該項研究第一次證實了D-InSAR技術可實現(xiàn)地表形變的高精度監(jiān)測。但是,D-InSAR技術在當時并未引起太大的反響,直至1993年,CLAUDIE等[13]使用雷達數(shù)據(jù)提取了美國加利福尼亞州蘭德斯地震的位移場,并且通過與實測數(shù)據(jù)比較證明了D-InSAR監(jiān)測結果比較準確,至此該技術開始受到學術界的重視。國內(nèi)在D-InSAR形變觀測方面的研究雖然起步比較晚,但也取得了很好的成績。在2000年前我國學者主要是針對合成孔徑雷達干涉測量(Interferometry Synthetic Aperture Radar,InSAR)技術的一些關鍵算法和應用進行了展望和試驗研究,具體的技術應用案例比較少。2001年,劉國祥等[14]使用衛(wèi)星雷達數(shù)據(jù)獲取了香港近?；靥畹鼗? a內(nèi)的地表沉降信息,其結果精度優(yōu)于1 cm。2011年,吳宏安等[15]使用ENVI ASAR影像數(shù)據(jù)獲取了太原市2003—2009年的地表形變,論證了差分干涉技術可以應用于城市地表形變監(jiān)測。2016年張永紅等[16]提出了多主影像相干目標小基線干涉技術,并使用該技術提取了京津冀地區(qū)22 a間的地面沉降信息,并從時間和空間角度分析了地表形變的變化特征。

研究發(fā)現(xiàn):對于礦區(qū)形變梯度比較小的區(qū)域,DInSAR能夠有效獲取其真實的地表形變量,但是當?shù)V區(qū)形變梯度過大而超過了D-InSAR的監(jiān)測能力時,會導致該技術無法獲取礦區(qū)真實的地表形變信息。這是因為當形變梯度過大時,D-InSAR在進行相位解纏時會使整周相位丟失,從而使得最終解算的形變量遠小于實際形變量。MASSONNET等[17]認為雷達數(shù)據(jù)分辨率、波長等都影響D-InSAR技術對地面沉降的監(jiān)測能力,當同一位置的像元在某一時段內(nèi)的形變量過大,且當該處像元在進行干涉測量時,由于相干性較低甚至出現(xiàn)失相干的現(xiàn)象,會導致該處像元無法直接通過D-InSAR技術獲取真實準確的形變值。蔣彌等[18]研究發(fā)現(xiàn),最大形變梯度和形變量不僅與相干性有關,多視的視數(shù)對可監(jiān)測最大形變梯度和最大形變量的能力也有影響。冷紅偉[19]對比分析了TerraSAR-X雷達數(shù)據(jù)和Sentinel-1A雷達數(shù)據(jù)分別在多視和全分辨率情況下監(jiān)測大梯度形變的能力。

對于上述存在的問題,不少學者從不同角度提出了解決方法,其中在1999年,MICHEL等[20]提出了像素偏移量估計(Offset-tracking)技術,該技術利用SAR影像強度信息,通過影像配準原理進行互相關計算從而提取了像元間的偏移量。像素偏移量估計技術對于雷達影像的相干性要求很低,可避免相位解纏錯誤的問題,在一定程度上可實現(xiàn)大梯度地表形變監(jiān)測,但是該技術監(jiān)測精度較低,無法滿足實際應用需求,并且在積水區(qū)域雷達的強度信息減弱,這些因素都導致了像素偏移量估計技術在地表形變監(jiān)測中的應用受到制約。蔣創(chuàng)等[21]結合開采沉陷地表動態(tài)移動變形機理,研究了基于動態(tài)概率積分法模型(Dynamic Probability Integral Method,DPIM)約束的InSAR相位解纏模型和地表三維變形提取方法,克服了礦區(qū)大梯度形變區(qū)域解纏失敗的不足,恢復了密集區(qū)纏繞相位,在較大程度上實現(xiàn)了礦區(qū)開采沉陷盆地三維形變提取。

針對D-InSAR因礦區(qū)形變梯度過大而無法獲取真實地表形變量的問題,實現(xiàn)礦區(qū)全尺度梯度地表形變監(jiān)測,本研究提出了一種融合水準和D-InSAR的礦區(qū)地表形變監(jiān)測方法(LEV-InSAR模型)。該模型結合D-InSAR地表形變監(jiān)測的基本原理,通過水準測量技術彌補D-InSAR相位解纏失敗的不足,獲取礦區(qū)完整連續(xù)且準確的地表形變信息。為了驗證該方法的監(jiān)測精度和適用性,針對山西省長治市霍爾辛赫礦區(qū)3501工作面,選取9景Sentinel-1A SAR影像,采用雙軌D-InSAR技術獲取2021年12月23日—2022年3月29日的礦區(qū)地表形變監(jiān)測結果,并結合同時期觀測線水準觀測數(shù)據(jù),成功解算了工作面內(nèi)準確的地表形變量,取得了良好的監(jiān)測效果。本研究提出的監(jiān)測方法對于礦區(qū)大梯度形變監(jiān)測具有一定的實用價值。

1 原理與方法

1.1 雙軌D-InSAR技術基本原理與流程

雙軌D-InSAR技術基本原理是利用同一區(qū)域不同時相的兩景雷達影像,依據(jù)雷達影像中包含的相位信息進行干涉處理,得到干涉相位圖。這些干涉相位中包含多種相位信息(式(1)),通過各種方法解算出形變相位,再將形變相位按照一定的關系式轉化成形變量。

式中,Φtopo為地形相位;Φflat為平地效應相位;Φdef為地表形變相位;Φatm為大氣延遲相位;Φnoise為噪聲相位。

得到真實地表形變相位Φdef后,進而計算雷達沿視線方向的地表形變ΔR,公式為

式中,λ為雷達波長,mm。

雙軌D-InSAR基本流程如圖1所示。

圖1 雙軌D-InSAR技術流程Fig.1 Flow of dual-track D-InSAR technique

1.2 LEV-InSAR方法基本原理與流程

在實際應用中,當影像相鄰像元的形變梯度過大時,D-InSAR的整周模糊數(shù)相位解纏結果不正確,導致丟失若干整周相位(解算的相位差與真實相位差之間相差2π的整數(shù)倍),從而使最終的D-InSAR監(jiān)測形變結果與真實量之間相差若干個整周相位所對應的形變量。理論上,通過D-InSAR得到的地表形變監(jiān)測值中包含的不足整周部分相位對應的形變量是較為準確的。

針對D-InSAR監(jiān)測方法的不足,提出了一種融合水準數(shù)據(jù)的D-InSAR礦區(qū)全尺度梯度地表形變監(jiān)測方法(LEV-InSAR模型)。該模型從觀測線水準數(shù)據(jù)出發(fā),借鑒概率積分法思想,應用觀測線擬合曲線方程得到工作面上方地表任意點的模擬形變值,由于模擬形變值與地表真實形變值的量級相當,實際上,根據(jù)模擬形變值可以提取出D-InSAR監(jiān)測像元的實際整周模糊數(shù),最后整合解算的整周相位形變量和去除整周相位后的D-InSAR監(jiān)測值,即可獲得開采工作面全尺度梯度地表形變監(jiān)測信息。具體思路是:首先在礦區(qū)工作面上方地表以走向線和傾向線為基準構建XY工作面坐標系,并將水準監(jiān)測點所在的大地坐標系轉換到工作面坐標系內(nèi);然后將水準測量值沿工作面走向和傾向方向擬合得到兩個擬合曲線方程,即傾向線擬合曲線方程和走向線擬合曲線方程,記為W(x)和W(y)。

根據(jù)概率積分法理論,對于工作面內(nèi)任意一點的形變值可表示為

式中,k為系數(shù),通過最小二乘原理擬合得到。

通過式(3)估算得到的任一點形變值W(x,y)并未進行相位纏繞,則可以通過該值計算出工作面任意一點真正的整周模糊數(shù),再由整周模糊數(shù)計算得到形變量WZ(x,y):

式中, 一個整周相位相當于1/2的雷達波長量級的形變;int·()為取整函數(shù)。

通過D-InSAR獲得的地表任意一點的形變值記為WD(x,y),當?shù)V區(qū)形變梯度過大時,由于影像失相干導致在相位解纏時丟失了若干整周相位,使得形變值WD(x,y)與真實值之間相差若干個整周相位量級的形變量,然而形變值WD(x,y)中包含的不足整周部分相位對應的形變量較為準確。一般來說,WD(x,y)中也可能含有部分整周相位形變量,為了不與WZ(x,y)中部分形變量重復,應從WD(x,y)中去除整周相位形變量,只保留不足整周部分相位對應的形變值,記為WD′(x,y)。

其中,n為非負整數(shù);WDi(x,y)為任意一點處第i期的D-InSAR監(jiān)測值,mm;m為D-InSAR監(jiān)測總期數(shù);WDi′(x,y)為經(jīng)過去除整周相位后的任意一點第i期D-InSAR監(jiān)測值,mm;WD′(x,y)為經(jīng)過去除整周相位后任意一點的D-InSAR累計監(jiān)測值,mm;“±”號由WDi(x,y)值正負決定,當WDi(x,y)為正時,取“-”號,當WDi(x,y)為負時,取“+”號。

將通過水準數(shù)據(jù)模擬得到的整周相位形變值WZ(x,y)與從D-InSAR監(jiān)測值中提取的不足整周部分相位對應的形變值WD′(x,y)求和,即可得到工作面內(nèi)任意位置的地表形變值W′(x,y):

LEV-InSAR模型的具體流程如圖2所示。

圖2 LEV-InSAR模型的具體流程Fig.2 Detailed flow of LEV-InSAR model

2 工程實例

2.1 研究區(qū)概況

選取山西省長治市長子縣霍爾辛赫礦區(qū)3501工作面作為研究區(qū)域。該工作面位于長子縣城區(qū)東北方向,距離長子縣主城區(qū)約3.4 km,地理坐標為東經(jīng)112°54′30″～112°55′10″,北緯36°7′40″～36°8′10″。工作面東西寬約300 m,南北長約800 m,占地面積約0.24 km2。3501工作面開采煤層平均厚度為5.6 m,煤層傾角為1°～12°,平均開采深度為576 m,工作面日推進速度為3～4刀,約為2.4 m/d。3501工作面于2022年2月初開始回采,2022年2月8日開始離層注漿,2022年6月15日該工作面基本停采,2022年9月15日注漿停止。工作面內(nèi)地形起伏較小,地表高程為935～939 m。該工作面內(nèi)EW向縣道橫穿工作面,道路南側有一焦化廠,其余部分為大面積的農(nóng)田。工作面內(nèi)共設置了兩條觀測線,一條為EW向的傾向觀測線,點號標記為Q01～Q22,另一條為SN向的走向觀測線,點號標記為Z01～Z27,共49個觀測站。工作面位置及觀測線布置如圖3所示。

圖3 3501工作面位置及觀測線布置Fig.3 3501 working face position and layout of observation lines

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

2.2.1 Sentinel-1A SAR影像數(shù)據(jù)

根據(jù)3501工作面的開采時間和進度,并且為了避免工作面上方夏秋季節(jié)農(nóng)田植被對D-InSAR干涉相位造成較大影響,選定2021年12月—2022年3月為研究時段。在此基礎上選取了覆蓋研究區(qū)域的9景C波段Sentinel-1A 升軌數(shù)據(jù),數(shù)據(jù)類型為單視復數(shù)影像(Single Look Complex,SLC),成像模式為干涉寬幅(Interferometric Wide swath,IW),極化方式為VV(垂直傳送和垂直接收),成像日期為2021年12月23日—2022年3月29日,影像的具體參數(shù)取值見表1,數(shù)據(jù)覆蓋范圍如圖4所示。DEM數(shù)據(jù)來源于美國航空航天局(National Aeronautics and Space Administration,NASA)和美國國家地理空間情報局(National Geospatial-Intelligence Agency,NGA)免費對外開放的航天飛機雷達地形測繪任務(Shuttle Radar Topography Mission,SRTM)數(shù)據(jù),SRTM DEM數(shù)據(jù)的空間分辨率為30 m,絕對高程精度為±16 m,DEM數(shù)據(jù)將用于去除差分干涉測量過程中由地形起伏因素引起的地形相位。

圖4 Sentinel-1A SAR影像覆蓋范圍Fig.4 Coverage of Sentinel-1A SAR image

表1 Sentinel-1A SAR影像參數(shù)Table 1 Sentinel-1A SAR image parameters

2.2.2 水準數(shù)據(jù)

本研究依據(jù)9景Sentinel-1A SAR影像成像時間,選取該工作面2021年12月22日—2022年3月29日走向線和傾向線上地表移動觀測站的觀測數(shù)據(jù)。該觀測數(shù)據(jù)來源于17期三等水準測量的實測數(shù)據(jù),第1期觀測時間為2021年12月22日,第17期觀測時間為2022年3月28日,平均每5 d觀測1次,測量儀器為萊卡sprinter100數(shù)字水準儀,標稱精度為每千米往返測高程中誤差0.4 mm,采用平差軟件進行水準網(wǎng)平差,平差后最弱點中誤差為2.93 mm。原定觀測方案中共布設了49個觀測站,但在后期觀測過程中,部分觀測站因外部因素影響而遭到破壞。為保證水準數(shù)據(jù)的可靠性和有效性,本研究在選取觀測站數(shù)據(jù)時,只保留了部分觀測站水準數(shù)據(jù),其中包括20個傾向線觀測站和26個走向線觀測站,共計46個觀測站的水準觀測數(shù)據(jù)。

2.3 結果與分析

2.3.1 D-InSAR地表形變監(jiān)測結果分析

對獲取的9景Sentinel-1A數(shù)據(jù)進行雷達差分干涉測量處理,按照時間序列構建8個干涉對,每個干涉對時間基線均為12 d。按照D-InSAR實用化數(shù)據(jù)處理流程對8個干涉對進行處理,得到相鄰影像之間的形變信息,通過統(tǒng)計礦區(qū)各時間段的沉降信息,得到該工作面2021年12月23日—2022年3月29日的累計形變圖以及觀測線下沉曲線如圖5所示。

圖5 D-InSAR地表形變監(jiān)測結果Fig.5 Monitoring results of ground deformation by D-InSAR

由圖5可知:該礦區(qū)在傾向線以北發(fā)生的形變較為明顯,沉降中心位于走向觀測點Z06附近,沉降分布范圍及形態(tài)類似于一個漏斗,在沉降中心區(qū)域沉降量最大,邊緣部分沉降量較小。2021年12月23日—2022年3月29日,D-InSAR所監(jiān)測的最大累計沉降量為172.5 mm,與水準實測數(shù)據(jù)相比,D-InSAR地表形變監(jiān)測結果的均方根誤差為112.0 mm。在沉降漏斗邊緣及外部區(qū)域,D-InSAR監(jiān)測到的累計沉降量小于20 mm。這些微量形變形成原因為:① 受到礦區(qū)開采影響而造成的微小形變;② 可能是衛(wèi)星接收雷達信號時產(chǎn)生的誤差以及衛(wèi)星本身固有的誤差;③ 數(shù)據(jù)處理過程中產(chǎn)生的誤差。

為了從數(shù)值上反映D-InSAR地表形變監(jiān)測效果,按照觀測點位置提取了若干個觀測點的D-InSAR地表形變監(jiān)測值,并與水準實測數(shù)據(jù)進行了對比,結果見表2。

表2 D-InSAR地表形變監(jiān)測誤差統(tǒng)計Table 2 Statistics of D-InSAR surface deformation monitoring errors

由表2可知:在礦區(qū)沉降盆地邊緣,地表形變量較小,D-InSAR監(jiān)測結果與水準監(jiān)測結果吻合程度較高,但隨著沉降程度的增加,在沉降嚴重的水準點上,D-InSAR監(jiān)測結果與水準監(jiān)測結果之間的絕對差值明顯增加。這一現(xiàn)象很好地說明了在地面形變梯度過大區(qū)域,差分干涉測量的相位解纏錯誤使得D-InSAR監(jiān)測結果與真實地表形變信息之間有較大差異,并且這種差異在沉降盆地中心表現(xiàn)得更為明顯。差分干涉測量技術雖然在礦區(qū)形變監(jiān)測中具有顯著優(yōu)勢,但從監(jiān)測結果來看,其更適用于礦區(qū)周邊緩慢微小形變監(jiān)測,對于形變梯度大的沉降中心區(qū)域監(jiān)測能力明顯不足。

2.3.2 LEV-InSAR模型監(jiān)測結果分析

基于本研究提出的融合D-InSAR和水準的地表形變監(jiān)測方法獲取3501工作面地表形變信息,按照LEV-InSAR模型的基本思路,首先需要構建LEVInSAR模型監(jiān)測方程。根據(jù)觀測站點位信息,構建基于剖面觀測線的XOY標系(圖6),賦予各觀測站點橫(x)縱(y)坐標。

圖6 XOY坐標系示意Fig.6 Schematic of XOY of coordinate system

按照觀測線水準觀測數(shù)據(jù),分別擬合得到走向和傾向觀測線擬合曲線方程:

式中,x0,xc,ω1,A1,y0,yc,ω2,A2均為刻畫觀測站點坐標與其形變值之間函數(shù)關系的參數(shù),無實際含義,當坐標信息和形變值確定時,這些參數(shù)均為實際值。

根據(jù)MATLAB非線性擬合得到的式(7)中各參數(shù)取值分別為:x0=-4.860 72,xc=543.543 81,ω1=160.062 82,A1=-136 062.176 7;y0=-32.115 81,yc=321.732 65,ω2=110.944 8,A2=-117.203 99。

由式(3)可得工作面內(nèi)任意一點處的形變值估算方程為

式中,k取-230。

根據(jù)式(5),分別對每一期D-InSAR監(jiān)測結果進行甄別,經(jīng)過去除整周相位形變量后累加,得到WD′(x,y),再根據(jù)式(4)和式(6)得到本研究工程實例的LEV-InSAR模型監(jiān)測方程為

其中,波長λ取56 mm。

按照式(9)模型,融合水準觀測數(shù)據(jù)和D-InSAR形變量解算得到3501工作面的LEV-InSAR模型地表形變監(jiān)測結果,如圖7所示。

圖7 LEV-InSAR模型地表形變監(jiān)測結果Fig.7 Surface deformation monitoring results by LEV-InSAR model

為驗證該方法的監(jiān)測精度,分別提取了兩條觀測線的下沉曲線,如圖8所示。此外,統(tǒng)計了各觀測點的LEV-InSAR模型地表形變監(jiān)測值相對于水準觀測數(shù)據(jù)的絕對誤差和相對誤差,見表3。

圖8 主斷面下沉曲線Fig.8 Subsidence curves of main section

表3 LEV-InSAR模型監(jiān)測誤差統(tǒng)計Table 3 Statistics of monitoring errors LEV-InSAR model

由圖8、表3可知:傾向線和走向線上觀測點的監(jiān)測結果與實際測量的下沉值基本吻合,傾向線和走向線的下沉曲線和水準曲線的形態(tài)特征趨近。與水準數(shù)據(jù)相比,傾向線上各觀測點誤差為-16.9～16.7 mm,最大監(jiān)測誤差為16.9 mm,約為最大沉降值的3.3%,最小監(jiān)測誤差為0.8 mm,傾向線均方根誤差為10.2 mm,平均相對誤差為14.4%;走向線上各觀測點誤差為-32.0～41.0 mm,且大部分觀測點的監(jiān)測誤差小于20 mm,最大監(jiān)測誤差為41.0 mm,約為最大沉降值的8%,最小監(jiān)測誤差為0.7 mm,走向線均方根誤差為15.8 mm,平均相對誤差為26.6%。從單點誤差來看,由于個別點實際下沉值較小,使得這些觀測點的相對誤差對于絕對誤差比較敏感,使得該類觀測點的相對誤差和整體的平均相對誤差較大。對于累計沉降量大的觀測點而言,相對誤差較小,由此可見,LEV-InSAR模型對于礦區(qū)大梯度地表形變監(jiān)測具有更明顯的優(yōu)勢。

3 結 語

針對D-InSAR無法實現(xiàn)礦區(qū)全尺度梯度地表形變監(jiān)測的不足,提出了一種融合D-InSAR和水準的礦區(qū)全尺度梯度地表形變監(jiān)測方法(LEV-InSAR模型)。研究表明:LEV-InSAR模型對于礦區(qū)大梯度地表形變監(jiān)測具有較好的適用性和可靠性,可以提取到礦區(qū)空間連續(xù)的真實地表形變信息,具有較好的應用價值。