基于合成孔徑雷達(dá)差分干涉測(cè)量技術(shù)的礦區(qū)地面形變分析①

張文展

(新鄉(xiāng)航空工業(yè)(集團(tuán))有限公司,河南 新鄉(xiāng)453049)

0 引 言

20世紀(jì)60年代末,合成孔徑雷達(dá)干涉測(cè)量技術(shù)(InSAR)作為合成孔徑雷達(dá)(SAR)技術(shù)的一個(gè)新的應(yīng)用領(lǐng)域發(fā)展起來(lái),尢其是利用合成孔徑雷達(dá)的相位信息提取地表的三維信息和高程變化信息的一項(xiàng)技術(shù)。由其拓展的合成孔徑雷達(dá)差分干涉測(cè)量技術(shù)(D-InSAR)具有全天候、全天時(shí)、覆蓋范圍廣等特點(diǎn),被廣泛應(yīng)用于火山、地震形變、滑坡等監(jiān)測(cè)[1-4],但是,將該技術(shù)應(yīng)用于煤礦開(kāi)采地面形變監(jiān)測(cè)還處于初步研究階段[5],對(duì)其進(jìn)行研究具有重要的理論和實(shí)際意義。利用雙軌法D-InSAR處理ENVISAT衛(wèi)星獲取4景覆蓋濟(jì)寧某礦區(qū)的先進(jìn)的合成孔徑雷達(dá)(ASAR)數(shù)據(jù),分別獲取了濟(jì)寧某礦區(qū)的地面形變相位圖、差分干涉圖、沉降分布以及沉降面積等,并對(duì)該礦區(qū)地面形變進(jìn)行了研究分析。

1 D-InSAR測(cè)量基本原理

圖1是差分干涉測(cè)量的成像幾何示意圖[6]。A1和A2是衛(wèi)星兩次同一地區(qū)成像的位(即天線的位置),R1和R2表示目標(biāo)點(diǎn)P到天線A1和A2距離,其中,R2=R1+ΔR.天線A1和A2對(duì)目標(biāo)點(diǎn)P的測(cè)量相位差為[7]

式中,λ為雷達(dá)信號(hào)的波長(zhǎng)。

由圖1中的幾何關(guān)系模型和余弦定理可得:

圖1 差分干涉測(cè)量成像幾何示意圖

式中:B為天線A1和 A2之間的距離,稱為基線距;θ為雷達(dá)入射角;α為水平線與基線的夾角。通常 ΔR?R1,則式(2)可化簡(jiǎn)為

式中,B∥是基線距沿雷達(dá)視線方向的分量。把式(3)代入式(1)可得:

在地表發(fā)生形變后,進(jìn)行觀測(cè)獲取的相位差不但與地形有關(guān),而且與沿雷達(dá)視線方向的形變量ΔRd有關(guān)。此時(shí)對(duì)目標(biāo)點(diǎn)P的測(cè)量相位差為

式中,B′∥表示地形未發(fā)生形變時(shí),獲取的干涉圖的基線距在雷達(dá)視線方向上的分量。

由式(4)和式(5),可得雷達(dá)視線方向變量Δ Rd所引起的相位為

式中左邊各量可由干涉紋圖的相位和軌道參數(shù)計(jì)算得到,進(jìn)而可確定圖像上每點(diǎn)的視線向形變量ΔRd.

2 數(shù)據(jù)源以及數(shù)據(jù)處理的流程

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

收集了4景覆蓋濟(jì)寧某礦區(qū)的ENVISAT ASAR數(shù)據(jù),這4景ASAR數(shù)據(jù)組合的干涉像對(duì)的具體參數(shù)如表1所示。本文所研究煤礦,礦井范圍南北長(zhǎng)4～4.5 km,東西寬4～6 km。礦區(qū)內(nèi)地形平坦,地勢(shì)東北略高,西南稍低,地面標(biāo)高為+37.04～+41.28 m,平均為+38 m,自然地形坡度為萬(wàn)分之七。氣候溫和,四季分明,氣候宜人;交通方便。ENVISAT衛(wèi)星[8]是由歐洲空間局在2002年3月1日通過(guò)阿里亞納5號(hào)火箭發(fā)射的一顆先進(jìn)的太陽(yáng)同步軌道地球環(huán)境監(jiān)測(cè)衛(wèi)星,該衛(wèi)星上配有的ASAR,是目前先進(jìn)的合成孔徑雷達(dá),具有多極化、多入射角、大寬度等特性。

表1 ASAR像對(duì)的基線列表

主圖像選取考慮到時(shí)間基線和空間基線對(duì)圖像相干性的影響以及成像期間植被的影響以及形成結(jié)果的可靠性等,利用D-InSAR技術(shù)中的雙軌差分干涉法對(duì)表1中的影像數(shù)據(jù)對(duì)進(jìn)行處理。采用雙軌差分的優(yōu)點(diǎn)是外部數(shù)字高程模型(DEM)和滿足干涉的影像對(duì)容易獲取,缺點(diǎn)是如果DEM精度不高,會(huì)對(duì)數(shù)據(jù)處理結(jié)果造成嚴(yán)重影響。本文選用分辨率為90 m的雷達(dá)地形測(cè)繪(SRTM)[9-10]數(shù)據(jù)作為外部DEM.

2.2 數(shù)據(jù)處理流程

采用雙軌差分干涉測(cè)量法處理數(shù)據(jù),該方法利用已知DEM反演的干涉相位進(jìn)行二次差分處理,即從干涉相位中減去地形相位,得到只表征地表形變信息的相位,從而可計(jì)算地表形變量。其具體數(shù)據(jù)處理流程如圖2所示[11]。

圖2 雙軌法差分干涉測(cè)量處理流程

3 D-InSAR結(jié)果及其礦區(qū)地面形變分析

3.1 D-InSAR處理結(jié)果

1)差分干涉圖

干涉圖是像對(duì)的對(duì)應(yīng)點(diǎn)復(fù)共軛相乘得到的,其相位等于原始兩SAR圖像的相位差,反映了不同觀測(cè)位置和天線到對(duì)應(yīng)的分辨單元中心的路程差的函數(shù)。在干涉圖中,2π的相位變化(圖中表現(xiàn)為一個(gè)紅藍(lán)綠顏色周期)稱為一個(gè)干涉條紋,由于復(fù)數(shù)對(duì)其相位的周期性,干涉得出的不是直接兩次測(cè)量的相位差原值,而是其被周期折疊后的主值,即真實(shí)相位 2π的模。研究數(shù)據(jù)通過(guò)雙軌法 D-In-SAR處理后獲得增強(qiáng)差分干涉圖和相干圖如圖3所示。

從圖3中可以看出,2004年 12月26日和2005年1月30日兩幅影像之間的相干性很差,只有極小的區(qū)域相干性可以,在差分干涉圖中,除去地形信息,也難以明確的看出干涉條紋。但2005年1月30日和2005年3月6日,兩幅影像之間的相干性較好,在差分干涉圖中可以可出明顯的沉降變化區(qū)域。

2)地面沉降圖

對(duì)差分干涉圖進(jìn)行相位解纏,并進(jìn)行相高轉(zhuǎn)換后,即可得到沿雷達(dá)視線方向的形變圖,即為沉降圖。

在圖4的沉降圖可以看到所研究區(qū)域的某些區(qū)域形成了明顯的“漏斗”沉降;某些區(qū)域發(fā)生了一些微小的沉降可能是由于開(kāi)采地下水或數(shù)據(jù)處理中不可避免的誤差造成;有些區(qū)域發(fā)生了上升現(xiàn)象,這可能與開(kāi)采礦物有關(guān),也可能與數(shù)據(jù)處理中有關(guān)參數(shù)設(shè)置不當(dāng)有關(guān)。

3.2 礦區(qū)地面形變分析

根據(jù)礦區(qū)范圍,從整體上分析該研究區(qū)域不同時(shí)間段的沉降量、沉降分布及其沉降面積統(tǒng)計(jì)。

1)2004年12月26日-2005年1月30日期間沉降分布及其沉降面積估算

2)2005年1月30日-2005年3月6日期間沉降分布及其沉降面積估算

從圖5和圖7中,可以看出整個(gè)礦區(qū)有不同沉降值的地面面積沉降,在這些地面沉降中,一些沉降程度不太明顯,這些沉降可能是由于地下水開(kāi)采過(guò)度所引起的,也可能是數(shù)據(jù)處理中的誤差所造成的。圖5和圖7的左上方,都存在著由于煤礦開(kāi)采而造成的明顯沉降,分布區(qū)域也大體相同。從圖6和圖8中,可以看出地面形變分布以及沉降程度有所不同,不同沉降值所對(duì)應(yīng)的沉降面積有所不同。

圖4 雙軌法D-InSAR沉降圖

圖8 2005年1月30日-2005年3月6日沉降面積曲線圖

3.3 精度分析

為了驗(yàn)證D-InSAR監(jiān)測(cè)礦區(qū)地面沉降的精度,收集了水準(zhǔn)監(jiān)測(cè)手段獲取2004.11.20-2005.12.31時(shí)間段內(nèi)所研究礦區(qū)地面沉降監(jiān)測(cè)值,與DInSAR監(jiān)測(cè)結(jié)果相比較。表2和圖9給出了DInSAR監(jiān)測(cè)結(jié)果與水準(zhǔn)結(jié)果的比較情況。

圖9 水準(zhǔn)結(jié)果與D-InSAR監(jiān)測(cè)結(jié)果比較

從表2和圖9中,可以看出水準(zhǔn)監(jiān)測(cè)的沉降量結(jié)果和D-InSAR監(jiān)測(cè)的沉降結(jié)果存在著一定差異,兩者之間的相符程度不太高。這可能是由于造成這種現(xiàn)象的原因有:①成像分辨率和波長(zhǎng)的制約;②空間上不一致;③時(shí)間上不一致;④ASAR數(shù)據(jù)的質(zhì)量和數(shù)量:ASAR成像波段容易受到植被影響,所研究礦區(qū)植被分布比較多;缺少2004.11時(shí)間段的ASAR數(shù)據(jù),造成了時(shí)間“漏洞”;⑤雙軌法D-InSAR數(shù)據(jù)處理中,引入外部DEM精度不高。

4 結(jié) 論

采用雙軌法D-InSAR處理覆蓋濟(jì)北某礦區(qū)的ENVISATASAR數(shù)據(jù),提取了增強(qiáng)差分干涉圖、相干圖、沉降圖,并從礦區(qū)地面沉降分析和精度分析等方面對(duì)D-InSAR監(jiān)測(cè)礦區(qū)地面沉降結(jié)果進(jìn)行了分析討論,可以得出:

1)D-InSAR技術(shù)可以從整體上得到礦區(qū)沉降的整體分布、不同時(shí)間段的沉降量以及估算不同沉降值的沉降面積。

2)由于受到ASAR分辨率和成像波長(zhǎng)等因素制約,造成D-InSAR技術(shù)監(jiān)測(cè)礦區(qū)地面沉降結(jié)果與水準(zhǔn)監(jiān)測(cè)沉降結(jié)果之間存在一定差異。

[1] 路 旭.用InSA R作地面沉降監(jiān)測(cè)的實(shí)驗(yàn)研究[J].大地測(cè)量與地球動(dòng)力學(xué),2002,22(4):66-70.

[2] Chris R.Tectonics of the Italian earthquake[EB/OL].[2009-04-06].http://scienceblogs.com./highlyalloch tho nous/2009/04/tectonics of the italicmearth.php.

[3] 龍四春,李 陶.D-InSAR中參考DEM誤差與軌道誤差對(duì)相位貢獻(xiàn)的靈敏度研究[J].遙感信息,2009(1):3-6.

[4] Didier M.The displacement field of the Landers earthquake mapped by radar interferometry[J].Nature 1993(364):138-142.

[5] Howard.On the derivation of coseismic displacement fields usingdifferential radarinterferometry:the Lander earthquake[J].Journal of Geophysical Research Solid Earth,1994,99(B10):617-634.

[6] 鄒積亭,劉遠(yuǎn)明.D-InSAR在地面形變監(jiān)測(cè)中的應(yīng)用研究[J].北京測(cè)繪,2009(2):19-20.

[7] 王 超,張 紅,劉 智.星載合成孔徑雷達(dá)干涉測(cè)量[M].北京:科學(xué)出版社,2002:59-60.

[8] Ramon F H.Radar interferometry data interpretation and error analysis[M].England:Kluwer Academic Publishers,2002:315-340.

[9] 郭華東,王長(zhǎng)林.全天候全天時(shí)三維航天遙感技術(shù)-介紹航天飛機(jī)雷達(dá)地形測(cè)圖計(jì)劃[J].遙感信息,2000(1):47-48.

[10] 邵裕森,戴先中.過(guò)程控制工程[M].北京:機(jī)械工業(yè)出版社,2006:139-175.

[11] 周金國(guó).InSAR數(shù)據(jù)處理方法與應(yīng)用研究[D].重慶:重慶大學(xué),2008:47-50.