全分辨率與多視法DInSAR礦區地面沉降監測研究

李賓， 陶秋香， 牛沖， 劉曉帥， 張亞鳳

1. 山東科技大學 測繪與空間信息學院,山東 青島 266590； 2. 山東省地質測繪院,濟南 250014

地面沉降主要是一種在自然和人為因素的作用下,使地下松散土層及巖層壓縮而導致地面高程逐漸降低的工程地質現象.隨著社會經濟飛速發展,人們在過度開采煤炭、石油、天然氣、地下水等資源以滿足自身需要的同時[1-2],地面沉降的程度和范圍正在進一步加深和加大[3-4],這已成為經濟社會可持續發展的重要制約因素.合成孔徑雷達差分干涉測量(Differential Interferometric Synthetic Aperture Radar,DInSAR ) 技術具有全天候、全天時、高分辨率和監測范圍廣等優勢,能夠應用于礦區地面沉降監測[5].1989年Grabriel等[6]利用Seasat衛星L波段的SAR影像首次論證了 DInSAR可用于厘米級的地表形變監測.而后,國內外研究人員開展了DInSAR監測礦區地面沉降的相關理論與應用研究,取得了一系列研究成果[7-12].現有研究大多在差分干涉測量之前進行了多視處理,以抑制SAR影像的斑點噪聲,進而獲取礦區地面沉降信息,即多視法DInSAR.然而,多視處理會降低影像的分辨率,不可避免地會降低InSAR監測到的最大沉降梯度和沉降量[13-14].全分辨率DInSAR則對SAR影像不進行多視處理,在保持原有SAR影像分辨率的基礎上進行差分干涉處理,從而提高DInSAR技術監測最大沉降梯度和沉降量的能力．

本研究以山東省某礦區為研究區,選取2018年10月18日-2019年1月22日時間段內6景哨兵數據為數據源,分別用全分辨率和多視法(距離向與方位向之比為5∶1)進行差分干涉處理,并結合礦區同時間段的工作面開采信息和水準數據監測,系統對比分析和研究了全分辨率和多視法DInSAR在礦區地面沉降監測中的應用效果.研究結果表明,由于受斑點噪聲的影響較為嚴重,全分辨率DInSAR在沉降梯度和沉降量都較大的沉降區域,其差分干涉條紋的連續性變差,相位解纏和形變解算的精度降低,對比多視法DInSAR在礦區沉降盆地中心的地面沉降監測結果,并未有明顯改善．

1 研究區概況

研究區位于華北平原的南緣,黃河北岸,地理坐標范圍為116°59′00″-117°06′00″E,36°47′00″-36°50′00″N.圖1給出了礦區的具體位置及范圍.礦區東西向長約10 km,南北向寬約5 km,面積約50 km2.礦區南距濟南市區8 km,北距濟陽縣城20 km,有國道、高速穿過井田,交通運輸便利.礦區地面標高+22～+27 m,地勢平坦,局部地勢低洼,雨季易積水.黃河從井田南部1 km處流過,井田中東部有邢家渡引黃總干渠,中西部有大寺河通過,為本井田主要地表水系．

審號圖：GS(2019)1822號．圖1 研究區地理范圍

2 數據源

2.1 SAR影像數據

在礦區形變監測實際應用中,夏季植被覆蓋茂密,失相干嚴重,對于礦區形變區域的定位和形變監測反演造成很大阻礙[15-16].研究中采用冬季的SAR影像數據進行DInSAR處理,獲取礦區地面沉降信息.本研究選取了覆蓋礦區的時間跨度為2018年10月18日-2019年1月22日的6景Sentinel-1A升軌SAR影像,組成5組干涉對進行差分處理.SAR影像參數如表1所示．

表1 DInSAR差分干涉對

2.2 DEM數據

航天飛機雷達地形測繪使命(Shuttle Radar Topography Mission,SRTM)由美國太空總署(National Aeronautics and Space Administration,NASA)和國防部國家測繪局(National Imagery and Mapping Agency,NIMA)以及德國與意大利航天機構共同合作完成聯合測量[17-18].2000年2月11日,美國發射“奮進”號航天飛機,其載SRTM系統,進行了11 d的數據采集工作,獲取了覆蓋地球80%以上的陸地數字高程模型(Digital elevation model,DEM).SRTM數據是迄今為止現勢性最高、分辨率最高、精度最高的全球數字地形數據.為了相位校正和地理編碼,本研究采用了地面分辨率約90 m的SRTM3 DEM,其絕對高程精度是±16 m ．

2.3 礦區部分工作面及其水準數據

2018年10月至2019年1月,研究區開采主要集中在第三、四采區,根據6景Sentinel-1A SAR影像的時間跨度(2018年10月18日-2019年1月22日),本研究收集、整理和分析了三、四采區5個工作面的開采情況和232個水準點監測數據,用于比較、分析和驗證全分辨率和多視法DInSAR礦區地面沉降監測的準確性和可靠性.圖2顯示了5個工作面和232個水準點的分布和范圍,表2提供了5個工作面的基本開采信息.其中,各水準點的編號、傾向線按自南到北的方向順序進行編號,走向觀測線按工作面推進的方向順序進行編號．

圖2 三、四采區部分工作面及其水準點布設

表2 三、四采區5個工作面的基本信息

3 礦區沉降監測原理與數據處理流程

首先選取覆蓋完整研究區形變前后的兩幅SAR影像,分別作為主、輔影像,根據二者之間的成像幾何關系進行復圖像配準處理,進而計算同名點偏移量,使主、輔影像的像元對應同一地面單元[19-20].對配準后的影像數據作復共軛相乘,獲取干涉條紋圖.此時的干涉圖中既包含目標地物的形變信息,又包括一系列干擾相位信息.為抑制各種噪聲對干涉圖質量的影響,對其進行濾波,繼而引入外部DEM,模擬地形相位,與SAR干涉結果進行差分,消除地形相位,得到由地表形變引起的相位.但由于三角函數的周期性,此時的差分干涉相位纏繞在[-π,π]之間,因此需要進行相位解纏,得到代表真實地表形變的相位?defo,進而計算雷達沿視線方向的地表形變Δr：

(1)

式中：λ為波長；Δr為雷達視線向上的地表形變位移．

雷達視線向上的地表形變位移Δr：

(2)

雷達在垂直方向上的地表形變位移Δrvertical：

(3)

式中：θ為雷達視線的入射角．

技術具體流程如圖3所示．

圖3 數據處理流程圖

4 結果及分析

將表1中5組干涉對的主、輔影像配準,利用外部DEM去除地形相位,再分別利用全分辨率和5∶1多視法進行干涉處理,得到差分干涉圖.然后對差分干涉圖做Goldstein濾波處理,得到濾波后如圖4所示的差分干涉圖．

圖4 研究區全分辨率、多視處理濾波后的差分干涉圖

分析圖4可以看出：

全分辨率與多視法DInSAR所得到的差分干涉條紋的位置和分布非常一致.但從視覺效果上來說,全分辨率的差分干涉條紋在距離向(干涉圖列方向上)被拉伸,可以更好地展示出地表形變的細節信息,從而也造成該方法得到的沉降中心的差分干涉條紋不如多視處理后的差分干涉條紋完整、連續．

圖4d和4e的左側均存在有規律的、類似系統誤差的干涉條紋,這是由SAR影像自身造成的,因為2018年12月17日獲取的SAR影像上相應位置就存在這種相同的條紋.全分辨率與多視法DInSAR均無法將其從差分干涉圖上去除．

用最小費用流法對圖4所示的差分干涉圖進行相位解纏,由差分干涉相位值恢復為真實相位,然后根據公式(2)和(3),計算雷達視線向和垂直向上的地表形變量,最后對其地理編碼,得到如圖5所示的礦區垂直向地面沉降圖.需要說明的是,為了便于對比和分析,圖5中的沉降圖都進行了距離向和方位向(列方向和行方向)上分辨率為15m×15m的重采樣處理．

圖5 全分辨率與多視法DInSAR監測的礦區地面沉降對比

結合表2,分析圖5可以看出：

2018年10月18日至2018年11月11日(圖5a),工作面4701的開采工作剛剛結束,全分辨率DInSAR在該工作面上監測到零星的-10～0 mm的沉降區域,多視法DInSAR 則監測到了較大的、較明顯的-10～0 mm的沉降區域； 工作面3107的開采工作已經結束2個月,而4702的開采工作進行了3個月,全分辨率與多視法DInSAR均監測到這兩個工作面上存在較大的、較明顯的-10～0 mm的沉降區域； 工作面4707的開采工作進行了一個月,全分辨率DInSAR在該工作面上監測到零星的-10～0 mm的沉降區域,多視法DInSAR 則監測到了較大的、較明顯的-10～0 mm的沉降區域； 工作面4705的開采工作進行了7個月,全分辨率與多視法DInSAR均監測到該工作面是主要的沉降區域,已形成沉降盆地,最大沉降的位置分別在a1和b1處,最大沉降值分別為-47 mm和-45 mm.但由于兩幅SAR影像的空間基線(106 m)和時間基線(24 d)都較大,在該沉降盆地,全分辨率與多視法DInSAR均未得到連續的干涉條紋,監測的沉降盆地也比較凌亂、不連續.總之,全分辨率與多視法DInSAR監測的礦區地面沉降與礦區工作面的開采進度相符,但二者監測的沉降值和沉降范圍并不相同．

2018年11月11日至2018年11月23日(圖5b),全分辨率與多視法DInSAR均監測到礦區主要的沉降區在5個工作面上方,礦區的地面沉降仍在持續加重,二者監測到的主要沉降區的位置、沉降量級也比較一致,例如,二者監測的最大沉降位置a1和b1非常一致,而且都監測到該處存在比較明顯的沉降盆地,但沉降量的大小略有不同； 全分辨率與多視法DInSAR都監測到了整個三、四采區都存在-20～0 mm的地面沉降,這有可能是由于DInSAR的兩幅SAR影像數據質量的影響,也有可能是數據處理過程中有些誤差難以去除造成的．

2018年11月23日至2019年1月22日(圖5c-5e),整體而言,全分辨率與多視法DInSAR監測的礦區地面沉降的分布、形態和變化趨勢等與礦區工作面的開采情況相符.但從細節上來看,二者監測到的最大沉降值、最大沉降的位置、范圍的大小等并不相同.分析圖5e發現,由于兩幅SAR影像的時間基線(36 d)較大,兩幅SAR影像的相干性較低,而且礦區沉降較為嚴重,全分辨率DInSAR所得到的差分干涉條紋的連續性較差,在工作面4702上方、4705右上方的沉降較為嚴重的沉降盆地,相位解纏和形變反演結果的連續性較差,沉降監測結果的精度變低.相比較而言,多視法DInSAR的形變監測結果較好,在4705右上方的沉降盆地得到了較為連續、精度較高的地面沉降監測結果.因此,全分辨率與多視法DInSAR監測到的最大沉降值出現的位置相差非常大．

為定量對比分析5個監測時間段內全分辨率與多視法DInSAR監測的礦區地面沉降情況,表3中對全分辨率和多視法DInSAR監測的各時段的最大沉降量、最大沉降位置進行了統計．

表3 全分辨率與多視法DInSAR監測的最大沉降值及其位置統計表

結合圖5,分析表3可以看出：前4個時間段內全分辨率與多視法DInSAR監測的最大地面沉降位置基本一致,全分辨率DInSAR監測的地面沉降最大值均大于多視法DInSAR監測的地面沉降最大值,最大絕對差值為12 mm.在第5個時間段內,全分辨率DInSAR在b5處的沉降盆地沒有得到連續的差分干涉相位和正確的相位解纏結果,從而導致了監測能力和精度降低,其監測到的最大沉降位置a5與b5并不在一個沉降盆地,二者相距甚遠,DInSAR監測的a5處的地面沉降最大值小于多視法DInSAR監測的b5處的地面沉降最大值.總體上看全分辨率方法監測到的沉降最大值更有效．

5 精度驗證與對比分析

為了驗證全分辨率與多視法DInSAR監測的礦區地面沉降結果的精度與可靠性,將5個時間段(2018年10月18日至2019年1月22日)全分辨率與多視法DInSAR監測的礦區地面沉降累計結果與濟陽礦區插值后同時間段經過精密水準測量所獲取的地面沉降結果進行了對比分析.疊加水準點的沉降圖如圖6所示,沿觀測線水準點的沉降量對比如圖7所示,其均方根誤差如表4所示．

圖6 水準點與2018年10月18日至2019年1月22日第三、四采區礦區地面沉降

圖7 沿觀測線沉降結果對比圖

表4 各觀測線上全分辨率與多視法DInSAR監測的地面沉降的均方根誤差

結合圖6,分析圖7和表4可以看出,全分辨率與多視法DInSAR監測的礦區三、四采區地面沉降與水準監測的地面沉降趨勢在空間上基本一致,其精度也非常接近.圖7a,7d中Q和L觀測線上全分辨率DInSAR監測的沉降曲線波動要比多視法DInSAR監測的沉降曲線的波動明顯,但均方根誤差較小,圖7b,7c中G和H觀測線上兩種DInSAR監測的沉降曲線比較平滑,全分辨率DInSAR比多視法DInSAR監測結果的均方根誤差大．

6 討論

全分辨率與多視法DInSAR監測的礦區地面沉降累計結果與濟陽礦區插值后同時間段的精密水準測量所獲取的地面沉降結果在沉降盆地中心位置有明顯偏移,最大沉降量差異非常顯著,這主要由InSAR相位干涉測量和側視成像的固有局限性引起的．

InSAR相位干涉測量的固有局限性主要是指InSAR干涉測量監測形變的能力和精度很大程度上取決于波長、入射角、地面分辨率、重訪周期等SAR成像固有參數.當同一像素單元內發生的形變量超出干涉波長的一半時,干涉測量結果無法反映地表真實形變[21].據此,Massonnet等[22]提出形變梯度的概念,并給出了差分InSAR監測到的最大形變梯度的理論公式,即雷達波長的一半與像素大小的比值.在實際應用中,由于地表形變信息極易被噪聲掩蓋,InSAR監測到的最大形變梯度真實值遠小于理論值[23]．

InSAR側視成像的固有局限性主要是指單軌道(升/降軌)DInSAR技術僅能監測得到監測區域一維地表形變信息,它僅僅代表了傳感器成像方向上的真實空間位移,即雷達視線(Line of Sight,LOS)方向上的投影[24].由于現有在軌SAR 衛星的飛行軌道一般為近極地軌道,且多為右側視成像,這種觀測方式導致了DInSAR技術在垂直方向和水平方向(東西向和南北向)的形變監測敏感度存在差異,一維視線LOS向形變場并不能完全反映監測對象的真實形變狀態,即存在視線向模糊問題[24].因此,直接由一維LOS向上的形變量轉換成常用垂直向的形變量時,就會存在由視線向模糊引起的轉換誤差的問題,當研究區水平形變較大時,由單軌道(升/降軌)DInSAR技術監測的垂直向形變量值與實際地表垂直向的形變量值的差異也就增大,而且會出現形變中心偏移的問題[25]．

7 結論

本研究利用全分辨率和多視法DInSAR分別對山東濟南某礦區5個Sentinel-1A SAR影像干涉對進行處理,結合礦區5個工作面開采信息和4條水準觀測路線的水準監測數據進行精度驗證和對比分析,得到如下主要結論：

2018年10月18日至2019年1月22日,全分辨率與多視法DInSAR均監測到礦區5個工作面上方的地面沉降逐漸加重,兩種方法監測得到的礦區地面沉降的分布、形態和變化趨勢比較一致,與礦區工作面的開采情況相符.但二者監測到的最大沉降值、最大沉降的位置、范圍的大小等并不相同．

理論上來說,全分辨率DInSAR由于沒有進行距離向上的多視濾波處理,其監測大梯度地面沉降的能力應該更強,但實際上,由于大梯度地面沉降區域的相干性較低,全分辨率DInSAR所得到的差分干涉條紋的連續性變差,從而使相位解纏和形變解算的精度受到較大的影響,監測得到的地面沉降最大值與多視法DInSAR監測的地面沉降最大值并沒有較大差異．

在水準觀測線上,全分辨率與多視法DInSAR監測得到的地面沉降分布曲線變化趨勢非常一致,沉降量級相差較小,與水準監測的地面沉降分布曲線基本一致,二者監測的沉降盆地中心位置和最大沉降量也比較一致,但與水準監測的沉降盆地中心位置有明顯偏移,最大沉降量也有明顯差異.全分辨率DInSAR并未明顯提高多視法DInSAR監測大梯度地面沉降的能力．