利用過采樣方法提高礦區(qū)沉降干涉圖的相干性

吳海翔,王 華,黃秀汝,Carolina Pagli,楊澤鑫

(1.廣東工業(yè)大學(xué) 土木與交通工程學(xué)院， 廣東 廣州 510006；2.比薩大學(xué) 地球科學(xué)系, 意大利 比薩 56126)

利用過采樣方法提高礦區(qū)沉降干涉圖的相干性

吳海翔1,王 華1,黃秀汝1,Carolina Pagli2,楊澤鑫1

(1.廣東工業(yè)大學(xué) 土木與交通工程學(xué)院， 廣東 廣州 510006；2.比薩大學(xué) 地球科學(xué)系, 意大利 比薩 56126)

采用過采樣方法提高SAR圖像相干性，使研究區(qū)域的干涉條紋更加清晰.以唐山礦區(qū)為例，采用過采樣方法獲得了6對礦區(qū)沉降干涉圖，并計算了該區(qū)域的地面沉降時間序列.結(jié)果表明，唐山部分礦區(qū)在2004年11月到2005年4月間存在快速的沉降，最大沉降速率達到2 398 mm/a.

雷達干涉測量； 礦區(qū)沉降； 過采樣； 相干性； 唐山

目前，利用雷達干涉測量(Interferometric Synthetic Aperture Radar, InSAR)進行地面沉降監(jiān)測多集中在自然地層變化或地下水開采等引起的地面沉降方面，這種地面沉降與礦產(chǎn)資源開采所引起的地面沉降的機理有所不同，前者一般發(fā)生在較大區(qū)域范圍內(nèi)且沉降速率比較緩慢平穩(wěn)，而后者一般發(fā)生在范圍相對較小的礦區(qū)工作面上，形成塌陷盆地.因此，監(jiān)測礦區(qū)地面沉降的 SAR 圖像要求具有比較高的質(zhì)量(即數(shù)據(jù)的充足性和連續(xù)性).由于礦區(qū)沉降中心的變形梯度經(jīng)常超過InSAR所能探測的最大值[1-2]，因此，經(jīng)常發(fā)生去相干現(xiàn)象[3]，致使干涉圖的條紋混亂，無法進行相位解纏.本文基于升余弦(raised cosine，RC)插值算法[4-5]，在形成干涉圖前對SAR圖像進行插值重采樣處理，提高了圖像的空間分辨率并顯著改善了相干性，從而最小化解纏錯誤，得到較可靠的相位估計.

1 處理策略

1.1 數(shù)據(jù)選擇與SLC圖像的生成

對InSAR數(shù)據(jù)進行處理時，根據(jù)垂直基線和時間基線來選取SAR圖像組成干涉像對，從而形成干涉圖網(wǎng)絡(luò)[4,6-7]，利用ROI_PAC軟件[8-9]將干涉像對源數(shù)據(jù)分別處理成SLC(Single Look Complex,單視復(fù)數(shù))圖像.本文所使用的數(shù)據(jù)是ENVISAT衛(wèi)星數(shù)據(jù)，波長5.6 cm，距離向和方位向分辨率分別為20 m和4 m.

1.2 過采樣

在形成干涉圖之前，對于SLC圖像的重采樣可應(yīng)用sinc型插值函數(shù)[10].Cho等人提出使用RC插值內(nèi)核.結(jié)合sinc型插值函數(shù)和升余弦函數(shù)，得到2維的脈沖響應(yīng)[11]

i(x,y)=

(1)

基于奈奎斯特準(zhǔn)則，Cho等人給出了α和β的最佳估值[12].

(2)

(3)

式(2)、(3)中，Xr和Xa分別是距離向和方位向的過采樣因子，即采樣頻率除以帶寬.這里取Xr= 1.200 5和Xa=1.158 8.

2 InSAR 數(shù)據(jù)處理

本實驗基于ENVISAT SAR數(shù)據(jù),利用ROI_PAC軟件處理ENVISAT衛(wèi)星175軌道2004年11月到2005年4月間獲取的唐山礦區(qū)5幅衛(wèi)星圖像，生成了6幅干涉圖，獲取時間及有關(guān)參數(shù)如表1所示.數(shù)據(jù)處理時采用二軌差分，利用美國的SRTM DEM去除地形相位[13].為了編碼更高分辨率的雷達幅度影像，這里采用的DEM是由原始分辨率為3弧秒的SRTM DEM(分辨率大約為90 m)插值成1弧秒所生成的.

表1 干涉圖像對

在InSAR 數(shù)據(jù)處理過程中，對于去相關(guān)現(xiàn)象，一般的處理方法是對圖像進行多視(multi-look)操作[14]，即對復(fù)干涉圖做空域平均，此操作在提高圖像信噪比的同時犧牲了空間分辨率[15].實驗表明，對圖像做N(如N=4)視處理對于城市的整體沉降監(jiān)測效果較好，但對于礦區(qū)的沉降監(jiān)測就顯得略為不足.如圖1(a)所示，對實驗圖像做4視處理，所生成的干涉圖距離向和方位向分辨率均為80 m.圖中可見礦區(qū)沉降中心外沿條紋較為清晰，但沉降中心由于沉降梯度大，去相干現(xiàn)象嚴(yán)重，造成相位信息缺失，直接影響后續(xù)相位解纏精度[16]，不利于礦區(qū)沉降中心沉降量的提取.

針對干涉圖中沉降中心條紋混疊的現(xiàn)象，本實驗對InSAR 數(shù)據(jù)進行1 視處理，生成的干涉圖如圖1(b)所示，礦區(qū)沉降中心的干涉條紋相對清晰，沉降區(qū)域外沿條紋及條紋分界更加明顯.比較圖1(a)和圖1(b)，可以發(fā)現(xiàn)，對于沉降中心，經(jīng)過 4 視處理的干涉圖中的條紋混疊，但在1視的干涉圖(不經(jīng)過多視操作)中，原來混疊的條紋被分開，即使分界仍不是很明顯.

圖1 纏繞干涉圖

對InSAR數(shù)據(jù)做1視處理雖然可以使條紋更加清晰，但如上所述，沉降中心去相干現(xiàn)象嚴(yán)重，相位信息缺失.因此，在對InSAR數(shù)據(jù)做 1 視處理的基礎(chǔ)上，利用RC插值算法對SLC像對進行過采樣處理.本實驗中，方位向過采樣因子為1，距離向過采樣因子為5.經(jīng)過過采樣處理后的圖像在方位向的分辨率不變，仍大約為4 m，距離向的分辨率由原來的大約20 m提高為大約4 m.將經(jīng)過過采樣處理后的SLC 像對進行干涉處理，可直觀發(fā)現(xiàn)，干涉圖中礦區(qū)沉降中心條紋及分界較于上述兩種方法清晰、明顯(如圖1(c))，從而實現(xiàn)沉降條紋(相位)的恢復(fù).

通過對圖1礦區(qū)沉降中心的干涉條紋進行對比，可得出結(jié)論：對SLC像對進行過采樣處理后形成的干涉圖具有更高的干涉條紋清晰度和更高的質(zhì)量.

3 礦區(qū)沉降時空演變分析

3.1 礦區(qū)沉降空間分布特征分析

形變干涉相位變化直接反映了地表的變化[17]，以像對20041121-20041226為例，如圖2(a)所示，圖中黑色線(即剖面線)穿過礦區(qū)的主要沉降中心，沿線提取其沉降量，見圖2(a).

圖2為干涉像對的二維變形圖，可大致確定變形區(qū)域.圖3中，提取的各剖面所在礦區(qū)在35 d(2004 年 11 月 21 日到 2004 年 12 月 26 日)中均發(fā)生了不同程度的沉降.其中，剖面 1-1′沉降變化較為明顯，其值基本分布在50～100 mm之間，沉降最大值約為120 mm；剖面2-2′所在礦區(qū)在剖面開始1 km范圍內(nèi)出現(xiàn)較大沉降，最大值達110 mm，后續(xù)沉降梯度較為平緩；剖面3-3′所在礦區(qū)是發(fā)生地表沉降較為顯著的地段，在剖面開始1.4 km處有一個最大沉降值為230 mm的沉降漏斗，沉降率高達2 398 mm/a，此后，沉降值集中在80～120 mm之間.

圖2 二維變形圖

圖3 各剖面變形圖(20041121-20041226)

3.2 礦區(qū)沉降時域演化分析

本實驗利用研究區(qū)的6幅干涉圖進行時域分析，提取圖中相對穩(wěn)定區(qū)的A點及礦區(qū)沉降中心的B，C點進行時域演化分析，如圖2(a)所示.從圖4中可發(fā)現(xiàn)，A點在2004年11月到2005年4月的整個研究期間，沉降增量接近于0，說明此點受礦區(qū)活動影響較小，為穩(wěn)定點；B點上，最大沉降增量在2004年11月到同年12月期間，達到約75 mm，此時期后到各個時域分析階段沉降增量逐漸變??；對于C點，各個時域分析階段沉降量相對大，最大沉降發(fā)生在2004年11月到同年12月期間，沉降量為92 mm，另外3個時刻的沉降增量在50 mm左右.

圖4 干涉圖時域分析

4 結(jié)論與展望

本文通過實驗表明，在形成干涉圖前對SLC圖像進行過采樣處理，可以使礦區(qū)沉降中心干涉條紋更加清晰、分界更加明顯，對于提高干涉圖的質(zhì)量有顯著的作用.將此方法應(yīng)用于唐山礦區(qū)，獲取了唐山礦區(qū)的沉降干涉圖，從圖中可看出唐山局部礦區(qū)在2004 年11月到 2005年4月間有不同程度的沉降.

[1] Baran I,Stewart M,Claessens S.A new func-tional model for determining minimum and maximum detectable deformation radient re-solved by satellite radar interferometry [J].IEEE Trans Geosci Remote Sens,2005,43:675-682.

[2] Ge L, Wang H, Chang H C et al. Linear conbination for differential radar interferometry[J].Observing our Changing Earth,2009,133: 825-829.

[3] 王華, 羅麗芳. 利用InSAR相干性提取青藏高原湖泊邊界[J].廣東工業(yè)大學(xué)學(xué)報,2014 ，31 (1): 118-120.

Wang H, Luo L F. Identifying lake boundaries in Tibet using InSAR coherence[J]. Journal of Guangdong University of Technology，2014，31(1): 118-120.

[4] Berardino P G, Fornaro, Lanari R, et al. A new algorithm for surface deformation monitoring based on small baseline differr-ential SAR interferograms[J]. IEEE Trans Geosci Remote Sens,2002,40(11):2375-2383.

[5] Lanari R, Lundgren P, Manzo M, et al.Sa-tellite radar interferometry time series analysis of surface deformation for Los Angeles ,California[J]. Geophys Res Lett, 2004,31 (23):L23613.

[6] Hooper A,Zebker H, Segall P, et al. A new method for measuring deformation on volca-noes and other natural terrains using InSAR persistent scatterers[J].Geophys Res Lett,2004,31: L23611.

[7] Wang H, Wright T J,Yu Y P. InSAR reveals coastal subsidence in the Pearl River De-lta,China[J].Geophysical Journal International,2012(3):1119-1128.

[8]王華, 彭佳卉. 基于ECMWF模型改正InSAR大 氣延遲誤差的可靠性實驗[J]. 廣東工業(yè)大學(xué)學(xué)報，2014，31(2):74-77.

Wang H, Peng J H. Test on InSAR atmospheric delay correction using ECMWF model[J]. Journal of Guangdong University of Technology，2014，31(2):74-77.

[9] Rosen P A,Hensley S,Peltzer G, et al. Updated repeat orbit interferometry packa-ge released[J]. Eos Trans,2004,85(5):47.

[10] Desai M D,Jenkins W K. Convolution back-projection image reconstruction for spot-light mode synthetic aperture radar [J]. IEEE Trans Geosci Remote Sens,1992(30): 505-515.

[11] Yun S H,Zebker H,Segall P,et al.Interf-erogram formation in the presence of comp-lex and large deformation [J].Geophys Res Lett,2007(34):L12305.

[12] Cho B L,KongY K, Kim Y S.Interpolation using optimum Nyquist filter for SAR interferometry [J].Electromagn Waves,2005, 19(1),129-135.

[13] Farr T G, Rosen P A, Caro E. The shuttle radar topography mission [J]. Rev Geophys, 2007,45(2):doi:10.1029/2005RG000183.

[14] Casu F, Buckley S M,Manzo M. Large scale InSAR deformation time series: Phoenix and Houston case studies [J]. IEEE Int Geosci Remote Sens,2005(7):5240-5243.

[15] Hanssen R F.Radar Interferometry: Data Interpretation and Error Analysis[M]. Dordrecht: Kluwer Academic Publishers, 2001.

[16] Wegmüller U, Walter D, Spreckels V, et al.Nonuniform ground motion monitoring with TerraSAR-X persistent scatterer interferometry [J]. IEEE Trans Geosci Remote Sens,2009,48(2):895-904.

[17] Wang H, Wright T J, Biggs J. Interseismic slip rate of the northwestern Xianshuihe fault from InSAR data[J]. Geophys Res Lett,2009(36):L03302.

Improving the Coherence of Interferograms of Mining Subsidence by Oversampling

Wu Hai-xiang1, Wang Hua1, Huang Xiu-ru1, Carolina Pagli2, Yang Ze-xin1

(1.School of Civil and Transportation Engineering, Guangdong University of Technology, Guangzhou 510006, China; 2.Department of Earth Science, University of Pisa, Pisa 56126, Italy)

The authors use oversampling method to improve the coherence of SAR images, which makes the fringes of the studied area clearer. By using this method, the authors process 6 pairs of interferograms and derive the time series of mining subsidence. The results show that some mines in Tangshan have been subsiding rapidly from November 2004 to April 2005 with the maximum rate of 2 398 mm/yr.

InSAR; mining subsidence; oversampling; coherence; Tangshan

2014- 04- 11

國家自然科學(xué)基金青年基金資助項目(41104016)

吳海翔(1988-)，男，碩士研究生，主要研究方向為雷達干涉測量.

10.3969/j.issn.1007- 7162.2015.03.023

TU196

A

1007-7162(2015)03- 0123- 04