StaMPS技術(shù)在區(qū)域地表沉降形變監(jiān)測(cè)中的應(yīng)用

郭海京，鄭慶章，王 斌

（1.廣東省國(guó)土資源測(cè)繪院，廣東 廣州 510500；2.惠東縣國(guó)土資源信息中心，廣東 惠州 516300）

Sentinel-1A是歐洲空間局哥白尼計(jì)劃的首顆環(huán)境監(jiān)測(cè)衛(wèi)星，于2014-04-03發(fā)射升空，衛(wèi)星系統(tǒng)采用12 d 的重訪周期進(jìn)行不間斷全球成像服務(wù)，成像模式主要有4種：條帶模式（SM）、干涉寬幅模式（IW）、極寬幅模式（EW）和波譜模式（Wave）[1]，已在地震形變監(jiān)測(cè)及災(zāi)害應(yīng)急響應(yīng)、區(qū)域地面沉降監(jiān)測(cè)等領(lǐng)域擁有優(yōu)勢(shì)[2-4]，地面沉降是一種典型的地質(zhì)災(zāi)害[5]。當(dāng)自然或人為因素誘發(fā)地面垂直形變超過(guò)一定量時(shí)，會(huì)造成滑坡、崩塌、地貌塌陷、泥石流等多種災(zāi)害破壞。傳統(tǒng)地面沉降監(jiān)測(cè)技術(shù)如水準(zhǔn)測(cè)量、GPS高精度似大地水準(zhǔn)面面測(cè)量、傳感器網(wǎng)絡(luò)監(jiān)測(cè)等技術(shù)手段，點(diǎn)位測(cè)量精度較高[6-7]，但對(duì)于大范圍地面沉降監(jiān)測(cè)，需要耗費(fèi)大量人力、物力資源，強(qiáng)度大，精度不均勻。

InSAR技術(shù)相較傳統(tǒng)地表形變監(jiān)測(cè)手段，具有監(jiān)測(cè)范圍廣、監(jiān)測(cè)效率高、監(jiān)測(cè)網(wǎng)絡(luò)空間密度高、提供整體時(shí)序形變等特點(diǎn)[8]。劉一霖等[9]利用短基線集INSAR技術(shù)（SBAS）分析黃河三角洲地表形變特征場(chǎng)和時(shí)間序列特征，對(duì)地下水抽取等人為因素進(jìn)行了相關(guān)分析，但難以有效的獲取地表INSAR信息，對(duì)沉積物固結(jié)壓實(shí)產(chǎn)生的地面沉降因素未納入相關(guān)分析，開展大范圍沉降監(jiān)測(cè)和預(yù)警，陳炳乾等[10]基于D-INSAR以及永久散射體雷達(dá)干涉技術(shù)（PS-InSAR）獲取和分析開采區(qū)沉陷影響和發(fā)展趨勢(shì)不需要先驗(yàn)知識(shí)，選取PS點(diǎn)，獲取沉陷區(qū)影響范圍和發(fā)展趨勢(shì)，但也存在有植被的地區(qū)提取的PS點(diǎn)目標(biāo)形變結(jié)果不精準(zhǔn)，在未聯(lián)合GPS連續(xù)運(yùn)行參考站數(shù)據(jù)解算情況下，還難以全面掌握精細(xì)的地表三維形變場(chǎng)，時(shí)序INSAR技術(shù)特點(diǎn)如表1所示。

表1 時(shí)序INSAR技術(shù)特點(diǎn)比較

針對(duì)山地和丘陵地區(qū)，傳統(tǒng)的PSInSAR選點(diǎn)技術(shù)難以保證點(diǎn)的選擇精度及密度，StaMPS技術(shù)[11-12]是基于統(tǒng)計(jì)方法篩選PS點(diǎn)，同時(shí)StaMPS利用三維解纏算法，首先在時(shí)間上計(jì)算每個(gè)PS點(diǎn)的相位差異，然后設(shè)置參考點(diǎn)用最小二乘法在空間上進(jìn)行解纏，相比于傳統(tǒng)的空間解纏算法準(zhǔn)確性更高。對(duì)于空間相關(guān)誤差去除，利用相位解纏使PS點(diǎn)的相位值恢復(fù)到絕對(duì)值，但同時(shí)除了包含形變相位，也包含了大氣相位、衛(wèi)星軌道誤差、DEM殘差相位、隨機(jī)噪聲相位，而這些相位在時(shí)間上是不相關(guān)的，采用時(shí)間域上的高通濾波分離出形變相位，并用空間域上的高通濾波分離出輔影像的大氣相位和軌道誤差，在PS點(diǎn)較少的地區(qū)具有明顯優(yōu)勢(shì)。本文采用StaMPS方法，以Sentinel-1A為主要數(shù)據(jù)源，利用30 m分辨率DEM數(shù)據(jù)去除高程相位，形成2015年12月～2018年2月惠東縣及城區(qū)重點(diǎn)區(qū)域沉降結(jié)果。

1 研究目的與區(qū)域概況

惠東縣地處廣東省東南部、惠州東部，陸地總面積約3500km2，屬沿海山區(qū)縣。蓮花山脈分布在縣境北部和東北部，地勢(shì)較高，南部沿海多為丘陵，地勢(shì)較低，整個(gè)地勢(shì)是由東北向西南傾斜（如圖1）。西枝江自東北向西南流經(jīng)縣境。北部、東北部是山區(qū)，中部是沿海平原和丘陵，地勢(shì)較平緩，是惠東的主要產(chǎn)糧區(qū)和經(jīng)濟(jì)作物區(qū)，占全縣總面積36%。南部是稔平半島，屬沿海丘陵區(qū)，是惠東的漁、鹽生產(chǎn)基地，占全縣總面積21%。海岸線218.3km，海域面積3200km2。全年氣象災(zāi)害造成的損失較嚴(yán)重，主要是季風(fēng)低壓環(huán)流和熱帶氣旋影響，初臺(tái)出現(xiàn)時(shí)間偏早，暴雨強(qiáng)度大，覆蓋范圍廣，汛期強(qiáng)，對(duì)流天氣頻繁，局部對(duì)流強(qiáng)烈，時(shí)空分布不均，是地質(zhì)災(zāi)害監(jiān)測(cè)重點(diǎn)區(qū)域。

圖1 惠東縣區(qū)域圖

據(jù)統(tǒng)計(jì)，全縣現(xiàn)有地質(zhì)災(zāi)害隱患點(diǎn)69處，受威脅人數(shù)5990位，潛在經(jīng)濟(jì)損失1.45億元，其中受威脅人數(shù)100～1000人的重要地質(zhì)災(zāi)害隱患點(diǎn)3處，1000人以上特大型地質(zhì)災(zāi)害隱患點(diǎn)3處，受威脅人數(shù)4739人。2019年消除12處地質(zhì)災(zāi)害隱患點(diǎn)，解除受威脅人員100人，搬遷與治理比例為14.8%[13]。傳統(tǒng)實(shí)地勘察方式需要投入大量人力資源，費(fèi)事耗力，獲取的數(shù)據(jù)較為零散、現(xiàn)勢(shì)性不足，難以獲取地質(zhì)災(zāi)害點(diǎn)附近一定區(qū)域形變狀況，對(duì)全縣區(qū)域范圍地表沉降變化缺乏全域數(shù)據(jù)支持。鑒于此，本文采用時(shí)序StaMPS方法[14-15]， 研究實(shí)現(xiàn)全縣域地表沉降變化情況，特別是城區(qū)開發(fā)建設(shè)重點(diǎn)區(qū)域監(jiān)測(cè)，持續(xù)跟蹤驗(yàn)證2019年已解除的 8處重點(diǎn)災(zāi)害隱患點(diǎn)穩(wěn)定狀態(tài)，為動(dòng)態(tài)開展地質(zhì)環(huán)境時(shí)序變化監(jiān)測(cè)預(yù)警，地面沉降監(jiān)測(cè)網(wǎng)絡(luò)建設(shè)以及地質(zhì)災(zāi)害綜合治理規(guī)劃提供決策參考信息（如表2）。

表2 惠東縣部分地質(zhì)災(zāi)害隱患點(diǎn)及形變參考點(diǎn)

2 數(shù)據(jù)預(yù)處理

現(xiàn)有SAR傳感器主要采用L、X、C三個(gè)波段，所形成的SAR圖像有各自特點(diǎn)，對(duì)于同一區(qū)域的SAR圖像干涉處理，L波段的圖像相干性高于X、C波段，但高程信息的敏感度較X、C波段較弱。針對(duì)植被地段，波長(zhǎng)越長(zhǎng)，其穿透能力越強(qiáng)。X波段波長(zhǎng)較短，穿透能力弱，雷達(dá)接受的是植被覆蓋頂端反射的信息，而在城區(qū)監(jiān)測(cè)方面有不錯(cuò)的表現(xiàn)。C波段的波長(zhǎng)短，數(shù)據(jù)相干性低，但是對(duì)于同一形變區(qū)域，對(duì)于較小的形變，C波段的干涉條紋更為敏感，也更能展現(xiàn)其形變趨勢(shì)，3種波段中C波段性質(zhì)更為均衡[16-17]。針對(duì)惠東縣及重點(diǎn)區(qū)域地形地貌特征，本文采用C波段Sentinel-1A IW模式下的VH數(shù)據(jù)，幅寬250km，分辨率為5×20m，時(shí)間跨度為2015年12月至2018年2月。

Sentinel-1A IW模式下的TOPS[18]沿用ScanSAR成像模式在距離向采用多個(gè)子條帶分時(shí)段進(jìn)行掃描方式，每個(gè)burst在進(jìn)行掃描時(shí)，天線波束沿方位向從后向前旋轉(zhuǎn)掃描，壓縮了方位向波束方向圖并且抑制了成像時(shí)所造成的扇貝效應(yīng)。對(duì)于粗配準(zhǔn)選擇利用SAR成像幾何、軌道數(shù)據(jù)及外部DEM的幾何配準(zhǔn)方式，獲得主副影像之間的偏移量，精配準(zhǔn)采用增強(qiáng)譜分集的方式，利用響應(yīng)脈沖函數(shù)，對(duì)主輔影像方位向上兩個(gè)相鄰的burst的重疊區(qū)域分別做干涉處理，得到重疊區(qū)域的干涉相位差，獲得方位向配準(zhǔn)誤差，如圖2所示。

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

選取2016-12-25獲取的影像為主影像，形成57幅干涉對(duì)，并據(jù)配準(zhǔn)情況和干涉對(duì)質(zhì)量，挑選準(zhǔn)確配準(zhǔn)且相干性較高的干涉對(duì)作為最終的輸入值。利用最小生成樹算法篩選用于配準(zhǔn)的基線對(duì)，對(duì)選取后的基線對(duì)進(jìn)行干涉處理，利用經(jīng)地理編碼后的DEM去除高程相位，形成差分干涉對(duì)，并進(jìn)行多視及濾波處理，如圖3所示。

圖3 相關(guān)性高干涉對(duì)

3 時(shí)序StaMPS三維解算方法

StaMPS是在PS-InSAR基礎(chǔ)上發(fā)展的時(shí)序技術(shù)，主要原理是通過(guò)振幅離差法和相干系數(shù)法篩選出相位穩(wěn)定的PS，然后將PS點(diǎn)的相位重新采樣到格網(wǎng)上，利用形變相位、大氣延遲相位、軌道誤差相位在時(shí)間域和空間域上頻譜特性的不同，通過(guò)一定的濾波手段分離出形變相位和其他誤差相位，獲得精度較高的地表形變場(chǎng)，如圖4所示。

圖4 利用StaMPS地表反演總體技術(shù)流程

3.1 三維解纏

StaMPS技術(shù)能降低對(duì)時(shí)間和空間基線的要求，進(jìn)行大量的SAR數(shù)據(jù)分析，當(dāng)研究區(qū)內(nèi)SAR圖像豐富且能選取足夠密度的PS點(diǎn)時(shí)，可測(cè)得0.1mm/y的視向移動(dòng)速度，同時(shí)降低了失相干性影響，有效的提取大氣、軌道和 DEM 等誤差因素的影響[10]。StaMPS中三維解纏分別體現(xiàn)在空間和時(shí)間上，基于最小化L-P范框架，將相位解纏看作最優(yōu)化問(wèn)題，找到目標(biāo)函數(shù)的最小值，從而獲取解纏后的相位值。目標(biāo)函數(shù)為：

式中，Δφ為解纏相位梯度；Δψ為纏繞相位梯度；x和y分別表示二維的兩個(gè)方向；z表示第三維方向；i和j表示點(diǎn)坐標(biāo)；w表示權(quán)重。和傳統(tǒng)L-P范分布不同之處在于，Hooper在StaMPS加入了第三維度的相位梯度信息，通常為時(shí)間序列信息。這里，參數(shù)P決定如何處理Δφ和Δψ的差值關(guān)系，當(dāng)P=0時(shí)，若僅考慮二維信息，上式為類似枝切解纏算法的L-0范數(shù)目標(biāo)函數(shù)；當(dāng)P=2時(shí)，即為最小二乘解纏算法。

3.2 最小二乘解算

1）觀測(cè)相位信息為：

2）和空間相關(guān)的DEM誤差項(xiàng)可表示為：

式中，Bi⊥為第i幅干涉圖垂直基線；Ri M為第i幅干涉圖主影像和目標(biāo)距離；θi為主影像成像時(shí)雷達(dá)和目標(biāo)傾角；Δh為高程估計(jì)誤差。

3）設(shè)第i幅干涉圖的形變量為ΔRi，RM為主影像傳感器與目標(biāo)的斜距，RS為副影像傳感器與目標(biāo)的斜距，Ci為主影像到副影像時(shí)間序列中SAR影像數(shù)，vk為形變速率，則有：

4）假設(shè)線性形變速率v，則在整個(gè)時(shí)間序列中有，所以線性形變相位可以表達(dá)為：

4 結(jié)果與分析

實(shí)驗(yàn)共選取了487 997個(gè)PS點(diǎn)，形成2015年12月至2018年2月期間全縣域地表形變速率，并估算其累計(jì)變化量。圖5～6可以看出惠東縣最大形變集中在縣東部山體；城區(qū)及沿海部分有少量抬升，區(qū)域地質(zhì)情況較為穩(wěn)定。

其中，圖5為2015年12月至2018年2月期間得到的平均形變速率圖，惠東縣東部山體部分有較明顯沉降，最大形變速率達(dá)到-10mm/y。城區(qū)部分及沿海部分有少量抬升，速率最大達(dá)到4mm/y。2015年 12月至2018年2月得到的累積時(shí)序形變序列，其中最大累積沉降量達(dá)-24 mm。

考慮城區(qū)經(jīng)濟(jì)建設(shè)開發(fā)重點(diǎn)形變監(jiān)測(cè)區(qū)，并持續(xù)跟蹤監(jiān)測(cè)驗(yàn)證2019年已解除的其中8處重點(diǎn)災(zāi)害隱患點(diǎn)，對(duì)惠東縣重點(diǎn)區(qū)域及災(zāi)害隱患點(diǎn)進(jìn)行解算，從 圖5～7可見，區(qū)域形變基本在0～-2 mm/y。因StaMPS選點(diǎn)的特點(diǎn)，會(huì)有少量的噪聲點(diǎn)被選為PS點(diǎn)，因此會(huì)有-8 mm/y的形變出現(xiàn)，該結(jié)果作為噪聲對(duì)待。從圖8形變序列中可以看到，該區(qū)域地表較穩(wěn)定，累積形變量大部分在0～-4 mm之間，已解除的其中8處重點(diǎn)災(zāi)害隱患點(diǎn)地質(zhì)區(qū)域地質(zhì)情況較為穩(wěn)定。

圖5 惠東縣平均形變速率

圖6 惠東縣重點(diǎn)區(qū)域平均形變速率

圖8 惠東縣重點(diǎn)區(qū)域形變序列（部分）

5 結(jié) 論

1）開展StaMPS反演技術(shù)對(duì)研究區(qū)全域特別是重點(diǎn)災(zāi)害危險(xiǎn)區(qū)域進(jìn)行地表形變監(jiān)測(cè)，可克服傳統(tǒng)InSAR技術(shù)的誤差影響，提高SAR影像的使用率和監(jiān)測(cè)精度，能夠較為快速的發(fā)現(xiàn)顯著變形區(qū)，有針對(duì)性的提前做好地質(zhì)災(zāi)害預(yù)警與防治[19-20]。本次研究實(shí)驗(yàn)發(fā)現(xiàn)惠東縣近海域有部分地區(qū)出現(xiàn)輕微的抬升，但總體不超過(guò)3 mm/y；C波段數(shù)據(jù)在多植被覆蓋區(qū)相干性較差，東部選點(diǎn)較少，城區(qū)部分地形較穩(wěn)定，基本無(wú)形變趨勢(shì)，山區(qū)有部分地區(qū)存在沉降。重點(diǎn)對(duì)已解除的8處重點(diǎn)災(zāi)害隱患點(diǎn)開展持續(xù)跟蹤監(jiān)測(cè)，獲得其時(shí)序沉降信息，總體形變速率在0～-2 mm/y，說(shuō)明其地質(zhì)情況較為穩(wěn)定。

2）本次研究主要是利用StaMPS技術(shù)，該技術(shù)在選點(diǎn)和解纏方面優(yōu)勢(shì)較大。但是StaMPS技術(shù)是以線性模型求解形變相位和沉降速率，若實(shí)際形變并不是以線性模型沉降，這會(huì)導(dǎo)致StaMPS技術(shù)得出的沉降值偏小，對(duì)滑坡體等突發(fā)自然災(zāi)害形變監(jiān)測(cè)還難以滿足要求。建議在選點(diǎn)和解纏方面可以利用該技術(shù)，后續(xù)的解算可以考慮與其他時(shí)序技術(shù)進(jìn)行結(jié)合處理，以提高在區(qū)域性地表形變監(jiān)測(cè)成效。

3）本次研究使用的數(shù)據(jù)源主要是C波段的Sentinel-1A數(shù)據(jù)，SAR衛(wèi)星影像較為單一，沒有使用不同衛(wèi)星對(duì)具有不同地表特征的區(qū)域分別進(jìn)行形變監(jiān)測(cè)分析。Sentinel-1A數(shù)據(jù)介于X波段和L波段之間，雖然均衡性較好，但是在城區(qū)監(jiān)測(cè)效果方面不如X波段數(shù)據(jù)；在高植被覆蓋區(qū)，C波段數(shù)據(jù)失相干嚴(yán)重，造成了選點(diǎn)較少，監(jiān)測(cè)效果不佳。建議對(duì)城區(qū)采用 X波段數(shù)據(jù)，對(duì)高植被區(qū)采用L波段數(shù)據(jù)能更加有效精確地反演地表形變特征。

圖7 惠東縣形變序列（部分）