基于合成孔徑雷達差分干涉測量技術的礦區地面形變分析①

張文展

(新鄉航空工業(集團)有限公司,河南 新鄉453049)

0 引 言

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

1 D-InSAR測量基本原理

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

式中,λ為雷達信號的波長。

由圖1中的幾何關系模型和余弦定理可得:

圖1 差分干涉測量成像幾何示意圖

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

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

在地表發生形變后,進行觀測獲取的相位差不但與地形有關,而且與沿雷達視線方向的形變量ΔRd有關。此時對目標點P的測量相位差為

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

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

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

2 數據源以及數據處理的流程

2.1 數據源

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

表1 ASAR像對的基線列表

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

2.2 數據處理流程

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

圖2 雙軌法差分干涉測量處理流程

3 D-InSAR結果及其礦區地面形變分析

3.1 D-InSAR處理結果

1)差分干涉圖

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

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

2)地面沉降圖

對差分干涉圖進行相位解纏,并進行相高轉換后,即可得到沿雷達視線方向的形變圖,即為沉降圖。

在圖4的沉降圖可以看到所研究區域的某些區域形成了明顯的“漏斗”沉降;某些區域發生了一些微小的沉降可能是由于開采地下水或數據處理中不可避免的誤差造成;有些區域發生了上升現象,這可能與開采礦物有關,也可能與數據處理中有關參數設置不當有關。

3.2 礦區地面形變分析

根據礦區范圍,從整體上分析該研究區域不同時間段的沉降量、沉降分布及其沉降面積統計。

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

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

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

圖4 雙軌法D-InSAR沉降圖

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

3.3 精度分析

為了驗證D-InSAR監測礦區地面沉降的精度,收集了水準監測手段獲取2004.11.20-2005.12.31時間段內所研究礦區地面沉降監測值,與DInSAR監測結果相比較。表2和圖9給出了DInSAR監測結果與水準結果的比較情況。

圖9 水準結果與D-InSAR監測結果比較

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

4 結 論

采用雙軌法D-InSAR處理覆蓋濟北某礦區的ENVISATASAR數據,提取了增強差分干涉圖、相干圖、沉降圖,并從礦區地面沉降分析和精度分析等方面對D-InSAR監測礦區地面沉降結果進行了分析討論,可以得出:

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

2)由于受到ASAR分辨率和成像波長等因素制約,造成D-InSAR技術監測礦區地面沉降結果與水準監測沉降結果之間存在一定差異。

[1] 路 旭.用InSA R作地面沉降監測的實驗研究[J].大地測量與地球動力學,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誤差與軌道誤差對相位貢獻的靈敏度研究[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] 鄒積亭,劉遠明.D-InSAR在地面形變監測中的應用研究[J].北京測繪,2009(2):19-20.

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

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

[9] 郭華東,王長林.全天候全天時三維航天遙感技術-介紹航天飛機雷達地形測圖計劃[J].遙感信息,2000(1):47-48.

[10] 邵裕森,戴先中.過程控制工程[M].北京:機械工業出版社,2006:139-175.

[11] 周金國.InSAR數據處理方法與應用研究[D].重慶:重慶大學,2008:47-50.