全分辨率與多視法DInSAR礦區(qū)地面沉降監(jiān)測研究

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

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

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

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

1 研究區(qū)概況

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

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

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

2.1 SAR影像數(shù)據(jù)

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

表1 DInSAR差分干涉對

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

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

2.3 礦區(qū)部分工作面及其水準(zhǔn)數(shù)據(jù)

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

圖2 三、四采區(qū)部分工作面及其水準(zhǔn)點(diǎn)布設(shè)

表2 三、四采區(qū)5個(gè)工作面的基本信息

3 礦區(qū)沉降監(jiān)測原理與數(shù)據(jù)處理流程

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

(1)

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

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

(2)

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

(3)

式中：θ為雷達(dá)視線的入射角．

技術(shù)具體流程如圖3所示．

圖3 數(shù)據(jù)處理流程圖

4 結(jié)果及分析

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

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

分析圖4可以看出：

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

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

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

圖5 全分辨率與多視法DInSAR監(jiān)測的礦區(qū)地面沉降對比

結(jié)合表2,分析圖5可以看出：

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

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

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

為定量對比分析5個(gè)監(jiān)測時(shí)間段內(nèi)全分辨率與多視法DInSAR監(jiān)測的礦區(qū)地面沉降情況,表3中對全分辨率和多視法DInSAR監(jiān)測的各時(shí)段的最大沉降量、最大沉降位置進(jìn)行了統(tǒng)計(jì)．

表3 全分辨率與多視法DInSAR監(jiān)測的最大沉降值及其位置統(tǒng)計(jì)表

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

5 精度驗(yàn)證與對比分析

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

圖6 水準(zhǔn)點(diǎn)與2018年10月18日至2019年1月22日第三、四采區(qū)礦區(qū)地面沉降

圖7 沿觀測線沉降結(jié)果對比圖

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

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

6 討論

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

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

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

7 結(jié)論

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

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

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

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