基于四叉樹濾波的InSAR礦區(qū)形變監(jiān)測(cè)與分析

楊小鴿，王 華，柏 俊，吳希文，李廷俊，柏玉超，皮廷亮

（1.廣東工業(yè)大學(xué) 土木與交通工程學(xué)院，廣東 廣州 510006；2.廣州海川信息科技有限公司，廣東 廣州 510300；3.云南華聯(lián)鋅銦股份有限公司，云南 文山 663700）

礦區(qū)開采導(dǎo)致的地表形變易造成礦區(qū)塌陷，破壞礦區(qū)建筑設(shè)施，影響工作人員的生命安全，因而礦區(qū)形變監(jiān)測(cè)是礦區(qū)開采的一項(xiàng)重要工作[1-3]。相比傳統(tǒng)礦區(qū)形變監(jiān)測(cè)手段，InSAR不僅能滿足毫米級(jí)的精度需求，同時(shí)憑借其大尺度和高分辨率的特點(diǎn)，InSAR還能同時(shí)獲得多個(gè)礦區(qū)米級(jí)空間分辨率的形變圖，對(duì)于礦區(qū)形變監(jiān)測(cè)具有明顯的優(yōu)勢(shì)[4-6]。時(shí)間序列InSAR能有效地監(jiān)測(cè)礦區(qū)形變過程，但礦區(qū)周圍區(qū)域植被茂密，容易導(dǎo)致時(shí)變?nèi)ハ喔?，形變結(jié)果中會(huì)存在大量的離散噪聲點(diǎn)。對(duì)于露天開挖的礦區(qū)，地表屬性變化較快，礦區(qū)內(nèi)長期穩(wěn)定點(diǎn)很少，難以獲得長期穩(wěn)定的永久散射體和分布式散射體。為了獲得高分辨率的InSAR形變時(shí)間序列，需要對(duì)原始短時(shí)間基線干涉圖的解纏結(jié)果進(jìn)行濾波，從而提高InSAR形變時(shí)間序列結(jié)果的質(zhì)量。由于能夠在降低數(shù)據(jù)量的同時(shí)保留形變場的主要特征，四叉樹采樣方法被廣泛地用在InSAR形變建模中[7-8]。在四叉樹采樣思想的基礎(chǔ)上，Peng等[9]提出了一種四叉樹濾波方法，能夠有效消除噪聲點(diǎn)，并保留形變特征，在2016年新西蘭凱庫拉地震的偏移量分析中得到了較好的應(yīng)用[9-10]。

本文研究的云南某礦區(qū)是一個(gè)露天開采礦，屬亞熱帶季風(fēng)氣候，年降水量大，礦區(qū)周圍植被較茂密。礦區(qū)范圍內(nèi)有采礦場、選礦車間、排石場3個(gè)主要區(qū)域,目前主要采用陡幫作業(yè)采剝工藝。該礦區(qū)形變范圍小、形變梯度大，與地震形變的震中區(qū)具有相似的特性，而對(duì)分辨率的要求比地震更高，在高分辨率下進(jìn)行相位解纏往往會(huì)帶來較多的噪聲點(diǎn)。圖1為20200620-20200702干涉對(duì)的解纏相位圖，其中黑線表示3個(gè)典型剖面的位置(見圖4)，黑框表示地面觀測(cè)點(diǎn)的分布范圍(見圖5)。從圖1可以看出該礦區(qū)干涉圖受時(shí)變?nèi)ハ喔捎绊懚休^多噪聲。為了提高InSAR形變結(jié)果的質(zhì)量，本文對(duì)該礦區(qū)2019年10月至2021年3月的InSAR解纏結(jié)果進(jìn)行四叉樹濾波，并通過時(shí)間序列分析改正軌道誤差、大氣延遲誤差，最終獲得研究區(qū)域干凈的時(shí)間序列形變結(jié)果。

圖1 20200620-20200702干涉對(duì)的解纏相位圖Fig.1 Unwrapped phase of the 20200620-20200702 interferogram

圖2 濾波效果圖和剔除的像素點(diǎn)(a-b分別為粗差閾值為2倍中誤差，最小劃分窗口為全局64×64的濾波后解纏相位圖和剔除的像素點(diǎn)；c-d為粗差閾值為2倍中誤差，最小窗口為全局8×8的濾波效果圖對(duì)應(yīng)圖；e-f為粗差閾值為2倍中誤差，近場8×8和遠(yuǎn)場64×64的濾波效果圖對(duì)應(yīng)圖；g-h為最小窗口為近場8×8和遠(yuǎn)場64×64，粗差閾值為1.5倍中誤差的濾波效果對(duì)應(yīng)圖；i-j為最小窗口為近場8×8和遠(yuǎn)場64×64，粗差閾值為3倍中誤差的濾波效果對(duì)應(yīng)圖)Fig.2 Performance of quadtree filtering.(a and b are the filtered unwrapped phase and the removed pixels, respectively, with the gross error threshold of 2 times of the standard deviation and the minimum division window of 64×64; c-d are the corresponding results with the gross error threshold of 2 times of the standard deviation and the window of 8×8; e-f are the corresponding results with the gross error threshold of 2 times of the standard deviation, but with different minimum window size in the near (8×8) and the far field (64×64) ; g-h and i-j are the corresponding results with the gross error threshold of 1.5 and 3 times of the standard deviation, respectively, and the minimum window size as used in e-f.)

圖3 礦區(qū)時(shí)間序列形變結(jié)果(遠(yuǎn)離衛(wèi)星為正)Fig.3 Time series deformation results of mining area (Positive indicates away from the satellite)

圖4 礦區(qū)剖面形變時(shí)空演化圖(遠(yuǎn)離衛(wèi)星為正)Fig.4 Spatio-temporalevolution map of deformation along profile(Positive indicates away from the satellite)

圖5 礦區(qū)形變結(jié)果驗(yàn)證圖(a)測(cè)量機(jī)器人監(jiān)測(cè)點(diǎn)位圖；(b)InSAR與32個(gè)測(cè)量機(jī)器人點(diǎn)的時(shí)間序列形變結(jié)果對(duì)比分析圖，顏色表示差值的大??；(c-e)為P1、P2、P3三個(gè)點(diǎn)的InSAR與測(cè)量機(jī)器人時(shí)間序列形變結(jié)果對(duì)比圖Fig.5 Validation between InSAR and total stations (a) Location of 32 total stations in black dots; (b) Scatter plot of the LOS displacement time series from InSAR and total stations, the color-coded circles represent LOS displacement difference of these two means; (c-e) Time series of three points P1, P2 and P3 from InSAR and total stations

1 四叉樹濾波

1.1 四叉樹濾波原理

四叉樹濾波能剔除影像中離散的噪聲點(diǎn)。該濾波方法首先將影像進(jìn)行遞歸分級(jí)，每次將區(qū)域分成4個(gè)大小相等的子區(qū)域，計(jì)算子區(qū)域內(nèi)像素值的中誤差或減去多項(xiàng)式模型后的中誤差，如果標(biāo)準(zhǔn)差超過某一給定的閾值，則繼續(xù)往下將各子區(qū)域進(jìn)行分割，直到各子區(qū)域的中誤差滿足閾值要求，或者該子區(qū)域已劃分至預(yù)設(shè)的最小劃分窗口。在四叉樹采樣中，當(dāng)各子區(qū)域劃分達(dá)到上述標(biāo)準(zhǔn)后則輸出該窗口的平均值(或中值)，而四叉樹濾波則是按照一定原則剔除該窗口殘差超過一定閾值的像素，并將剩余的所有像素值輸出[10]。通常，最小劃分窗口和粗差閾值對(duì)四叉樹濾波的影響較為顯著。下文將比較不同的最小劃分窗口和粗差閾值對(duì)濾波效果的影響。

1.2 最小劃分窗口

通常情況下，較大的最小劃分窗口有利于剔除離散的噪聲。考慮到礦區(qū)形變范圍小且形變量大，如果對(duì)整個(gè)解纏圖采用同樣的最小劃分窗口，則會(huì)導(dǎo)致礦區(qū)形變中心的有效觀測(cè)數(shù)據(jù)損失。采用粗差閾值為2倍中誤差，全局最小劃分窗口為64×64，對(duì)圖1進(jìn)行四叉樹濾波的效果如圖2(a-b)所示，其中圖2(a)為濾波后的結(jié)果，圖2(b)為濾除的點(diǎn)。可見在近場由于礦區(qū)形變梯度大，而最小劃分窗口過大，在濾波后的解纏圖中，礦區(qū)中心只留下小部分像素點(diǎn)(見圖2(a))，礦區(qū)形變中心的大量有效形變點(diǎn)被濾除(見圖2(b))。

當(dāng)最小劃分窗口較小時(shí)，各子區(qū)域劃分得比較小，各子區(qū)域內(nèi)的形變梯度相對(duì)較小，子區(qū)域內(nèi)像素值的中誤差會(huì)比大窗口時(shí)小，導(dǎo)致近場濾除的點(diǎn)比較少，但此時(shí)遠(yuǎn)場的噪聲也會(huì)隨著被劃分至不同子區(qū)域，導(dǎo)致各子區(qū)域的中誤差相對(duì)較小，濾波后仍然會(huì)存在大量的噪聲。同樣采用粗差閾值為2倍中誤差，圖2(c-d)顯示了全局最小劃分窗口縮小為8×8時(shí)對(duì)圖1進(jìn)行濾波的效果。由圖2(d)可見濾除的點(diǎn)很少，近場的形變點(diǎn)基本保留了，但由圖2(c)可見濾波后在遠(yuǎn)場仍有許多噪聲。

由上述比較可知，合理調(diào)節(jié)最小劃分窗口是獲得理想濾波效果的一個(gè)關(guān)鍵步驟。理論上近場的形變梯度比較大，應(yīng)采用較小的濾波窗口，而遠(yuǎn)場形變梯度小，則應(yīng)采用較大的濾波窗口。將礦區(qū)中心2 km范圍內(nèi)的區(qū)域作為近場形變中心，將2 km以外區(qū)域視為遠(yuǎn)場，同樣采用粗差閾值為2倍標(biāo)準(zhǔn)差，將近場的最小劃分窗口設(shè)置為8×8，而遠(yuǎn)場則設(shè)置為64×64，濾波效果如圖2(e-f)所示。通過對(duì)比濾波后的相位解纏圖(見圖2(e))與原始的相位解纏圖(見圖1) ，可以看出對(duì)近場與遠(yuǎn)場設(shè)置不同的最小劃分窗口既能有效剔除噪聲點(diǎn)，又不損失近場的形變數(shù)據(jù)。

1.3 粗差閾值

上節(jié)對(duì)最小窗口的分析中均采用2倍中誤差為粗差閾值，但在實(shí)際數(shù)據(jù)處理中，粗差閾值需根據(jù)影像噪聲的整體離散程度進(jìn)行設(shè)置，閾值過高會(huì)導(dǎo)致影像粗差點(diǎn)濾除不干凈，而閾值過低則會(huì)將大量非粗差點(diǎn)濾除。將最小劃分窗口均設(shè)置為遠(yuǎn)場64×64和近場8×8，圖2(e-f)、圖2(g-h)、圖2(i-j)分別表示粗差閾值為2倍、1.5倍與及3倍中誤差的濾波效果圖。對(duì)比發(fā)現(xiàn)，采用1.5倍中誤差作為粗差閾值雖然能將幾乎所有的噪聲點(diǎn)剔除，但是遠(yuǎn)場與近場的一些非噪聲點(diǎn)也被濾除(見圖2(g-h))；采用3倍中誤差作為粗差閾值雖然能保留形變中心的像素，但是濾波后遠(yuǎn)近場均存在較多的殘余噪聲點(diǎn)(見圖2(i-j)) ；而采用2倍中誤差(見圖2(e-f)) 則能在基本保留形變中心有效像素的基礎(chǔ)上，遠(yuǎn)場的濾波效果與采用1.5倍中誤差(見圖2(g-h)) 閾值濾波效果差異不大。

2 InSAR時(shí)間序列解算

2.1 基線選擇

本文選取45景2019年10月至2021年3月的Sentinel-1A影像數(shù)據(jù)。礦區(qū)為一露天礦區(qū)，礦區(qū)中心范圍幾乎無植被，且在選擇時(shí)間范圍內(nèi)空間基線均小于150 m，因此在組成干涉對(duì)時(shí)優(yōu)先考慮時(shí)間基線。由于礦區(qū)中心的形變速率快、范圍小，且Sentinel-1A為C波段的數(shù)據(jù)，波長比較短，如果時(shí)間基線過長，則會(huì)導(dǎo)致干涉圖條紋混疊，無法解纏，因此本文的所有時(shí)間基線都選取為12 d，即為相鄰的圖像組成干涉對(duì)，一共組成44個(gè)干涉對(duì)。

2.2 DInSAR數(shù)據(jù)處理

首先使用ISCE軟件[11]進(jìn)行DInSAR數(shù)據(jù)處理，干涉圖采用Goldstein濾波方法進(jìn)行濾波[12]，根據(jù)形變中心的干涉圖條紋密集情況以及遠(yuǎn)場條紋的清晰情況將Goldstein濾波的強(qiáng)度在0.3～0.8之間進(jìn)行調(diào)整，地形相位去除以及解纏結(jié)果的地理編碼都采用1角秒的美國SRTM DEM[13]。為使最終方位向與距離向地面分辨率一致，使用方位向與距離向?yàn)?:4對(duì)數(shù)據(jù)進(jìn)行多視處理。干涉圖的解纏方法采用枝切法[14]，并通過閉合環(huán)自動(dòng)檢測(cè)解纏誤差[15]。對(duì)由于去相干導(dǎo)致的孤島采用人工橋接相干性良好的區(qū)域進(jìn)行解纏[15]，對(duì)形變中心條紋密集區(qū)域進(jìn)行人工檢查解纏結(jié)果的正確性。

2.3 InSAR時(shí)間序列計(jì)算

在解纏完成后對(duì)解纏結(jié)果采用四叉樹濾波方法進(jìn)行濾波，其中近場和遠(yuǎn)場分別采用2×2和8×8的最小劃分窗口，粗差閾值為2倍中誤差。接著，將濾波后的解纏圖進(jìn)行地理編碼，采用 π ?RATE軟件計(jì)算礦區(qū)的時(shí)間序列結(jié)果[15]。時(shí)間序列計(jì)算中的軌道誤差采用網(wǎng)絡(luò)法的雙線性擬合模型改正。與地形相關(guān)的InSAR大氣延遲誤差改正與軌道誤差改正類似，采用線性擬合模型擬合干涉相位與地形的關(guān)系，將擬合后的大氣延遲誤差從解纏后相位中減掉。由于礦區(qū)范圍比較小，在遠(yuǎn)場對(duì)與地形不相關(guān)的大氣延遲誤差采用空間域低通濾波，獲得遠(yuǎn)場大氣延遲部分低通分量，然后將其進(jìn)行空間插值獲得礦區(qū)近場的大氣延遲誤差，然后將插值得到的與地形不相關(guān)的大氣延遲誤差減掉。在時(shí)域上，各個(gè)干涉圖的相干性不一致，且礦區(qū)形變的速度快，干涉相位正負(fù)交替，在時(shí)域上無規(guī)律，因此本文對(duì)InSAR時(shí)間序列不進(jìn)行時(shí)域?yàn)V波。

3 InSAR結(jié)果與分析

3.1 礦區(qū)形變時(shí)間序列結(jié)果

本文采用四叉樹濾波方法對(duì)相位解纏結(jié)果進(jìn)行濾波后，采用 π ?RATE軟件處理得到礦區(qū)最終的時(shí)間序列形變結(jié)果，如圖3所示。由圖可見從2019年10月至2021年3月，礦區(qū)地表形變劇烈，尤其是2020年6月至2020年9月時(shí)間范圍內(nèi)更為顯著。礦區(qū)形變?cè)跁r(shí)域上整體呈波浪狀，時(shí)沉降時(shí)抬升，無明顯規(guī)律，總體與礦區(qū)的開挖與修復(fù)工作進(jìn)度相關(guān)。以2019年11月5日為例，從圖3可清晰看到，在空間上，礦區(qū)形變顯著區(qū)域大致可分為3個(gè)區(qū)域，礦區(qū)中心區(qū)域既有抬升也有沉降。

3.2 礦區(qū)形變時(shí)空演化分析

礦區(qū)在研究時(shí)段內(nèi)發(fā)生了劇烈的地表形變，整個(gè)礦區(qū)間隔12 d的最大沉降增量為67.3 mm，發(fā)生在2020年8月31日，最大抬升增量為79.4 mm，發(fā)生在2019年11月5日。

本文選取3個(gè)形變中心區(qū)域作剖面線以進(jìn)行時(shí)域分析，剖面線位置如圖1所示，沿剖面線每30 m取一個(gè)點(diǎn)。圖4為3個(gè)形變中心A、B、C區(qū)域的剖面時(shí)空演化結(jié)果，其中上圖的折線表示每個(gè)時(shí)段該區(qū)域沿剖面所有點(diǎn)的平均形變量，下圖各點(diǎn)的顏色表示不同時(shí)段沿剖面所有點(diǎn)的形變值。由這3個(gè)剖面圖可見，3個(gè)形變中心區(qū)域的形變?cè)隽吭跁r(shí)域上呈波浪狀，時(shí)沉降時(shí)抬升，A區(qū)域和C區(qū)域在時(shí)域上均發(fā)生了嚴(yán)重的形變，形變起伏波動(dòng)大于B區(qū)域。A區(qū)域的形變范圍最大，長約2 km，寬約1 km。對(duì)比3個(gè)剖面圖的平均沉降量折線圖，在大多數(shù)情況下B區(qū)域與C區(qū)域的形變趨勢(shì)相同，而與A區(qū)域的趨勢(shì)相反。這是由于在礦區(qū)交替出現(xiàn)挖方和填方的施工，導(dǎo)致礦區(qū)抬升與沉降交替出現(xiàn)。

A剖面的最大抬升量為52.7 mm，最大平均抬升量為35.0 mm,發(fā)生在2019年11月5日，最大沉降量為36.1 mm，最大平均沉降量為27.0 mm，發(fā)生在2020年12月29日；B剖面的最大抬升量為20.3 mm，最大平均抬升量為10.2 mm，發(fā)生在2020年10月18日，最大沉降量為19.6 mm，最大平均沉降量為13.2 mm，發(fā)生在2020年8月31日；C剖面的最大抬升量為36.5 mm，最大平均抬升量為15.2 mm，發(fā)生在2020年10月18日,最大沉降量為53.3 mm，最大平均沉降量為31.7 mm，發(fā)生在2020年8月31日。

3.3 礦區(qū)形變結(jié)果驗(yàn)證

該礦區(qū)開采區(qū)北部布設(shè)了一些測(cè)量機(jī)器人，范圍如圖1的黑框所示。本文將32個(gè)點(diǎn)位41個(gè)時(shí)段測(cè)量機(jī)器人的測(cè)量結(jié)果轉(zhuǎn)換到InSAR視線(line-of-sight，LOS)方向，點(diǎn)位如圖5(a)中黑點(diǎn)所示。圖5(b)中彩色的圓圈代表這32個(gè)點(diǎn)所有時(shí)段的結(jié)果與InSAR時(shí)間序列形變結(jié)果的差值。圖5(c)～(e)為選取其中3個(gè)點(diǎn)P1、P2、P3(點(diǎn)位見圖5(a))的測(cè)量機(jī)器人與InSAR形變結(jié)果的對(duì)比，兩種手段的形變結(jié)果趨勢(shì)一致，大部分時(shí)段差異均小于20 mm。對(duì)所有點(diǎn)所有時(shí)段，兩種測(cè)量手段形變差值的平均值為0 mm，標(biāo)準(zhǔn)差為10.9 mm，這與測(cè)量機(jī)器人所能達(dá)到的測(cè)量精度一致，驗(yàn)證了本文InSAR時(shí)間序列形變結(jié)果的可靠性。

4 結(jié)論

本文將四叉樹濾波應(yīng)用于云南某礦區(qū)，并計(jì)算了該礦區(qū)的地表形變時(shí)間序列。研究結(jié)果表明，當(dāng)形變中心的形變梯度較大時(shí)，近場與遠(yuǎn)場采用不同最小劃分窗口的四叉樹濾波，既能保留近場的形變點(diǎn)，又能去除大量的噪聲，從而獲得干凈的形變結(jié)果。通過對(duì)該礦區(qū)形變時(shí)間序列的時(shí)空特征分析顯示，該礦區(qū)在2019年10月至2021年3月地表形變劇烈，形變?cè)跁r(shí)域上總體呈波浪狀，時(shí)沉降時(shí)抬升，在空間上，礦區(qū)形變顯著區(qū)域大致可分為3個(gè)區(qū)域，并且A區(qū)域和B、C區(qū)域形變趨勢(shì)相反，這是由于在礦區(qū)開采區(qū)與排土區(qū)同時(shí)存在挖方和填方，導(dǎo)致礦區(qū)抬升與沉降交替出現(xiàn)。整個(gè)礦區(qū)相隔12 d的最大沉降增量達(dá)到67.3 mm，最大抬升增量達(dá)到79.4 mm。通過對(duì)比測(cè)量機(jī)器人的形變監(jiān)測(cè)結(jié)果，驗(yàn)證了本文InSAR結(jié)果的可靠性。鑒于InSAR的技術(shù)優(yōu)勢(shì)和該礦區(qū)還有很多區(qū)域沒有布設(shè)地面監(jiān)測(cè)點(diǎn)，建議今后進(jìn)一步采用InSAR技術(shù)加強(qiáng)對(duì)礦區(qū)的監(jiān)測(cè)，確保礦區(qū)安全施工。