羅 莉 ， 王 斌

（1. 廣東省測繪工程公司， 廣州 510670； 2. 廣東省國土資源測繪院， 廣州 510500）

0 引言

地面沉降是一種典型的地質(zhì)災(zāi)害[1]。 當(dāng)自然或人為因素誘發(fā)地面垂直形變超過一定量時， 會造成滑坡、 崩塌、 地貌塌陷、 泥石流等多種災(zāi)害破壞[2]。 傳統(tǒng)地面沉降監(jiān)測技術(shù)如水準(zhǔn)測量、 傳感器網(wǎng)絡(luò)監(jiān)測[3]等技術(shù)手段， 點位測量精度優(yōu)于毫米級， GPS 高精度似大地水準(zhǔn)測量精度達到或者優(yōu)于厘米級[4-5]， 在建筑變形測量、 邊坡監(jiān)測、 區(qū)域地表沉降[6]等領(lǐng)域發(fā)揮重要作用。 對于大范圍地面沉降監(jiān)測， 需耗費大量人力、 物力資源， 勞動強度大， 精度不均勻。 InSAR 是一種新興的空間大地測量技術(shù)[7]， 目前已被廣泛用于地表沉降監(jiān)測， 自然災(zāi)害 （如地震、 火山、 滑坡、 水土流失等）監(jiān)測及環(huán)境監(jiān)測（如海上溢油、 土地利用變化等）等民用、 軍用領(lǐng)域。 相較傳統(tǒng)監(jiān)測手段， 具有監(jiān)測范圍廣、 監(jiān)測效率高、 監(jiān)測網(wǎng)絡(luò)空間密度高、 提供整體時序形變等特點[8]。 劉一霖等利用短基線集INSAR 技術(shù)（SBAS）分析黃河三角洲地表形變特征場和時間序列特征， 對地下水抽取等人為因素進行了相關(guān)分析[9]， 但難以獲取有效地表INSAR 信息， 對沉積物固結(jié)壓實產(chǎn)生的地面沉降因素未納入相關(guān)分析。 開展大范圍沉降監(jiān)測和預(yù)警， 陳炳乾等基于D-INSAR 和永久散射體雷達干涉技術(shù)（PS-INSAR）獲取和分析開采區(qū)沉陷影響和發(fā)展趨勢[10-12]， 不需要先驗知識獲取沉陷區(qū)影響范圍和發(fā)展趨勢， 但也存在在有植被的地區(qū)提取的PS 點目標(biāo)形變結(jié)果可靠不足， 在未聯(lián)合GPS 連續(xù)運行參考站數(shù)據(jù)解算情況下， 還難以全面掌握精細(xì)的地表三維形變場。

本文采用StaMPS-PS（Stanford Method for Persistent Scatters PS-INSAR）技術(shù)[13-14]， 分析研究區(qū)地表沉降變化。 通過振幅離差法和相干系數(shù)法篩選出相位穩(wěn)定的PS， 然后將PS 點的相位重新采樣到格網(wǎng)上， 利用形變相位、 大氣延遲相位、 軌道誤差相位在時間域和空間域上頻譜特性的不同， 通過一定的濾波手段分離出形變相位和其他誤差相位，獲得精度較高的地表形變場。

1 研究區(qū)與數(shù)據(jù)預(yù)處理

1.1 研究區(qū)域

惠州市地勢總體北高南低， 地貌類型復(fù)雜多樣， 有山地、 丘陵、 臺地、 平原和湖泊， 北部多為山地和高丘陵， 南部則為平原和臺地。 構(gòu)成惠州市各類地貌的基巖巖石以花崗巖最為普遍， 砂巖和變質(zhì)巖也較多， 此外局部還有紅色巖系地貌，沿海沿河地區(qū)多優(yōu)質(zhì)沙灘及復(fù)雜形態(tài)的珊瑚礁， 地質(zhì)構(gòu)造多為第四紀(jì)沉積層。 這些不同的地貌類型及地質(zhì)構(gòu)造導(dǎo)致惠州市地表沉降變化復(fù)雜， 尤其是連續(xù)性地表沉降可以發(fā)展為一種嚴(yán)重的地質(zhì)災(zāi)害， 能持續(xù)對道路、 房屋建筑、 地下管線等造成破壞。 尤其位于珠江三角洲平原區(qū)的惠州市是廣東省經(jīng)濟發(fā)達地區(qū)之一， 人口稠密， 經(jīng)濟繁榮， 但近幾年出現(xiàn)的地面沉降給當(dāng)?shù)毓贰?橋梁、 水利設(shè)施、 地下管網(wǎng)設(shè)施等造成了不同程度的破壞。 如圖1 所示。

表1 Sentinal-1A 主要參數(shù)Table 1 Main parameters of Sentinal-1A

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

1.2.1 SAR 影像數(shù)據(jù)

歐空局Sentinel-1A 作為歐洲空間局哥白尼計劃的首顆環(huán)境監(jiān)測衛(wèi)星， 于2014 年4 月3 日發(fā)射升空。 在 2015 年 4 月至 5 月期間， 衛(wèi)星系統(tǒng)穩(wěn)定運行后開始采用12 d 的重訪周期進行不間斷全球成像任務(wù)。 成像模式主要有4 種： 條帶模式（Strip Map， SM）、 干涉寬幅模式（Interferometric Wide，IW）、 極寬幅模式（Extra-Wide Swath， EW）和波譜模式（Wave）。 本項目采用 IW 模式下的 VH 數(shù)據(jù)， 具體參數(shù)如表1。

選取合適的SLC 影像作為原始SAR 影像。 IW模式下 Sentinel-1A 采用 TOPS 技術(shù)， 一景 SAR 影像可以成像三個子條帶及多個Burst， 根據(jù)惠州市范圍對原始影像進行裁剪及拼接， 形成最終符合要求的整幅單視復(fù)數(shù)影像。 如圖2 所示。

圖1 惠州市轄區(qū)分布圖Fig.1 Distribution map of Huizhou City

圖2 惠州市單視復(fù)數(shù)影像Fig.2 The single-view complex image of Huizhou City

1.2.2 DEM 數(shù)據(jù)

采用 ALOS Global Digital Surface Model， 30 m分辨 率 DEM 高程 數(shù) 據(jù)。 AW3D30（ALOS World 3D－30 m）是日本宇宙航空研究開發(fā)機構(gòu)（JAXA）于2016 年5 月發(fā)布的全球數(shù)字表面模型數(shù)據(jù)集， 水平分辨率約30 m。 由陸地觀測衛(wèi)星ALOS 上搭載的PRISM 立體相機對全球陸地±80°緯度地區(qū)進行覆蓋觀測。

圖3 惠州市地理編碼后DEMFig.3 The geocoding DEM of Huizhou City

AW3D30 存在一定范圍內(nèi)的無效數(shù)據(jù)值， 在處理過程中， 對原始的DEM 進行缺失數(shù)據(jù)的探測， 沿海范圍有少量的數(shù)據(jù)缺失， 將其高程置為0值。 然后對預(yù)處理后的DEM 數(shù)據(jù)進行拼接及地理編碼， 使其與SAR 影像統(tǒng)一坐標(biāo)系。 如圖3 所示。

2 時序StaMPS 模型及處理方法

2.1 時序StaMPS

StaMPS 技術(shù)是基于統(tǒng)計方法篩選PS 點， 利用三維解纏算法， 首先在時間上計算每個PS 點的相位差異， 然后設(shè)置參考點用最小二乘法在空間上進行解纏。 對于空間相關(guān)誤差去除， 利用相位解纏使PS 點的相位值恢復(fù)到絕對值， 但同時除了包含形變相位， 也包含了大氣相位、 衛(wèi)星軌道誤差、DEM 殘差相位、 隨機噪聲相位， 而這些相位在時間上是不相關(guān)的， 采用時間域上的高通濾波分離出形變相位， 并用空間域上的高通濾波分離出輔影像的大氣相位和軌道誤差。 其具體算法流程如下：

（1）在N 幅影像中選取1 幅作為主影像， 其余影像均與主影像配準(zhǔn)， 并進行輻射定標(biāo)， 將副影像的幾何信息和幅度信息與主影像匹配。

（2）作干涉及差分處理， 差分主要去除平地相位， 并利用外部DEM 去除地形相位， 得到N-1 幅干涉圖和差分干涉圖。

（3）利用 M 幅影像的幅度信息， 設(shè)置閾值， 進行PS 點的選取。

（4）進行相位解纏， PSInSAR 較為常用的空間域解纏方法為稀疏格網(wǎng)解纏。 估計相鄰PS 點的纏繞相位梯度， 對梯度進行積分， 空間搜索， 進行解纏。

（5）計算相應(yīng)的形變量、 DEM 誤差， 并從殘差相位中分離大氣擾動和非線性形變量。

（6）得到各個分量之后， 進行后處理工作， 如克里金插值， 最終得到整個實驗區(qū)的形變場。

2.2 三維解纏及最小二乘解算

StaMPS 中三維解纏分別體現(xiàn)在空間（二維）和時間（一維）上， 基于最小化L-P 范框架， 將相位解纏看做最優(yōu)化問題， 找到目標(biāo)函數(shù)的最小值，從而獲取解纏后的相位值。 目標(biāo)函數(shù)為：

圖4 StaMPS 處理流程Fig.4 The processing flow of StaMPS

式（1）中， ΔΦ 為解纏相位梯度， Δψ 為纏繞相位梯度， x 和y 分別表示二維的兩個方向， z 表示第三維方向， i 和 j 表示點坐標(biāo)， w 表示權(quán)重。 和傳統(tǒng)L-P 范分布不同之處在于， Hooper 在StaMPS加入了第三維度的相位梯度信息， 通常為時間序列信息。 這里， 參數(shù)P 決定如何處理ΔΦ 和Δψ 的差值關(guān)系， 當(dāng)P=0 時， 若僅考慮二維信息， 上式為類似枝切解纏算法的L-0 范數(shù)目標(biāo)函數(shù)； 當(dāng)P=2 時，即為最小二乘解纏算法。 在解纏完成并去除空間不相關(guān)的誤差后， 構(gòu)建觀測方程， 加入最小二乘準(zhǔn)則， 求取線性形、 主影像的大氣延遲和軌道誤差以及空間相關(guān)的 DEM 誤差

參數(shù)估值及協(xié)方差陣可以解得：

3 試驗結(jié)果及分析

3.1 平均形變速率

惠州市覆蓋范圍跨度較大， 若在不做多視的情況下整體處理， 對內(nèi)存的要求很高且時間效率低。 本文對原始影像裁剪分為龍門縣、 博羅縣、惠城、 惠陽區(qū)和惠東縣分別進行處理。 每個區(qū)域的控制點選取穩(wěn)定的水準(zhǔn)點作為參考， 惠州市整體形變速率圖根據(jù)各區(qū)域的形變結(jié)果進行拼接，選取了龍門 277 092、 博羅 974 064、 惠城 794 967、 惠陽 1 001 252、 惠東 487 997 個 PS 點， 獲取其形變結(jié)果。 試驗發(fā)現(xiàn)： 2015 年 12 月至2018年2 月期間， 惠州市地表平均形變速率整體形變較小， 大部分區(qū)域在-3～3 mm/y 之間， 龍門縣南部地區(qū)和博羅縣城區(qū)有部分下沉（紅圈范圍），最大形變速率達到-12 mm/y （如圖 5 所示）。 惠城、 惠陽區(qū)及惠東縣近海域有部分地區(qū)出現(xiàn)輕微的抬升，但不超過3 mm/y （藍(lán)圈范圍）。 由于C 波段數(shù)據(jù)在多植被覆蓋區(qū)相干性較差， 惠東縣東部選點較少，城區(qū)部分地形較穩(wěn)定， 基本無形變趨勢， 山區(qū)有部分地區(qū)存在沉降。 博羅縣和惠陽區(qū)穩(wěn)定性較差，相對沉降速率較高； 惠東縣近海域地區(qū)以及惠城區(qū)有少量抬升情況， 如圖6 所示。

圖5 惠州區(qū)域平均形變速率Fig.5 The average deformation rate of Huizhou

3.2 形變序列

2015 年 12 月至 2018 年 2 月期間， 惠州市地表整體形變趨勢隨時間變化總體較穩(wěn)定。 累積沉降量 （年平均形變速率與監(jiān)測期乘積）最大為-24 mm， 抬升量最大為 8 mm， 如圖 7 所示。

3.3 結(jié)果討論與分析

圖6 惠州市形變速率Fig.6 The deformation rate of Huizhou

（1）利用StaMPS 技術(shù)對惠州市地表沉降進行分區(qū)監(jiān)測， 獲得了沉降速率、 時序累計沉降值。根據(jù)監(jiān)測結(jié)果， 龍門縣南部地區(qū)和博羅縣城區(qū)有部分下沉情況， 最大形變速率達到-12 mm/y， 應(yīng)該予以重視。 博羅縣最大累積形變量量達-24 mm、惠陽區(qū)北部萬里工業(yè)區(qū)附件形變速率達到-6 mm/y，穩(wěn)定性較差， 相對沉降速率較高， 對于這些地區(qū)的地表狀況應(yīng)予以關(guān)注； 惠東縣近海域地區(qū)以及惠城區(qū)速率最大達到4 mm/y 有少量抬升情況。

（2）重點監(jiān)測惠州市整體地表時序沉降， 并分析整個目標(biāo)區(qū)內(nèi)的沉降時空特征及變化， 對惠州市地表沉降做出預(yù)判， 反映地表形變的動態(tài)變化信息， 為城市地表沉降防治提供依據(jù)。 為惠州市經(jīng)濟社會發(fā)展提供科學(xué)可靠的地理國情信息服務(wù)。

（3）本項目使用的SAR 衛(wèi)星影像較為單一， 主要是Sentinel-1A 的C 波段數(shù)據(jù)， 沒有使用不同衛(wèi)星對具有不同地表特征的區(qū)域分別進行形變監(jiān)測分析比較。 Sentinel-1A 數(shù)據(jù)介于X 波段和L 波段之間， 雖然均衡性較好， 但是在城區(qū)監(jiān)測方面不如X 波段數(shù)據(jù)； 另一方面， 惠州市龍門縣和惠東縣地貌多變， 植被覆蓋范圍較廣， C 波段數(shù)據(jù)失相干嚴(yán)重， 造成了選點過少因此無法有效監(jiān)測的問題。 建議對城區(qū)采用X 波段數(shù)據(jù)， 對山區(qū)采用L波段數(shù)據(jù)能更加有效精確地反演地表形變特征。

（4）本項目主要使用的StaMPS 的技術(shù)， 該技術(shù)在選點和解纏方面優(yōu)勢較大， 但是StaMPS 技術(shù)是以線性模型求解形變相位和沉降速率， 若實際形變并不是以線性模型沉降， 這會導(dǎo)致StaMPS 技術(shù)得出的沉降值偏小。 建議在選點和解纏方面可以利用該技術(shù)， 后續(xù)的解算可以考慮與其他時序技術(shù)聯(lián)合處理。

圖7 惠州市時序形變速率圖（部分）Fig.7 The time series deformation rate diagram of Huizhou City （part）