余景波,劉國(guó)林,曹振坦

(山東科技大學(xué)測(cè)繪科學(xué)與工程學(xué)院,山東青島266510)

0 引 言

采礦引起的地面沉降和塌陷是煤礦開(kāi)采地區(qū)經(jīng)常發(fā)生的一種破壞性災(zāi)難,這不僅僅是個(gè)環(huán)境問(wèn)題,而且會(huì)影響到煤礦附近地區(qū)社會(huì)和諧與可持續(xù)發(fā)展[1]。因此,有必要利用先進(jìn)的技術(shù)進(jìn)行監(jiān)測(cè)和控制礦區(qū)地表沉降引起的破壞,從而保證礦區(qū)及礦區(qū)附近地區(qū)的安全和穩(wěn)定發(fā)展。以InSAR技術(shù)為基礎(chǔ)進(jìn)行地面沉降監(jiān)測(cè),是國(guó)內(nèi)外研究的熱點(diǎn)之一。Gabriel等[2]首次論證了差分干涉測(cè)量技術(shù)可用于探測(cè)厘米級(jí)的地表形變。Wegm等[3]利用DInSAR對(duì)德國(guó) Ruhrgebiet地區(qū)地面沉降進(jìn)行研究,獲取了因采礦而導(dǎo)致的地表沉降信息。Ge等[4]將GPS和InSAR技術(shù)結(jié)合在一起進(jìn)行礦區(qū)地表形變測(cè)量,準(zhǔn)確獲取了澳大利亞Appin地區(qū)煤礦開(kāi)采造成的地面沉降信息。吳立新等[5]利用多時(shí)相D-InSAR結(jié)合角反射器方法進(jìn)行了工礦區(qū)地表沉降監(jiān)測(cè)。

采用雙軌法和三軌法兩種差分干涉測(cè)量方法對(duì)6景覆蓋濟(jì)寧某礦區(qū)的ENVISAT ASAR數(shù)據(jù)進(jìn)行差分干涉處理,并對(duì)實(shí)驗(yàn)處理結(jié)果進(jìn)行討論和分析。

1 基本原理

1.1 InSAR及差分干涉測(cè)量(D-InSAR)的基本原理

InSAR技術(shù)是以合成孔徑雷達(dá)復(fù)數(shù)據(jù)提取的相位信息為信息源獲取地表三維信息和變化信息的一項(xiàng)技術(shù),其基本思想是:利用兩副天線同時(shí)成像或一副天線相隔一定時(shí)間重復(fù)成像,獲取同一區(qū)域的復(fù)雷達(dá)圖像對(duì),由于兩副天線與地面某一目標(biāo)之間的距離不等,形成干涉圖,干涉圖中的相位值即為兩次成像的相位差測(cè)量值,根據(jù)兩次成像的相位差與地面目標(biāo)的三維空間位置之間存在的幾何關(guān)系,利用飛行軌道的參數(shù),即可測(cè)定地面目標(biāo)的三維坐標(biāo),它可以用來(lái)提供大范圍的高精度數(shù)字高程模型(DEM),并用于探測(cè)地表形變[6]。

差分干涉測(cè)量(D-InSAR)技術(shù)是以雷達(dá)干涉測(cè)量為出發(fā)點(diǎn),通過(guò)星載或機(jī)載的雷達(dá)天線(單天線或雙天線)主動(dòng)向地面目標(biāo)發(fā)射電磁波,并接收地面目標(biāo)的后向散射回波,同時(shí)記錄各自回波的相位信息。由于在所獲取的相位信息中,含有大氣延遲、軌道誤差、平地效應(yīng)、地形起伏、目標(biāo)兩次成像過(guò)程中的微小形變及噪聲信息,因此,D-InSAR的處理過(guò)程就是通過(guò)外部DEM(“雙軌法”)或其它干涉影像對(duì)獲取的相位(“三軌法”或“四軌法”)進(jìn)行差分處理,以便去除無(wú)用相位信息而獲取形變相位的一個(gè)過(guò)程[7]。

1.2 雙軌法和三軌法數(shù)據(jù)處理流程

差分干涉測(cè)量數(shù)據(jù)處理,經(jīng)過(guò)配準(zhǔn)、重采樣、濾波、平地效應(yīng)消除、干涉條紋濾波和相位解纏,相高轉(zhuǎn)換,得到差分干涉圖、增強(qiáng)干涉圖、相干圖、以及地面形變圖,通過(guò)進(jìn)一步后處理得到地面沉降分布圖。InSAR數(shù)據(jù)處理的基本流程如圖1所示[8]:

圖1 InSAR數(shù)據(jù)處理基本流程

雙軌法是Massonnet等[9]提出的,該方法是利用監(jiān)測(cè)區(qū)域地表變化前后兩幅影像生成干涉紋圖,再利用事先獲取的DEM數(shù)據(jù)模擬紋圖,從干涉紋圖中去除地形信息以得到地表變化信息。無(wú)需進(jìn)行相位解纏,處理工作量少是該方法的優(yōu)點(diǎn)。但是,對(duì)于無(wú)DEM 數(shù)據(jù)的地區(qū),無(wú)法使用該方法,同時(shí)引入DEM數(shù)據(jù)的同時(shí)有可能帶入新誤差,這些是該方法的缺點(diǎn)。雙軌法數(shù)據(jù)處理基本流程如圖2所示:

圖2 雙軌法差分干涉數(shù)據(jù)處理流程

三軌法是Zebker等[10]人提出的一種數(shù)據(jù)處理方法,該方法利用三幅影像生成兩幅干涉紋圖,一副作為地形干涉對(duì)反映地形信息,另一幅作為形變干涉對(duì)反映地表形變信息,進(jìn)行平地效應(yīng)消除后,分別進(jìn)行相位解纏,最后利用差分干涉原理計(jì)算得到地形信息。三軌法的主要優(yōu)點(diǎn)是無(wú)需地面信息、數(shù)據(jù)間的配準(zhǔn)比較容易實(shí)現(xiàn);其主要缺點(diǎn)是相位解纏質(zhì)量的好壞將對(duì)最終結(jié)果產(chǎn)生影響,但在相干性較高地區(qū)可以獲得較好的結(jié)果。三軌法數(shù)據(jù)處理的基本流程如圖3所示。

圖3 三軌法差分?jǐn)?shù)據(jù)處理簡(jiǎn)單流程圖

2 實(shí)驗(yàn)研究數(shù)據(jù)選取

研究的區(qū)域位于濟(jì)寧市的某煤礦,該煤礦位于濟(jì)寧煤田西北部,礦井范圍南北長(zhǎng)4～4.5 km,東西寬4～6 km。礦區(qū)內(nèi)地形平坦,地勢(shì)東北略高,西南稍低,自然坡度為萬(wàn)分之七。氣候溫和,四季分明,交通方便。

利用ENVISAT衛(wèi)星獲取的6景ASAR數(shù)據(jù)進(jìn)行差分干涉測(cè)量處理,ASAR數(shù)據(jù)參數(shù)見(jiàn)表1所示。

表1 所選用的 ASAR數(shù)據(jù)列表

根據(jù)研究的需要,這6景ASAR數(shù)據(jù)可以組成不同的干涉像對(duì),所需的干涉像對(duì)垂直基線和時(shí)間間隔參數(shù)情況見(jiàn)表2和表3。

對(duì)表2中的干涉像對(duì)進(jìn)行雙軌法差分干涉處理,對(duì)表3的干涉像對(duì)進(jìn)行三軌法差分干涉處理。

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

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

雙軌法利用外部DEM數(shù)據(jù)來(lái)消除干涉相位圖中的地形因素影響,在實(shí)驗(yàn)中采用分辨率為30 m Aster數(shù)據(jù)作為雙軌法差分干涉處理數(shù)據(jù)的外部DEM。

3 差分干涉處理和結(jié)果分析

采用雙軌法處理表1中的干涉像對(duì),把2009年4月3日和2009年6月12日獲取的影像為主圖像,2009年5月8日和2009年7月17日獲取的影像為輔圖像,同時(shí)引入分辨率為30 m的Aster數(shù)據(jù)作為外部DEM,進(jìn)行雙軌法差分干涉處理,處理后分別生成差分干涉圖、增強(qiáng)干涉圖、相干圖、地面形變圖,如圖4、5所示。

在圖4和圖5的地面形變圖中,可以看到盡管所選取的實(shí)驗(yàn)數(shù)據(jù)時(shí)間跨度不同,但是差分干涉處理生成地面形變圖中沉降分布范圍基本上大體一致,但是由于干涉選用不同時(shí)期的數(shù)據(jù)進(jìn)行差分處理,數(shù)據(jù)相干性存在差異,得出結(jié)果與實(shí)際情況會(huì)有誤差。這說(shuō)明雙軌法可以從整體上確定礦區(qū)的沉降分布情況。

在圖4和圖5的差分干涉圖中,可以看到在圖4的差分干涉圖中沒(méi)有明顯的條紋,在圖5差分干涉圖的中央位置有一明顯的條紋,在其他地方都沒(méi)有表現(xiàn)出來(lái)。從表2中,可以看出生成圖4結(jié)果的干涉像對(duì)垂直基線較長(zhǎng),基線長(zhǎng)度是相干一個(gè)重要的限制因素,垂直基線過(guò)長(zhǎng)可能導(dǎo)致失相干。造成圖4中差分干涉圖沒(méi)明顯干涉條紋的一個(gè)重要原因可能是該干涉像對(duì)垂直基線過(guò)長(zhǎng),造成該干涉像對(duì)相干性不強(qiáng)。圖4的增強(qiáng)干涉圖,相干效果不好,表現(xiàn)出雜亂無(wú)章的彩色斑點(diǎn);圖5增強(qiáng)干涉圖中,有一相干性較好的地方,可能這個(gè)地方是該時(shí)間段,地面沉降最大的所在地,在其它地方基本上沒(méi)有干涉出來(lái),只是一些噪聲表現(xiàn)出的雜亂無(wú)章的彩色斑點(diǎn)。上面是對(duì)雙軌法差分干涉處理數(shù)據(jù)生成干涉圖的分析。由于雙軌法需要引入外部DEM,外部數(shù)據(jù)帶來(lái)誤差可能也是導(dǎo)致上面干涉圖現(xiàn)象的其中一個(gè)原因。

將2008年12月19日獲取的影像作為主圖像,分別以2009年 2月 27日和 2009年5月 8日獲取的影像為輔影像,其中,把2008年12月19日和2009年2月27日獲取的影像作為地表形變前獲取的兩幅影像,作為地形干涉對(duì)進(jìn)行處理;同樣,把2009年5月8日獲取的影像作為地表形變后獲取的影像,將其與主圖像作為形變干涉對(duì)進(jìn)行處理;同樣,以2009年5月8日獲取的影像作為主影像,將2009年5月8日和 2009年6月12日獲取的影像作為地形干涉對(duì)進(jìn)行處理;把2009年5月8日和2009年7月17日獲取的影像作為形變干涉對(duì)進(jìn)行處理;最后,進(jìn)行三軌法差分干涉處理,獲取差分干涉圖、增強(qiáng)干涉圖相干圖和地面形變圖,如圖6、7所示。

在圖6和圖7中的地面形變圖中,也可看到沉降分布的相同之處,在形變圖中形成了一些“漏斗區(qū)”。這可以說(shuō)明,三軌法可以檢測(cè)到區(qū)域地表形變發(fā)生的大致位置和范圍,并且從得到的結(jié)果,基本上可以確定地表形變大致位置和范圍。這表明三軌法可以整體上確定礦區(qū)沉降分布情況。

在圖6和圖7中,圖6中差分干涉圖比圖7中差分干涉圖干涉條紋多些,除了垂直基線長(zhǎng)度的原因,季節(jié)可能也是一個(gè)原因,即時(shí)間去相干。生成圖6結(jié)果的干涉像對(duì)所處的時(shí)間段正是研究區(qū)域植被稀疏、雨水較少、天氣干燥的時(shí)期,這有助于提高獲取的ASAR數(shù)據(jù)差分干涉處理精度,提高干涉效果;而生成圖7結(jié)果的干涉像對(duì),正處于研究區(qū)的夏季,植被茂盛、雨水較多,這可能對(duì)干涉效果有影響;再者,生成圖6和圖7結(jié)果,需要先進(jìn)行地形干涉像對(duì)差分干涉處理,從表3中,可以看出兩個(gè)作為地形干涉像對(duì)的時(shí)間間隔相差不是很大,前者干涉像對(duì)時(shí)間處于研究區(qū)冬季,后則處于夏季;并且,由于季節(jié)不同,生成含地形信息干涉圖效果不同,這就造成了最終結(jié)果的差異;再者,作為地形干涉像對(duì)的垂直基線長(zhǎng)度對(duì)DEM的生成具有影響,DEM的質(zhì)量又會(huì)影響到差分干涉處理的結(jié)果。在圖6和圖7的增強(qiáng)干涉圖中,干涉條紋明顯程度也不一樣,前者較多、較明顯;后者則較差一些。除了上面的分析外,相位解纏的質(zhì)量也有差別,這也是造成上面現(xiàn)象出現(xiàn)的一個(gè)原因。以上是對(duì)三軌法數(shù)據(jù)處理生成干涉圖的分析。

從差分干涉圖到增強(qiáng)干涉圖有時(shí)需要除去殘余地形相位,主要有梯度擬合法和基線向量調(diào)整法進(jìn)行殘余地形相位去除兩種方法。在數(shù)據(jù)處理中因干涉圖條紋不明顯,可以看作沒(méi)有殘余地形相位,沒(méi)有進(jìn)行殘余地形相位去除操作,因而實(shí)際存在的殘余地形相位會(huì)對(duì)增強(qiáng)干涉圖質(zhì)量帶來(lái)影響,從而影響最終差分干涉處理的結(jié)果。

相干圖表示的是兩幅影像的相關(guān)性,在相干性低的地方,可以通過(guò)相干圖的對(duì)比檢測(cè)出礦區(qū)地表沉降的大致位置和范圍。垂直基線長(zhǎng)度,對(duì)相干圖質(zhì)量影響很大。從圖4～圖7的相干圖中可以看出圖4相干圖質(zhì)量很差,其他三幅相干圖質(zhì)量較好,能夠看出明顯地物的特征。

對(duì)圖4～圖7中地面形變圖進(jìn)一步處理,可以獲取該研究區(qū)域在各時(shí)間段的沉降分布,從而可以從宏觀尺度分析該區(qū)域在不同時(shí)間段的沉降量和沉降分布。

從圖8～圖9中,可以看出選取的研究區(qū)域都存在著不同程度的沉降,一些不明顯的沉降可能是開(kāi)采地下水等原因引起的,也可能是數(shù)據(jù)處理過(guò)程不可避免的誤差引起的。在數(shù)據(jù)處理中,選用的差分干涉處理方法不同,最后產(chǎn)生的結(jié)果會(huì)有差別。這說(shuō)明差分干涉測(cè)量技術(shù)可以得到礦區(qū)的沉降分布。

以雙軌法處理的2009年6月12日-2009年7月17日期間獲取的干涉像對(duì)和三軌法處理2009年5月8日-2009年6月12日-2009年7月17日期間或取得干涉像對(duì)為例,進(jìn)行沉降面積估算。這兩個(gè)時(shí)間段的沉降面積估算如表4、5所示。

表4 2009年6月12日-2009年7月17日沉降面積估算(雙軌法)

表5 2008年12月19日-2009年2月27日-2009年5月8日沉降面積估算(三軌法)

從表4和表5,可以看出雙軌差分干涉測(cè)量和三軌差分干涉測(cè)量都可對(duì)不同沉降值的沉降面積進(jìn)行估算。這說(shuō)明差分干涉測(cè)量技術(shù)可以對(duì)沉降面積進(jìn)行估算。

4 結(jié) 論

隨著社會(huì)經(jīng)濟(jì)的快速發(fā)展,我國(guó)對(duì)煤炭資源的開(kāi)采需求將進(jìn)一步擴(kuò)大,開(kāi)采深度和范圍也將進(jìn)一步拓展,勢(shì)必會(huì)造成礦區(qū)地表沉降災(zāi)害的發(fā)生,因此必須對(duì)礦區(qū)地表沉陷進(jìn)行監(jiān)測(cè)。差分干涉測(cè)量技術(shù)具有其自身的優(yōu)勢(shì),越來(lái)越受到人們的重視。

采用雙軌法和三軌法差分干涉處理方法處理了6景覆蓋濟(jì)寧某礦區(qū)ENVISAT ASAR數(shù)據(jù),并對(duì)數(shù)據(jù)處理結(jié)果進(jìn)行了分析和討論,從而證明了差分干涉測(cè)量(D-InSAR)技術(shù)在礦區(qū)地面沉降監(jiān)測(cè)應(yīng)用中的可行性。但是,應(yīng)用差分干涉測(cè)量技術(shù)獲取礦區(qū)地表形變信息受多種因素影響,因而需要進(jìn)一步將其測(cè)量結(jié)果與水準(zhǔn)測(cè)量或GPS測(cè)量結(jié)果進(jìn)行比較,以評(píng)價(jià)其測(cè)量精度,這是我們今后進(jìn)一步的研究工作。

[1]侯建國(guó),楊成生,王亞男,王宏宇.煤礦開(kāi)采沉降監(jiān)測(cè)-基于合成孔徑雷達(dá)差分干涉測(cè)量的研究[J].自然災(zāi)害學(xué)報(bào),2009,18(6):71-72.

[2]Gabriel A K,Goldstein R M,Zebker H A.Mapping Small Elevation Changes Over Large Areas:Differential Radar Interferometry[J].Journal of Geophysiacal Research,1989,94(B7):9183-9191.

[3]Wegm U,Strozzi T,Werner C,et al.Monitoring of mining-Induced Surface Deformation in the Ruhrgebiet(Germany)with SAR Interferometry[C]∥IGARSS'00,Honolulu,USA,2000.

[4]Ge L L,Chen H Y,Han S,et al.Intergrated GPS and Interferometric SAR Techniques for Highly Dense Crustal Deformation Monitoring[C]∥The 14th Internal Technology meeting of the Satellite Division of the U.S.Inst.Of Navigation,Salt Lake City,Utah,2001.

[5]吳立新,高均海,葛大慶,等.基于D-InSAR的煤礦區(qū)開(kāi)采沉降遙感監(jiān)測(cè)技術(shù)分析[J].地理學(xué)與國(guó)土研究,2004,20(2):22-25.

[6]陶秋香,劉國(guó)林,馬飛虎,等.軌道差分干涉測(cè)量技術(shù)在地面沉降監(jiān)測(cè)中的應(yīng)用[J].測(cè)量信息與工程,2008,33(4):1-2.

[7]李晶晶,郭增長(zhǎng).基于D-InSAR技術(shù)的煤礦區(qū)開(kāi)采沉降監(jiān)測(cè)[J].河南理工大學(xué)學(xué)報(bào),2006,25(4):306-309.

[8]李平湘,楊 杰.雷達(dá)干涉測(cè)量原理與應(yīng)用[M].北京:測(cè)繪出版社,2006(44):107-108.

[9]Assonnnet D M.The displacement field of the Lander earth-quake mapped by radar interferometry[J].Nature,1993(364):138-142.

[10]Zerker H A,Rosen P A,Goldstein R M,etal.On the derivation of coseismic displacement fields using Differential Radar Interferometry of the Landers[J].J Geophys Res,1994,99(B10):19617-19634.