基于雙軌D-InSAR技術(shù)的礦區(qū)沉降監(jiān)測(cè)

朱文德，王長(zhǎng)栓

（1.廣西壯族自治區(qū)地理信息測(cè)繪院，廣西 柳州 545006）

D-InSAR 是以合成孔徑雷達(dá)復(fù)數(shù)據(jù)提取的相位信息為信息源，利用雷達(dá)圖像的相位差信息來(lái)精確測(cè)量出圖像上每一點(diǎn)的三維位置和提取地面目標(biāo)點(diǎn)微小地形變化的技術(shù)[1]。該技術(shù)不同于水準(zhǔn)、GPS 測(cè)量等傳統(tǒng)的地表沉降監(jiān)測(cè)方法，其具有覆蓋范圍廣、觀測(cè)連續(xù)、高效率、高分辨率、低成本、無(wú)須布設(shè)大量觀測(cè)點(diǎn)等優(yōu)點(diǎn)。在地震位移測(cè)量、監(jiān)測(cè)火山運(yùn)動(dòng)、冰川漂移、地表沉降及山體滑坡等方面展現(xiàn)了其獨(dú)特的優(yōu)越性。1989 年，Gabriel[2]等首次利用D-InSAR 技術(shù)實(shí)現(xiàn)cm 級(jí)的地表形變監(jiān)測(cè)。1993 年，Massonnet[3-4]等在Nature 上發(fā)表了震驚國(guó)際地震界的成果（利用D-InSAR 技術(shù)成功獲取了1992 年Landers 地區(qū)7.3 級(jí)大地震的同震形變場(chǎng)，與GPS 測(cè)量所獲得的結(jié)果吻合度非常高）。二十世紀(jì)九十年代以來(lái)，國(guó)內(nèi)的眾多學(xué)者也在合成孔徑雷達(dá)干涉測(cè)量領(lǐng)域取得顯著的成果。2001 年單新建[5]等基于InSAR 技術(shù)提取了西藏瑪尼地區(qū)的DEM，所獲結(jié)果優(yōu)于1∶20 萬(wàn)DEM；2002 年利用D-InSAR 技術(shù)通過(guò)三通差分干涉處理獲取西藏瑪尼地震同震形變場(chǎng)[6]。2001 年，劉國(guó)祥[7]等利用D-InSAR 技術(shù)獲取了香港赤臘角機(jī)場(chǎng)近乎1 a 內(nèi)的非均勻沉降場(chǎng)，所獲結(jié)果與離散水準(zhǔn)監(jiān)測(cè)結(jié)果相比吻合較好。2004 年，葛大慶[8]等將合成孔徑雷達(dá)干涉測(cè)量技術(shù)應(yīng)用于礦山沉降監(jiān)測(cè)，為煤礦區(qū)地表時(shí)空演變過(guò)程研究和開(kāi)采沉陷實(shí)時(shí)動(dòng)態(tài)監(jiān)測(cè)提供了新的技術(shù)方法。這些研究在很大程度上推動(dòng)了我國(guó)D-InSAR 技術(shù)的發(fā)展。

論文闡述了利用雙軌D-InSAR 技術(shù)進(jìn)行礦區(qū)沉降監(jiān)測(cè)的基本原理和技術(shù)方法，選用覆蓋山東某礦區(qū)2017-02-17、2017-03-11、2017-04-02 三景雷達(dá)SAR 影像，獲取該時(shí)段礦區(qū)的沉降信息，并針對(duì)沉降情況進(jìn)行精細(xì)化分析。

1 D-InSAR技術(shù)原理及方法

1.1 D-InSAR 測(cè)量地表形變基本原理

D-InSAR 技術(shù)是在單天線模式下，進(jìn)行重軌干涉測(cè)量，即在大致相同的軌道上，于2 個(gè)時(shí)刻獲取同一地區(qū)的數(shù)據(jù)。技術(shù)方法包括二軌法（Two-Pass）、三軌法（Three-Pass）和四軌法（Four-Pass）[9]。本文以二軌法為例，簡(jiǎn)要介紹D-InSAR 測(cè)量地表形變的基本原理。

二軌法的基本思想是對(duì)同一地區(qū)變形前、變形后的兩幅SAR 圖像進(jìn)行干涉處理，再引入已知的外部DEM，模擬地形相位，與干涉結(jié)果進(jìn)行差分以去除干涉圖中的地形因素，獲取地表形變信息。該方法所需的SAR 數(shù)據(jù)量小，且通過(guò)引入外部DEM 消除地形相位，減少數(shù)據(jù)處理的工作量。但外部DEM 本身存在的誤差會(huì)傳遞到最終的形變結(jié)果中，此外，利用外部DEM 模擬的SAR 圖像需要與主影像進(jìn)行配準(zhǔn)，配準(zhǔn)的精度也會(huì)影響最終形變監(jiān)測(cè)精度。二軌法監(jiān)測(cè)地表形變的基本原理可以表示為[10]:

式中，φ d表示研究區(qū)域發(fā)生形變前后兩幅SAR 影像生成的干涉圖；φsim,t表示采用外部DEM 模擬生成的地形相位；ΔRtow即為去除地形因素后的形變結(jié)果。二軌法數(shù)據(jù)處理流程如圖1 所示。

圖1 二軌法數(shù)據(jù)處理流程

1.2 D-InSAR 沉降監(jiān)測(cè)數(shù)據(jù)處理流程

1）復(fù)圖像配準(zhǔn)。復(fù)圖像配準(zhǔn)的主要目的是使主、輔影像的像素對(duì)應(yīng)同一地面單元?；静僮靼ㄍc(diǎn)偏移量計(jì)算、偏移量擬合及輔影像重采樣[11]。主要配準(zhǔn)方法有:①相干系數(shù)法，選擇搜索窗口后，逐一計(jì)算每個(gè)像元上的相干系數(shù)，取其最大值對(duì)應(yīng)的點(diǎn)作為配準(zhǔn)點(diǎn)；②平均波動(dòng)函數(shù)法，計(jì)算每個(gè)像元的干涉相位變化梯度，計(jì)算其平均值，值最小的點(diǎn)為配準(zhǔn)點(diǎn)。

2）生成干涉圖及干涉質(zhì)量分析。主、輔影像完成數(shù)據(jù)配準(zhǔn)后，進(jìn)行共軛相乘，生成干涉圖，此時(shí)的相位是纏繞的，數(shù)值范圍為[-π,π]。為確保最終結(jié)果的精度，需要保證干涉的質(zhì)量，通常以相干系數(shù)作為評(píng)定干涉質(zhì)量的指標(biāo):

式中，M×N代表目標(biāo)窗口大??；u1(m,n)和u2(m,n)表示目標(biāo)窗口中(m,n)位置上的2 個(gè)復(fù)型數(shù)據(jù)；u*2(m,n)表示u2(m,n)的共軛復(fù)數(shù)。相干系數(shù)的值越大，代表干涉的質(zhì)量越好。

3）去平地效應(yīng)及相位濾波?；€的存在致使沒(méi)有高程變化的平地產(chǎn)生干涉相位，為了降低解纏難度，將每一像素初始相位中斜距方向上的平地相位剔除。此外，干涉相位會(huì)受到很多噪聲的影響，噪聲的存在會(huì)降低相位解纏的精度，甚至?xí)?duì)最終形變結(jié)果產(chǎn)生影響。因此，在相位解纏前必須進(jìn)行濾波處理。

4）相位解纏。差分干涉圖中的相位被纏繞在[-π,π] 范圍內(nèi)，必須進(jìn)行相位解纏來(lái)獲取反映地面高程的真實(shí)相位。常用的相位解纏方法有3 種:①Goldstein 枝切法；②最小費(fèi)用流法；③加權(quán)解纏法。

5）形變圖生成及地理編碼。通過(guò)公式ΔRtow=λ/(4π×φreal)，將解纏后的真實(shí)相位φreal轉(zhuǎn)換成雷達(dá)視線向上的形變?chǔ)tow。為了更為直觀的分析研究區(qū)域的形變，通常將雷達(dá)視線向上的形變圖進(jìn)行地理編碼，轉(zhuǎn)換到地圖投影坐標(biāo)系下。

2 D-InSAR在礦區(qū)沉降監(jiān)測(cè)中的應(yīng)用實(shí)例

2.1 實(shí)驗(yàn)區(qū)概況

本文選用礦產(chǎn)資源豐富的山東某礦區(qū)作為實(shí)驗(yàn)區(qū)域，位置及范圍如圖2 所示。該礦區(qū)位于山東省西南部，含煤面積約4 826 km2，是全國(guó)重點(diǎn)開(kāi)發(fā)的煤炭基地之一。煤炭開(kāi)采會(huì)導(dǎo)致礦區(qū)地表出現(xiàn)沉降、房屋坍塌、道路裂縫等問(wèn)題，為研究礦區(qū)的沉降情況，本文選取覆蓋研究區(qū)域的三景SAR 影像，基本參數(shù)如表1 所示。

圖2 研究區(qū)域示意圖

表1 干涉像對(duì)參數(shù)情況列表

2.2 數(shù)據(jù)處理關(guān)鍵步驟及成果

引入外部90 m 精度的DEM 數(shù)據(jù)，采用二軌法對(duì)上述3 個(gè)相對(duì)進(jìn)行數(shù)據(jù)配準(zhǔn)及差分處理，得到差分干涉圖（如圖3 所示）。

由圖3 可以看出干涉相對(duì)1 所得的干涉效果較清晰，沉降區(qū)域較明顯；干涉相對(duì)2、3 存在明顯的噪聲，影響干涉結(jié)果。為過(guò)濾噪聲，提取精確的相位信息，采用Goldstein 濾波對(duì)上述差分干涉圖進(jìn)行濾波處理，濾波后的干涉圖如圖4 所示。

圖3 差分干涉圖

濾波后的干涉圖，去除了噪聲的影響，圖像更加清晰，標(biāo)記位置的干涉條紋發(fā)生突變，條紋密集，色度跨越大，可見(jiàn)此區(qū)域由于煤礦長(zhǎng)期開(kāi)采，發(fā)生了地表沉降，為了獲取真實(shí)相位值，進(jìn)一步做相位解纏處理。利用式ΔRtow=λ/(4π×φreal)，將解纏后的真實(shí)相位φreal轉(zhuǎn)換成雷達(dá)視線向上的形變?chǔ)tow，再通過(guò)地理編碼轉(zhuǎn)換為地圖投影坐標(biāo)系下，圖5 為所得形變圖。

需要說(shuō)明的是圖5 中的形變與圖4 表示的區(qū)域是一致的，由于所采用的SAR 數(shù)據(jù)為升軌數(shù)據(jù)，所以未經(jīng)地理編碼的干涉圖的經(jīng)緯度與實(shí)際地理經(jīng)緯度相反，經(jīng)過(guò)地理編碼后的形變圖所表示的經(jīng)緯度與實(shí)際地理經(jīng)緯度一致。

圖4 濾波后的干涉圖

圖5 形變結(jié)果

3 礦區(qū)沉降結(jié)果分析

由圖5 可以得知，2017-02-17～2017-04-02 期間，研究區(qū)域內(nèi)出現(xiàn)A、B、C、D 四處明顯的沉降盆地，是由于煤礦的長(zhǎng)期開(kāi)采，導(dǎo)致地面出現(xiàn)沉降，且對(duì)周邊地區(qū)造成一定程度的影響。其中D 處的煤礦沉降尤為嚴(yán)重，為更直觀的分析礦區(qū)沉降情況，以該煤礦為例，作圖示剖面線，進(jìn)行精細(xì)化分析，如圖6 所示。

圖6 剖面線沉降圖

通過(guò)對(duì)比2017-02-17 ～2017-03-11、2017-03-11～2017-04-02、2017-02-17～2017-04-02 三期成像的剖面圖可以得到:D 煤礦由于長(zhǎng)期開(kāi)采及工作面間的相互影響，形成2 個(gè)沉降漏斗。3 個(gè)干涉相對(duì)均呈現(xiàn)這一特征，且沉降位置及整體沉降趨勢(shì)一致。2017-02-17 ～2017-03-11 期間2 個(gè)沉降漏斗最大沉降量分別達(dá)到3.8 cm、4.7 cm，2017-03-11 ～2017-04-02 期間2 個(gè)沉降漏斗最大沉降量分別達(dá)到3.0 cm、4.8 cm，2017-02-17 ～2017-04-02 期間2 個(gè)沉降漏斗最大沉降量分別為3.0 cm、4.5 cm。干涉相對(duì)3 的沉降結(jié)果應(yīng)為干涉相對(duì)1、2 之和，但由于時(shí)間基線過(guò)長(zhǎng)，加上各種誤差的影響，該時(shí)段失相干嚴(yán)重，沉降漏斗周邊出現(xiàn)大量不連續(xù)的點(diǎn)，造成3 個(gè)成像期間沉降量有所差異，說(shuō)明時(shí)間基線過(guò)長(zhǎng)會(huì)降低雙軌D-InSAR 的監(jiān)測(cè)精度。

4 結(jié) 語(yǔ)

本文論述了基于雙軌道D-InSAR 技術(shù)進(jìn)行礦區(qū)沉降監(jiān)測(cè)的基本原理和技術(shù)方法，采用2017-02-17、2017-03-11、2017-04-02 三景SAR 影像提取該時(shí)段礦區(qū)的沉降信息，并對(duì)沉降情況作進(jìn)一步的精細(xì)化分析。結(jié)果表明，在2017-02-17、2017-03-11、2017-04-02成像期間，出現(xiàn)4 個(gè)明顯的沉降區(qū)域，其中D 處煤礦沉降最為嚴(yán)重，形成2 個(gè)連續(xù)的沉降漏斗，并呈現(xiàn)沉降區(qū)域擴(kuò)大，沉降情況加重的趨勢(shì)。實(shí)驗(yàn)證明D-InSAR 技術(shù)可以反映可靠的礦區(qū)沉降情況，能夠用于礦區(qū)沉降的實(shí)時(shí)監(jiān)測(cè)，但由于外界條件制約，各種系統(tǒng)誤差的影響，還需通過(guò)限制時(shí)間基線，結(jié)合水準(zhǔn)數(shù)據(jù)或GPS 等測(cè)量結(jié)果，以提高成果的精度。