劉一霖， 張勤， 黃海軍， 楊成生， 趙超英

(1.中國(guó)科學(xué)院海洋研究所海洋地質(zhì)與環(huán)境重點(diǎn)實(shí)驗(yàn)室，青島 266071； 2.長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院，西安 710054； 3.中國(guó)科學(xué)院大學(xué)，北京 100049)

礦區(qū)地表大量級(jí)沉陷形變短基線集InSAR監(jiān)測(cè)分析

劉一霖1,2,3， 張勤2， 黃海軍1,3， 楊成生2， 趙超英2

(1.中國(guó)科學(xué)院海洋研究所海洋地質(zhì)與環(huán)境重點(diǎn)實(shí)驗(yàn)室，青島 266071； 2.長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院，西安 710054； 3.中國(guó)科學(xué)院大學(xué)，北京 100049)

針對(duì)礦區(qū)地表大量級(jí)形變導(dǎo)致的InSAR影像配準(zhǔn)精度低、可監(jiān)測(cè)性差、探測(cè)量級(jí)小、地表沉陷前后完整形變信息難以獲取等問(wèn)題，研究了相應(yīng)的偏移量追蹤法、FFT過(guò)采樣法、濾波技術(shù)與基線精化等數(shù)據(jù)處理方法，并利用短基線集(small baseline subset,SBAS)技術(shù)，使距離向配準(zhǔn)精度、最大累積探測(cè)量級(jí)得到明顯提高，礦區(qū)地表形變可監(jiān)測(cè)性有了很大改善。研究結(jié)果表明，該方法不僅獲得了2008—2011年間研究區(qū)開(kāi)采進(jìn)程中地表大量級(jí)沉陷的完整形變時(shí)間序列，而且其監(jiān)測(cè)結(jié)果與外業(yè)實(shí)測(cè)數(shù)據(jù)以及采礦進(jìn)程資料具有良好的一致性； 通過(guò)對(duì)礦區(qū)地表形變剖線的統(tǒng)計(jì)分析，得到了開(kāi)采工作面地表形變的時(shí)空演變規(guī)律。

SBAS InSAR； 偏移量追蹤； 大量級(jí)形變； 采工圖； 時(shí)空演變

0 引言

礦山被開(kāi)采后，巖體應(yīng)力平衡狀態(tài)遭到破壞引起的地表連續(xù)和非連續(xù)的移動(dòng)、形變稱為開(kāi)采沉陷[1-2]。礦山資源開(kāi)采在給社會(huì)帶來(lái)巨大經(jīng)濟(jì)效益的同時(shí)，也引起礦區(qū)地形和水文變化，特別是采空區(qū)地表沉陷，破壞耕地、建筑物及其他基礎(chǔ)設(shè)施，誘發(fā)山體開(kāi)裂、崩塌、滑坡、泥石流、甚至地(礦)震等次生地質(zhì)災(zāi)害，嚴(yán)重影響礦區(qū)生態(tài)環(huán)境，制約礦區(qū)城市的可持續(xù)發(fā)展[3-4]。及時(shí)、準(zhǔn)確地監(jiān)測(cè)地表沉陷及其發(fā)展過(guò)程，研究其與礦山開(kāi)采的關(guān)聯(lián)性，對(duì)指導(dǎo)礦產(chǎn)合理開(kāi)發(fā)利用，預(yù)防并減少礦區(qū)開(kāi)采所導(dǎo)致的災(zāi)害發(fā)生，保護(hù)礦區(qū)生態(tài)環(huán)境具有重要意義。

然而，無(wú)論是傳統(tǒng)大地測(cè)量技術(shù)還是高精度全球定位技術(shù)(global positioning system, GPS)，在礦區(qū)地表形變監(jiān)測(cè)中都存在勞動(dòng)強(qiáng)度大、監(jiān)測(cè)費(fèi)用高、離散布點(diǎn)且大量級(jí)地表沉陷易導(dǎo)致監(jiān)測(cè)點(diǎn)遭到破壞等問(wèn)題，難以甚至無(wú)法對(duì)礦區(qū)地表沉陷進(jìn)行全面監(jiān)測(cè)[3-4]。差分合成孔徑雷達(dá)干涉測(cè)量技術(shù)(differential interferometric synthetic aperture radar，DInSAR)，以其全天時(shí)、全天候、大覆蓋、高分辨率、高精度以及周期性觀測(cè)等特點(diǎn)成為空間對(duì)地觀測(cè)的有效手段。國(guó)內(nèi)外諸多學(xué)者用DInSAR技術(shù)在礦區(qū)小量級(jí)(cm～dm)緩慢地表形變監(jiān)測(cè)的研究中取得良好的效果[3-8]，證明了該技術(shù)在礦區(qū)地表形變監(jiān)測(cè)應(yīng)用中的可行性。但由于礦區(qū)地表沉陷具有空間尺度小、形變梯度大、突發(fā)性強(qiáng)等特點(diǎn)，因而會(huì)導(dǎo)致DInSAR時(shí)間與空間去相干嚴(yán)重； 同時(shí)由于山區(qū)的數(shù)字高程模型(digital elevation model，DEM)往往精度偏低，使得DInSAR技術(shù)無(wú)法獲取礦區(qū)地表沉陷發(fā)生前后完整形變信息[9]，從而難以實(shí)現(xiàn)對(duì)礦區(qū)大量級(jí)(m)突發(fā)地表沉陷形變的及時(shí)有效監(jiān)測(cè)。

針對(duì)這一問(wèn)題，本文運(yùn)用GAMMA軟件[10]開(kāi)展了一系列相應(yīng)技術(shù)研究，包括利用InSAR中的短基線干涉圖集技術(shù)(small baseline subset,SBAS)以增加時(shí)空分辨率[11-13]，采用偏移量追蹤方法迭代處理以提高配準(zhǔn)精度，研究快速傅里葉變換(fast fourier transformation，F(xiàn)FT)過(guò)采樣較長(zhǎng)波長(zhǎng)數(shù)據(jù)以提高影像分辨單元增大探測(cè)量級(jí)，綜合多窗口迭代自適應(yīng)濾波技術(shù)與多方法組合的基線精化以減少形變監(jiān)測(cè)誤差。為此獲得了礦區(qū)地表完整的大量級(jí)形變時(shí)間序列，在研究區(qū)首次從時(shí)間序列上完整跟蹤監(jiān)測(cè)到采礦引起的大量級(jí)(m)地表沉陷全過(guò)程，探討了礦區(qū)開(kāi)采工作面內(nèi)地表時(shí)空演變規(guī)律，并且通過(guò)分析形變監(jiān)測(cè)結(jié)果與實(shí)地調(diào)查結(jié)果以及實(shí)際采工圖的相關(guān)性，驗(yàn)證了形變監(jiān)測(cè)結(jié)果與實(shí)地調(diào)查的一致性，揭示了與開(kāi)采過(guò)程的對(duì)應(yīng)關(guān)系。

1 短基線干涉圖集技術(shù)(SBAS)

為解決時(shí)間與空間去相干等問(wèn)題，提高InSAR技術(shù)監(jiān)測(cè)形變的時(shí)空分辨率和精度，Berardino和Lanari等[11-12]提出了SBAS技術(shù)。該技術(shù)在一定程度上克服了DEM誤差以及大氣相位的不一致性。

在(t0,∧,tN)時(shí)間內(nèi)獲取同一區(qū)域N幅SAR影像，根據(jù)干涉條件組合，生成M幅干涉圖，N/2≤M≤[N(N-1)]/2，假設(shè)用tA,tB2個(gè)時(shí)間獲得的SAR影像生成第幅干涉圖，在去除地形效應(yīng)后，像元處的解纏相位為

(1)

式中:λ為SAR影像的中心波長(zhǎng)；d(tB,x)和d(tA,x)分別為相對(duì)于參考時(shí)間t0的LOS方向累計(jì)形變；φ(tB,x)和φ(tA,x)分別為d(tB,x)和d(tA,x)引起的形變相位。利用最小二乘(least square, LS)或奇異值分解 (singular value decomposition, SVD) 對(duì)如下線性模型求取形變相位φ。模型為

Aφ=△φ,

(2)

式中:φ為待求點(diǎn)上N個(gè)時(shí)刻的SAR影像上未知形變相位組成的矩陣； △φ為M幅干涉圖上相位值組成的矩陣； 系數(shù)矩陣A[M×N]每行對(duì)應(yīng)一幅干涉圖。用式Dv=△φ將形變相位φ轉(zhuǎn)化為平均形變速率，將SVD分解應(yīng)用于矩陣D，可以得到速度矢量v的最小范數(shù)解，并求出不同SAR影像獲取期間的形變速率。

2 InSAR高精度大量級(jí)時(shí)序監(jiān)測(cè)

本文引入偏移量追蹤法迭代處理以提高影像精確配準(zhǔn)，采用較長(zhǎng)波長(zhǎng)數(shù)據(jù)迭加FFT過(guò)采樣提高可監(jiān)測(cè)量級(jí)，以及基于多窗口的迭代自適應(yīng)濾波和多方法組合的基線精化法，提高InSAR技術(shù)對(duì)礦區(qū)地表沉陷監(jiān)測(cè)的可監(jiān)測(cè)性。

2.1 偏移量追蹤方法精配準(zhǔn)

配準(zhǔn)是整個(gè)處理過(guò)程中最為基礎(chǔ)和關(guān)鍵的一步，關(guān)系到后續(xù)結(jié)果的好壞[14]。礦區(qū)特殊的地理環(huán)境以及大量級(jí)形變的失相干導(dǎo)致進(jìn)行常規(guī)的影像配準(zhǔn)容易造成嚴(yán)重配準(zhǔn)誤差，降低可監(jiān)測(cè)性，針對(duì)此問(wèn)題本文采用偏移量追蹤方法(offset tracking)[15-20]迭代處理進(jìn)行影像精確配準(zhǔn)。該方法的最大優(yōu)點(diǎn)是當(dāng)SAR影像失相干時(shí)也可以進(jìn)行配準(zhǔn)且具有較高配準(zhǔn)精度。

偏移量場(chǎng)是對(duì)窗口中實(shí)值SAR強(qiáng)度影像進(jìn)行歸一化互相關(guān)處理得到的，用大小為M1×M2(距離向方位向)的窗口對(duì)影像進(jìn)行探測(cè)與計(jì)算，并從二維互相關(guān)函數(shù)峰值所在的位置得到所需的影像距離向與方位向偏移量，采用大小為16 像元×16像元，32像元×32像元，64像元×64像元以及128像元×128像元的窗口分別求取偏移量值，并對(duì)各結(jié)果取平均以減少誤差，提高精配準(zhǔn)精度，二維相關(guān)函數(shù)R(n1,n2)定義為

(3)

(4)

采用偏移量追蹤法迭代法進(jìn)行精確配準(zhǔn)迭代處理，本文影像精配準(zhǔn)流程如圖1所示。

圖1 影像精配準(zhǔn)流程

2.2 最大探測(cè)量級(jí)

去除平地相位后并不是所有的形變都可以用DInSAR方法來(lái)探測(cè)的，如果形變的量級(jí)超過(guò)了一定形變梯度，該形變就無(wú)法被準(zhǔn)確探測(cè)到。Chen等[21]給出了最大可探測(cè)形變量的公式，即

(5)

式中：Smax,LOS為指視線方向上最大可探測(cè)沉降量;gresolution為SAR傳感器的地面分辨單元;w為形變梯度的半徑;λ為傳感器的波長(zhǎng)。

針對(duì)礦區(qū)地表大量級(jí)沉陷形變，可以通過(guò)采用波長(zhǎng)較大的數(shù)據(jù)或者提高影像地面分辨率來(lái)提高最大形變探測(cè)量級(jí)。同一區(qū)域的地表沉陷形變梯度半徑w相同。本文選用L波段波長(zhǎng)為23.5 cm的ALOS/PALSAR影像數(shù)據(jù)，從數(shù)據(jù)源上采用較長(zhǎng)波長(zhǎng)的影像數(shù)據(jù)提高探測(cè)量級(jí)； 同時(shí)對(duì)原始像元距離向×方位向大小為7.5 m×3.1 m的影像利用2.1中所述FFT采樣方法，將處理生成的SLC影像數(shù)據(jù)進(jìn)行采樣因子為2的過(guò)采樣，得到距離向像元大小為3.8 m，方位向像元大小為1.6 m的SLC影像，用于后續(xù)數(shù)據(jù)處理。對(duì)上述處理過(guò)的SLC影像，在形變梯度半徑為100 m的情況下，Smax,LOS約為1.5 m，在配準(zhǔn)精度較高時(shí)，該探測(cè)量級(jí)，對(duì)于礦區(qū)大部分地表沉陷區(qū)域都可以完整監(jiān)測(cè)。

2.3 迭代自適應(yīng)濾波

數(shù)據(jù)處理中的諸多因素都會(huì)引起干涉圖中各種噪音，從而導(dǎo)致形變結(jié)果產(chǎn)生誤差。尤其對(duì)于礦區(qū)特殊的地理環(huán)境條件以及不連續(xù)形變的監(jiān)測(cè)，采用濾波技術(shù)來(lái)抑制干涉圖的噪音是非常必要的。為此本文采用自適應(yīng)濾波進(jìn)行迭代濾波處理，以達(dá)到較好的濾波效果，基于待濾波的干涉圖本身統(tǒng)計(jì)計(jì)算得到偽相干圖(pseudo-coherent，pc)，確定濾波參數(shù)，并在迭代濾波中進(jìn)行更新，偽相干圖定義為

(6)

式中：φ為干涉圖的復(fù)相位值；N為用于計(jì)算偽相干值的窗口大小。

濾波參數(shù)確定為

(7)

2.4 基線精化

文中除運(yùn)用傳統(tǒng)的干涉條紋頻率法與基于地面控制點(diǎn)(ground control point, GCP)方法進(jìn)行基線精化處理外，對(duì)于ALOS/PALSAR衛(wèi)星影像數(shù)據(jù)，由于其軌道精度偏低，臨界基線非常長(zhǎng)，同時(shí)礦區(qū)地表沉陷量級(jí)大，解纏精度低，因此精化后的基線仍存在殘差，為此引入二次多項(xiàng)式進(jìn)行基線殘差相位的擬合和剔除，模型為

φunw=a0+a1x+a2y+a3xy+a4x2+a5y2,

(8)

式中，a0,…,a5為多項(xiàng)式系數(shù)；x,y為SAR影像的像元坐標(biāo)。而后在非礦區(qū)沒(méi)有形變的區(qū)域內(nèi)GCP提取改正后的解纏相位以及已知的高程，根據(jù)式(8)迭代計(jì)算出不同位置的基線參數(shù)B和α，達(dá)到基線精化目的，即

(9)

式中，φunw,i為提取的解纏相位；H為衛(wèi)星高度；hi為GCP高程；Ri為衛(wèi)星到各GCP距離。

綜合應(yīng)用以上處理技術(shù)，同時(shí)基于短基線干涉圖集技術(shù)(SBAS)的礦區(qū)地表大量級(jí)沉陷形變處理流程如圖2所示。

圖2 數(shù)據(jù)處理流程圖

3 實(shí)例分析

3.1 研究區(qū)與數(shù)據(jù)源

研究區(qū)位于我國(guó)陜北黃土高原北側(cè)與毛烏素沙漠東北邊緣接合的區(qū)域，海拔1 060～1 332 m之間，位于N 38.9°～39.6°，E 109.9°～110.9°之間(見(jiàn)圖3)。相關(guān)資料顯示，區(qū)域內(nèi)礦區(qū)開(kāi)采面積達(dá)到數(shù)百km2，導(dǎo)致地表發(fā)生嚴(yán)重形變，引起諸多地質(zhì)災(zāi)害，嚴(yán)重破壞了當(dāng)?shù)厣鷳B(tài)環(huán)境，威脅到該區(qū)域人民生命財(cái)產(chǎn)和工農(nóng)業(yè)生產(chǎn)的安全。

(圖中紅色實(shí)線內(nèi)為影像覆蓋范圍，紅色虛線內(nèi)為本文研究區(qū)，左下角小圖中的黑方塊為研究區(qū)位置示意)

結(jié)合研究區(qū)的地理位置與環(huán)境條件，考慮到區(qū)內(nèi)地表形變的特殊性，運(yùn)用SBAS技術(shù)對(duì)該研究區(qū)進(jìn)行地表形變監(jiān)測(cè)，選擇L波段HH極化方式，距離向與方位向像素大小分別為4.7 m與3.1 m的ALOS/PALSAR影像作為數(shù)據(jù)源。該數(shù)據(jù)具有較好的空間分辨率和較長(zhǎng)的波長(zhǎng)，受植被影響小，可以用于監(jiān)測(cè)較大量級(jí)的形變，數(shù)據(jù)詳細(xì)信息如表1所示。

表1 ALOS/PALSAR影像參數(shù)

3.2 數(shù)據(jù)處理

文中運(yùn)用偏移量追蹤法按圖1所示流程，將表1中2008年9月26日至2011年1月2日之間的11景SLC影像分別與2009年12與30日的影像精確配準(zhǔn)。為分析對(duì)比效果，將本文方法得到的距離向與方位向偏移量配準(zhǔn)精度與常規(guī)的最大干涉頻譜法處理得到的結(jié)果進(jìn)行比較，并對(duì)各影像對(duì)配準(zhǔn)精度進(jìn)行統(tǒng)計(jì)，對(duì)比結(jié)果如圖4所示。

(a) 距離向精度 (b) 方位向精度

從圖4可以看出，最大干涉頻譜法(藍(lán)色柱狀圖)配準(zhǔn)精度相對(duì)較低，特別在方位向上。偏移量追蹤方法的配準(zhǔn)精度高(綠色柱狀圖)且穩(wěn)定性好，并無(wú)論是在方位向還是在距離向都要高出最大干涉頻譜法1/10個(gè)像元，總體配準(zhǔn)精度可以提高1/10～1/5像元。對(duì)比分析可知，在由地表沉陷造成的小范圍大量級(jí)形變引起失相干的區(qū)域，采用傳統(tǒng)最大干涉頻譜法無(wú)法實(shí)現(xiàn)精確配準(zhǔn)，而本文的偏移量追蹤方法在此區(qū)域中不受相干性限制可進(jìn)行相對(duì)精確配準(zhǔn)。結(jié)合礦山開(kāi)采沉陷理論還可知，地表沉陷主要表現(xiàn)為垂向形變與位移。由于ALOS影像數(shù)據(jù)對(duì)垂向形變相對(duì)較敏感[22]，因而偏移量追蹤法特別在距離向有較高的配準(zhǔn)精度。

考慮到礦區(qū)特殊地理環(huán)境以及大量級(jí)地表形變，在對(duì)數(shù)據(jù)集內(nèi)精配準(zhǔn)后的影像進(jìn)行干涉組合時(shí)，選取時(shí)間基線小于138 d并且空間基線小于1.2 km的影像對(duì)進(jìn)行組合，最終形成了17個(gè)干涉影像對(duì)。然后利用美國(guó)SRTM-3 DEM去除與地形起伏相關(guān)的地形相位。

3.3 結(jié)果與分析

對(duì)礦區(qū)大量級(jí)沉陷形變影像數(shù)據(jù)按照本文流程進(jìn)行數(shù)據(jù)處理，得到地表大量級(jí)形變的完整監(jiān)測(cè)結(jié)果，形成的地表累積形變圖如圖5(a)所示，圖5(b)為僅采用SBAS技術(shù)常規(guī)處理方法得到的結(jié)果。

(a) 本文方法 (b) 常規(guī)處理方法

從圖5可以看出，2種處理方法監(jiān)測(cè)到的礦區(qū)地理位置與分布范圍大致相同。在獲取地表監(jiān)測(cè)量級(jí)方面，以圖中紅色橢圓所示的補(bǔ)連塔礦區(qū)一開(kāi)采工作面為例分析。僅采用SBAS方法處理，由于該開(kāi)采工作面地表存在失相干區(qū)域，部分地表沉陷形變監(jiān)測(cè)值缺失，不能得到完整的地表沉陷形變信息，監(jiān)測(cè)到的地表累積最大沉陷形變量?jī)H約70 cm； 而運(yùn)用本文處理方法不但可以監(jiān)測(cè)到更多地表沉陷形變區(qū)域，而且還監(jiān)測(cè)到補(bǔ)連塔礦區(qū)開(kāi)采工作面內(nèi)完整的地表沉陷形變信息，監(jiān)測(cè)到其最大累積沉陷形變量達(dá)到約1.2 m，與圖6所示實(shí)地考察得到的地表實(shí)際沉陷量級(jí)相一致。另外，本文還將補(bǔ)連塔礦區(qū)地表累積形變過(guò)程與礦區(qū)實(shí)際開(kāi)采進(jìn)程進(jìn)行疊加對(duì)比，從對(duì)比圖(圖7)中可以看出，監(jiān)測(cè)到礦區(qū)的位置、范圍與采礦進(jìn)程資料相一致，進(jìn)一步證實(shí)了監(jiān)測(cè)結(jié)果的可靠性。

(a) 開(kāi)采沉陷區(qū)警示牌 (b) 礦區(qū)m級(jí)沉陷 (c) 礦區(qū)階梯沉陷

圖7 監(jiān)測(cè)結(jié)果與實(shí)際采工圖對(duì)比與剖線示意

為了對(duì)礦區(qū)地表大量級(jí)沉陷形變時(shí)空演變規(guī)律進(jìn)行研究分析，特別提取出研究區(qū)中補(bǔ)連塔礦區(qū)與大柳塔礦區(qū)完整地表形變時(shí)間序列圖(圖8)。各形變圖均以2008年9月26日的影像為參考進(jìn)行地表形變累積。從圖8中可以精確得到開(kāi)采工作面內(nèi)不同時(shí)間煤層開(kāi)采進(jìn)程范圍與地表形變量級(jí)等信息。隨著時(shí)間的增加，礦區(qū)各開(kāi)采工作面內(nèi)地表沉陷形變的范圍與量級(jí)都在逐漸增加，統(tǒng)計(jì)可得補(bǔ)連塔礦區(qū)最終的地表沉陷面積約10 km2。為了更好地分析礦區(qū)地表形變的時(shí)間序列過(guò)程，將圖5中橢圓所示補(bǔ)連塔礦區(qū)開(kāi)采工作面地表形變進(jìn)行剖線分析，剖線1與剖線2位置見(jiàn)圖7，2條剖線形變時(shí)間序列曲線如圖9所示。

(a) 大柳塔礦區(qū) (b) 補(bǔ)連塔礦區(qū)

(a) 剖線1 (b) 剖線2

從圖9中可以看出，地表伴隨著緩慢形變范圍增大的同時(shí)大量級(jí)沉陷形變也逐漸出現(xiàn)； 相同時(shí)間間隔內(nèi)，地表不同區(qū)域，煤層開(kāi)采位置不同，地表形變與位移不同； 相同區(qū)域，不同時(shí)間間隔，煤層開(kāi)采時(shí)間不同，前后地表形變也存在較大差別，有緩慢形變也有突然劇烈沉陷形變。

圖9中剖線1在經(jīng)度為E110.1°附近(由圖7中可知此位置煤層開(kāi)采時(shí)間為2010年2月)的紅線位置的形變時(shí)間統(tǒng)計(jì)見(jiàn)表2。

表2 沉陷量時(shí)間序列

從表2中可以看出，在20090629—20100214時(shí)間段內(nèi)，巷道內(nèi)開(kāi)采煤層逐漸接近該區(qū)域，地表形變量相對(duì)較小，基本表現(xiàn)為緩慢形變； 在2010年2月巷道內(nèi)煤層開(kāi)采完畢后，地表發(fā)生大量級(jí)沉陷形變， 2010年7月—2010年11月間，地表又轉(zhuǎn)化為緩慢形變，2010年11月后巷道內(nèi)無(wú)煤層開(kāi)采，地表基本趨于穩(wěn)定狀態(tài)。因此，基于剖線位置地表形變規(guī)律，按時(shí)間序列可以將煤礦開(kāi)采巷道上方的地表形變大致分為4個(gè)階段： 緩慢形變、劇烈形變、緩慢形變和穩(wěn)定狀態(tài)，與礦山開(kāi)采理論中傳統(tǒng)大地測(cè)量技術(shù)監(jiān)測(cè)分析結(jié)果[1-2]相一致，由此跟蹤監(jiān)測(cè)到開(kāi)采進(jìn)程中地表沉陷過(guò)程的完整信息。

4 結(jié)論

常規(guī)InSAR技術(shù)難以實(shí)現(xiàn)對(duì)礦區(qū)地表大量級(jí)突發(fā)沉陷形變的全面監(jiān)測(cè)，針對(duì)這一問(wèn)題，本文研究了相應(yīng)的數(shù)據(jù)處理技術(shù)，可有效提高礦區(qū)地表形變的可監(jiān)測(cè)性，取得如下結(jié)論：

1)偏移量追蹤方法迭代配準(zhǔn)可以很好地提高礦區(qū)地表小范圍大量級(jí)沉陷形變?cè)斐傻氖喔蓞^(qū)域配準(zhǔn)精度； FFT過(guò)采樣方法能夠較好地提高影像分辨單元，增加礦區(qū)地表沉陷探測(cè)量級(jí)； 多窗口迭代自適應(yīng)濾波技術(shù)在保存大量級(jí)地表形變的同時(shí)較好地去除了噪聲信息，改善了礦區(qū)地表形變結(jié)果的監(jiān)測(cè)質(zhì)量； 基于基線殘差擬合與剔除的基線精化方法對(duì)地理環(huán)境特殊的礦區(qū)能夠有效改善ALOS影像基線精度，進(jìn)而提高監(jiān)測(cè)結(jié)果質(zhì)量及效果。

2)礦區(qū)開(kāi)采沉陷的實(shí)際監(jiān)測(cè)結(jié)果表明，融合以上技術(shù)的短基線集(SBAS)技術(shù)，可以監(jiān)測(cè)完整礦區(qū)地表沉陷形變，其結(jié)果與外業(yè)調(diào)查以及實(shí)際采工圖等資料具有良好一致性。

3)結(jié)合文中InSAR監(jiān)測(cè)結(jié)果分析可得到礦區(qū)開(kāi)采工作面內(nèi)地表位移與形變的時(shí)空演變規(guī)律為： 緩慢形變、劇烈形變、緩慢形變和穩(wěn)定狀態(tài)，與礦山開(kāi)采理論相吻合。

本文研究雖可以提高礦區(qū)地表沉陷的可監(jiān)測(cè)性，監(jiān)測(cè)得到與開(kāi)采進(jìn)程相對(duì)應(yīng)的完整形變信息，但對(duì)于開(kāi)采工作面內(nèi)中心位置地表沉陷的最大量級(jí)值的精確監(jiān)測(cè)還存在一定誤差甚至粗差，這將是今后礦區(qū)地表大量級(jí)沉陷InSAR監(jiān)測(cè)的重點(diǎn)與難點(diǎn)。

[1] 何國(guó)清,楊倫,凌賡娣,等.礦山開(kāi)采沉陷學(xué)[M].徐州:中國(guó)礦業(yè)大學(xué)出版社,1991. He G Q,Yang L,Ling G D,et al.Mining Subsidence[M].Xuzhou:China University of Mining and Technology Press,1991.

[2] 劉寶琛,廖國(guó)華.煤礦地表移動(dòng)的基本規(guī)律[M].北京:中國(guó)工業(yè)出版社,1965. Liu B C,Liao G H.Basic Rule of Coal Surface Movement[M].Beijing:China Industry Press,1965.

[3] 吳立新,高均海,葛大慶,等.基于D-InSAR的煤礦區(qū)開(kāi)采沉陷遙感監(jiān)測(cè)技術(shù)分析[J].地理與地理信息科學(xué),2004,20(2):22-25,37. Wu L X,Gao J H,Ge D Q,et al.Technical analysis of the remote sensing monitoring for coal-mining subsidence based on D-InSAR[J].Geography and Geo-Information Science,2004,20(2):22-25,37.

[4] 吳立新,高均海,葛大慶,等.工礦區(qū)地表沉陷D-InSAR監(jiān)測(cè)試驗(yàn)研究[J].東北大學(xué)學(xué)報(bào):自然科學(xué)版,2005,26(8):778-782. Wu L X,Gao J H,Ge D Q,et al.Experimental study on surface subsidence monitoring with D-InSAR in mining area[J].Journal of Northeastern University:Natural Science,2005,26(8):778-782.

[5] 楊成生.基于D-InSAR技術(shù)的煤礦沉陷監(jiān)測(cè)[D].西安:長(zhǎng)安大學(xué),2008. Yang C S.Based on D-InSAR Technology in Coal Mine Subsidence Monitoring[D].Xi’an:Chang’an University,2008.

[6] 朱建軍,邢學(xué)敏,胡俊,等.利用InSAR技術(shù)監(jiān)測(cè)礦區(qū)地表形變[J].中國(guó)有色金屬學(xué)報(bào),2011,21(10):2564-2576. Zhu J J,Xing X M,Hu J,et al.Monitoring of ground surface deformation in mining area with InSAR technique[J].The Chinese Journal of Nonferrous Metals,2011,21(10):2564-2576.

[7] Ge L L,Chang H C,Rizos C.Mine subsidence monitoring using multi-source satellite SAR images[J].Photogrammetric Engineering & Remote Sensing,2007,73(3):259-266.

[8] Ng A H M,Ge L L,Yan Y G,et al.Mapping accumulated mine subsidence using small stack of SAR differential interferograms in the Southern coalfield of New South Wales, Australia[J].Engineering Geology,2010,115(1/2):1-15.

[9] 劉一霖.礦區(qū)開(kāi)采沉陷大量級(jí)形變監(jiān)測(cè)與反演分析[D].西安:長(zhǎng)安大學(xué),2013. Liu Y L.The Monitoring and Inverse of Large Scale Deformation of Minging Subsidence[D].Xi’an:Chang’an University,2013.

[10]Werner C,Wegmüller U,Strozzi T,et al.Gamma SAR and interferometric processing software[C]//Proceedings of the ERS-ENVISAT Symposium.Gothenburg,Sweden:[s.n.],2000.

[11]Berardino P,Fornaro G,Lanari R,et al.A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms[J].IEEE Transactions on Geoscience and Remote Sensing,2002,40(11):2375-2383.

[12]Lanari R,Mora O,Manunta M,et al.A small-baseline approach for investigating deformations on full-resolution differential SAR interferograms[J].IEEE Transactions on Geoscience and Remote Sensing,2004,42(7):1377-1386.

[13]Lanari R,Casu F,Manzo M,et al.An overview of the small baseline subset algorithm:A DInSAR technique for surface deformation analysis[J].Pure and Applied Geophysics,2007,164(4):637-661.

[14]廖明生,林琿.雷達(dá)干涉測(cè)量:原理與信號(hào)處理基礎(chǔ)[M].北京:測(cè)繪出版社,2003:45-52. Liao M S,Lin H.Synthetic Aperture Radar Interferometry-Principle and Signal Processing[M].Beijing:Surveying and Mapping Press,2003:45-52.

[15]Strozzi T,Luckman A,Murray T,et al.Glacier motion estimation using SAR offset-tracking procedures[J].IEEE Transactions on Geoscience and Remote Sensing,2002,40(11):2384-2391.

[16]Strozzi T,Kouraev A,Wiesmann A,et al.Estimation of Arctic glacier motion with satellite L-band SAR data[J].Remote Sensing of Environment,2008,112(3):636-645.

[17]Werner C,Strozzi T,Wiesmann A,et al.Complimentary measurement of geophysical deformation using repeat-pass SAR[C]//Proceedings of 2001 IEEE International Geoscience and Remote Sensing Symposium.Sydney,NSW:IEEE,2001,7:3255-3258.

[18]Werner C,Wegmüller U,Strozzi T,et al.Precision estimation of local offsets between pairs of SAR SLCs and detected SAR images[C]//Proceedings of 2005 IEEE International Geoscience and Remote Sensing Symposium.Seoul:IEEE,2005:4803-4805.

[19]Wegmüller U,Werner C,Strozzi T,et al.Automated and precise image registration procedures[J].Analysis of Multi-temporal Remote Sensing Images,2002,2:37-49.

[20]Wegmüller U,Werner C,Strozzi T,et al.Ionospheric electron concentration effects on SAR and INSAR[C]//Proceedings of 2006 IEEE International Conference on Geoscience and Remote Sensing Symposium.Denver,CO,USA:IEEE,2006:3731-3734.

[21]Chen C W,Zebker H A.Phase unwrapping for large SAR interferograms:Statistical segmentation and generalized network models[J].IEEE Transactions on Geoscience and Remote Sensing,2002,40(8):1709-1719.

[22]Lu Z,Wicks Jr C.Characterizing 6 August 2007 crandall canyon mine collapse from ALOS PALSAR InSAR[J].Geomatics,Natural Hazards and Risk,2010,1(1):85-93.

(責(zé)任編輯： 李瑜)

Monitoring and analyzing large scale land subsidence over the mining area using small baseline subset InSAR

LIU Yilin1,2,3, ZHANG Qin2, HUANG Haijun1,3, YANG Chengsheng2, ZHAO Chaoying2

(1.KeyLaboratoryofMarineGeologyandEnvironment,InstituteofOceanology,ChineseAcademyofSciences,Qingdao266071,China; 2.CollegeofGeologyEngineeringandGeomatics,Chang’anUniversity,Xi’an710054,China; 3.UniversityofChineseAcademyofSciences,Beijing100049,China)

Due to large scale earth surface deformation, the application of conventional InSAR technique to monitor land subsidence over the mining area has many limitations, such as low image co-registration accuracy and monitoring capability, small detection scale and unavailable complete mining subsidence information. In view of such a situation, the small baseline subset (SBAS) InSAR technique combined with offset tracking method, fast fourier transformation oversampling technique, filter technique and baseline refine method was studied in this paper to overcome the limitations. On such a basis, the co-registration accuracy, monitoring capability and the accumulative detection scale could be improved considerably. Meanwhile, the complete large scale time series deformation over the mining area from 2008 to 2011 was generated, which is well consistent with field and mining processing data. Furthermore, spatial and temporal evolution law of earth surface over the mining area was obtained by analyzing the cross-section time series deformation.

SBAS InSAR; offset tracking; large scale deformation; mining working diagram; spatial and temporal evolution law

10.6046/gtzyyg.2017.02.21

劉一霖,張勤,黃海軍,等.礦區(qū)地表大量級(jí)沉陷形變短基線集InSAR監(jiān)測(cè)分析[J].國(guó)土資源遙感,2017,29(2)：144-151.(Liu Y L,Zhang Q,Huang H J，et al.Monitoring and analyzing large scale land subsidence over the mining area using small baseline subset InSAR[J].Remote Sensing for Land and Resources,2017,29(2)：144-151.)

2015-11-20；

2015-12-30

國(guó)家自然科學(xué)基金項(xiàng)目“多分辨率雷達(dá)干涉融合技術(shù)用于礦區(qū)塌陷災(zāi)害研究”(編號(hào)： 41072266)、國(guó)家自然科學(xué)基金項(xiàng)目“黃河三角洲地面沉降監(jiān)測(cè)與形成機(jī)理研究”(編號(hào)： 41276082)、國(guó)家自然科學(xué)基金項(xiàng)目“渤海西南岸末次冰期以來(lái)古海岸線重建及影響因素分異研究”(編號(hào)： 41106041)、中國(guó)科學(xué)院知識(shí)創(chuàng)新工程重要方向性項(xiàng)目“現(xiàn)代黃河口三角洲垂向形變驅(qū)動(dòng)機(jī)制與環(huán)境效應(yīng)研究”(編號(hào)： KZCX2-EW-207)及國(guó)土資源大調(diào)查項(xiàng)目“汾渭盆地重點(diǎn)地區(qū)地面沉降地裂縫InSAR與GPS監(jiān)測(cè)”(編號(hào)： 1212011120069)共同資助。

劉一霖(1986-)，男，博士生，主要從事地質(zhì)遙感InSAR研究與應(yīng)用，Email： lyilin@msn.com。

張勤(1958-)，女，教授，博士生導(dǎo)師，主要從事GPS、InSAR及地質(zhì)災(zāi)害監(jiān)測(cè)方面的研究。 Email: zhangqinle@263.net.cn。

TP 79

A

1001-070X(2017)02-0144-08