鄧曉景，曲國(guó)慶，張建霞，席換，王暉

(山東理工大學(xué) 建筑工程學(xué)院, 山東 淄博 255049)

地面沉降是指由于自然因素和人類活動(dòng)因素共同作用下引起地面土層壓縮，從而發(fā)生的地質(zhì)災(zāi)害[1]。傳統(tǒng)地面沉降監(jiān)測(cè)方法，如水準(zhǔn)測(cè)量、GPS測(cè)量等，受設(shè)備、外部環(huán)境等因素影響，存在監(jiān)測(cè)時(shí)間長(zhǎng)、監(jiān)測(cè)范圍小、維護(hù)成本高等缺點(diǎn)，難以實(shí)現(xiàn)連續(xù)、準(zhǔn)確的地面沉降監(jiān)測(cè)。隨著雷達(dá)衛(wèi)星監(jiān)測(cè)技術(shù)的不斷發(fā)展，合成孔徑雷達(dá)差分干涉測(cè)量技術(shù)(Differential Interferometric SAR, D-InSAR)已經(jīng)成為一種全天時(shí)、全天候、高精度監(jiān)測(cè)地面沉降的有效方法[2]。但D-InSAR技術(shù)易受空間失相干、時(shí)間失相干、大氣延遲等因素影響，使其監(jiān)測(cè)精度降低，無(wú)法進(jìn)行長(zhǎng)時(shí)間序列的地面沉降監(jiān)測(cè)[3]。永久散射體干涉測(cè)量技術(shù)(Permanent Scatterer InSAR,PS-InSAR)通過(guò)獲取高相干性目標(biāo)點(diǎn)，采用時(shí)空濾波的方法，能有效的削弱基線誤差、大氣延遲誤差等因素的影響，得到高精度時(shí)序地面沉降監(jiān)測(cè)結(jié)果[4]。

在利用PS-InSAR技術(shù)進(jìn)行大區(qū)域地面沉降監(jiān)測(cè)時(shí)，由于單一軌道覆蓋范圍有限，無(wú)法滿足研究區(qū)需要，往往需要將升軌和降軌沉降監(jiān)測(cè)結(jié)果進(jìn)行拼接。由于在同步觀測(cè)時(shí)間下，相同研究區(qū)域內(nèi)多種觀測(cè)模式下獲取的監(jiān)測(cè)結(jié)果是對(duì)相同地面點(diǎn)的多角度采樣，這為升降軌兩組解算結(jié)果的精度互檢提供了條件[5]。由于監(jiān)測(cè)結(jié)果中垂直方向分量受位移影響較小，將雷達(dá)視線向形變量轉(zhuǎn)換為垂直向形變量后再進(jìn)行精度互檢，這樣可以提高監(jiān)測(cè)結(jié)果內(nèi)部互檢的準(zhǔn)確性。

升降軌PS-InSAR數(shù)據(jù)融合的基本條件是提取同名相干目標(biāo)的形變量。然而，受雷達(dá)采樣方式差異的影響，升降軌數(shù)據(jù)中PS點(diǎn)的分布位置不相同，導(dǎo)致在升降軌解算結(jié)果中，同名相干目標(biāo)點(diǎn)數(shù)量極為有限，嚴(yán)重影響升降軌解算結(jié)果融合的可靠性[6]。針對(duì)上述問(wèn)題，本文提出了利用抗差最小二乘擬合法對(duì)升降軌PS-InSAR解算結(jié)果進(jìn)行精度內(nèi)部互檢及數(shù)據(jù)融合的方法。對(duì)升降軌數(shù)據(jù)進(jìn)行融合，可以加密單軌條件下的觀測(cè)序列，能夠更加豐富地顯示研究區(qū)地面沉降監(jiān)測(cè)時(shí)序變化規(guī)律。

1 研究區(qū)概況及數(shù)據(jù)選取

研究區(qū)為山東省東營(yíng)市，其大部分位于黃河三角洲沖積平原上，地勢(shì)平坦，地形以平原為主[7]。該地區(qū)海洋資源豐富，沿海地區(qū)分布大量曬鹽和養(yǎng)殖場(chǎng)，且石油、天然氣等自然資源豐富，其地層覆蓋主要為第四系，巖性由粉砂、細(xì)砂、黃色黏土和沙礫層組成[8]。東營(yíng)市近年來(lái)由于礦產(chǎn)資源開采以及人類活動(dòng)等因素影響，開始出現(xiàn)地面沉降現(xiàn)象。

圖1 研究區(qū)概況及數(shù)據(jù)范圍Fig.1 Overview of study area and data range

哨兵1號(hào)(Sentinel-1)衛(wèi)星是歐空局哥白尼計(jì)劃(GMES)中的地球觀測(cè)衛(wèi)星，發(fā)射于2014年，由兩顆衛(wèi)星組成，載有C波段合成孔徑雷達(dá)，本研究選擇覆蓋東營(yíng)市的Sentinel_1A雷達(dá)衛(wèi)星數(shù)據(jù)，該數(shù)據(jù)采樣模式為IW寬幅干涉模式，在該模式下采用TOPS (terrain observation with progressive scans)掃描方式，可以同時(shí)獲取3個(gè)子條帶，其幅寬為250 km,重訪周期為12 d，空間分辨率為5 m×20 m[9]。實(shí)驗(yàn)數(shù)據(jù)選取分別為降軌數(shù)據(jù)(Track-76)左側(cè)子條帶和升軌數(shù)據(jù)(Track-69)中間子條帶，數(shù)據(jù)參數(shù)見(jiàn)表1。參考DEM數(shù)據(jù)為美國(guó)宇航局提供的SRTM(shuttle radar topography mission)數(shù)據(jù)，其分辨率為90 m。

表1 Sentinel-1數(shù)據(jù)參數(shù)Tab.1 Sentinel-1 data parameters

2 研究方法

2.1 PS-InSAR技術(shù)

PS-InSAR技術(shù)將數(shù)據(jù)處理分為數(shù)據(jù)預(yù)處理、差分干涉計(jì)算和形變量計(jì)算三個(gè)步驟。首先，選取覆蓋研究區(qū)的N+1幅SAR影像，根據(jù)基線參數(shù)和多普勒效應(yīng)，選取一幅SAR影像作為公共主影像，其余影像為輔影像，分別與公共主影像配準(zhǔn)進(jìn)行干涉處理，生成N幅干涉圖；通過(guò)精密軌道文件去除平地效應(yīng)，同時(shí)利用參考DEM數(shù)據(jù)去除地形相位，對(duì)N幅干涉圖進(jìn)行差分處理，生成差分干涉圖[10]。其中，任意差分干涉圖像元的相位Δφ可以表示為：

Δφ=φdef+φtopo+φatm+φnoise

(1)

式中：φdef為由于地形位移而產(chǎn)生的形變相位,φtopo為外部DEM引起的殘差地形相位,φatm為大氣因素引起的延遲相位,φnoise為其他因素引起的噪聲相位。

然后，通過(guò)設(shè)定合理振幅離差閾值，提取PS候選點(diǎn)，采用Delaunay三角網(wǎng)法構(gòu)建PS網(wǎng)絡(luò)；根據(jù)建立相鄰PS點(diǎn)相位差分模型，獲得地表形變信息；利用空間搜索法，將PS點(diǎn)中線性形變和高程誤差去除[11]；通過(guò)時(shí)間域上高通濾波，空間域上低通濾波，獲取并分離非線性形變相位和大氣延遲相位，最終得到研究區(qū)時(shí)序地面形變量[12]。

2.2 升降軌數(shù)據(jù)融合簡(jiǎn)化模型

PS-InSAR技術(shù)獲取的地表形變量可分為東西方向和垂直方向的地表形變量，對(duì)于同一研究區(qū)的升降軌數(shù)據(jù)，PS-InSAR技術(shù)獲取的觀測(cè)結(jié)果東西方向上水平形變量呈現(xiàn)相反的特點(diǎn)，而垂直方向上垂直形變量基本相同[13]。雷達(dá)視線向(line of sight,LOS)累計(jì)形變量dLOS為：

dLOS=dNsinφsinθ-dEcosφsinθ+dVcosθ

(2)

式中:dN為南北方向的形變分量,dE為東西方向的形變分量,dV為垂直方向的形變分量,φ為方位角,θ為入射角。由于南北方向的水平位移對(duì)垂直向和東西向的影響非常小[5]，所以雷達(dá)視線向累計(jì)形變量簡(jiǎn)化為：

dLOS≈dEsinθ+dVcosθ

(3)

地表形變量以垂直方向的垂直形變量為主，東西方向的水平形變很小，即dE?dLOS,因此，垂直方向的垂直形變量可近似表示為：

dV≈dLOS/cosθ

(4)

針對(duì)升降軌數(shù)據(jù)融合中同名目標(biāo)點(diǎn)數(shù)量稀少問(wèn)題，傳統(tǒng)方法采用選取一條軌道解算結(jié)果作為主軌道，另一條軌道的解算結(jié)果為輔軌道，采用插值方法提取升降軌道相同位置的形變量，根據(jù)公式(5)進(jìn)行升降軌道解算結(jié)果的數(shù)據(jù)融合[14]。

(5)

在實(shí)際數(shù)據(jù)處理過(guò)程中，由于受到升降軌數(shù)據(jù)入射角以及采樣方式的影響，即使相同位置的垂直形變量偏差也不盡相同，只用主副軌道垂直形變量的平均偏差進(jìn)行整體修正，會(huì)嚴(yán)重影響修正后的輔軌道垂直形變量的可靠性。因此，通過(guò)對(duì)輔軌道進(jìn)行插值處理，獲得主副軌道相同位置的整體偏差，再利用抗差最小二乘方法對(duì)整體偏差進(jìn)行曲面擬合，用偏差擬合曲面矯正輔軌道，從而進(jìn)行整體形變場(chǎng)的解算。

2.3 抗差最小二乘曲面擬合

根據(jù)研究區(qū)內(nèi)升降軌解算結(jié)果的覆蓋情況，將降軌作為主軌道，升軌為輔軌道，以升降軌共同覆蓋區(qū)域?yàn)檠芯繉?duì)象，對(duì)升軌解算結(jié)果進(jìn)行插值處理，獲得降軌相干目標(biāo)點(diǎn)同名位置的偏差，采用二次曲面方程進(jìn)行偏差擬合,設(shè)二次曲面函數(shù)為：

Z=a0+a1x+a2y+a3xy+a4x2+a5y2

(6)

式中：Z為已知降軌相干目標(biāo)點(diǎn)偏差值；x，y為已知降軌相干目標(biāo)點(diǎn)的坐標(biāo)值；a1，a2，a3，a4，a5，a6為待求系數(shù)。

其誤差方程矩陣形式為:

V=KM-Z

(7)

式中:V為觀測(cè)值的殘差向量，K為待求系數(shù)矩陣；M為設(shè)計(jì)矩陣，Z為觀測(cè)向量。按照最小二乘原理,VTPV=min的原則，求解得：

K=(MTPM)-1MTPZ

(8)

當(dāng)偏差數(shù)據(jù)含有粗差時(shí),如果使用初始權(quán)值,會(huì)使擬合曲面差生彎曲，影響擬合結(jié)果精度[15]。利用抗差估計(jì)理論，構(gòu)造極值函數(shù)為：

(9)

式中:Pi為第i個(gè)觀測(cè)值的權(quán)重，ρ為凸函數(shù)。利用等價(jià)權(quán)原理，由上式可得：

(10)

單位權(quán)中誤差：

(11)

(12)

式中, |vi/σ|為第i個(gè)觀測(cè)值經(jīng)k次迭代后的標(biāo)準(zhǔn)化殘差；根據(jù)經(jīng)驗(yàn)c0一般為1～1.5,c1為3.0左右。通過(guò)抗差最小二乘曲面擬合獲得偏差擬合曲面，將輔軌道整體矯正，使其與主軌道參考基準(zhǔn)統(tǒng)一，最終實(shí)現(xiàn)升降軌PS-InSAR地面沉降監(jiān)測(cè)結(jié)果數(shù)據(jù)融合。

3 升降軌解算結(jié)果互檢與融合

3.1 升降軌垂直形變量解算結(jié)果

將38景降軌數(shù)據(jù)和46景升軌數(shù)據(jù)，利用PS-InSAR技術(shù)時(shí)序分析，獲取了研究區(qū)內(nèi)2018 年1月至2019年7月間地面沉降序列，同時(shí)利用式(4)將獲取雷達(dá)視線向的形變量轉(zhuǎn)換為垂直方向上的累計(jì)形變量，如圖2所示。可以看出，研究區(qū)內(nèi)分布著多個(gè)沉降中心，最大沉降量為-459 mm，除沉降中心外，其他區(qū)域的沉降量均小于50 mm。

3.2 升降軌解算結(jié)果內(nèi)部互檢

通過(guò)獲得升降軌公共覆蓋區(qū)域，提取該區(qū)域內(nèi)升降軌模式下各自的PS點(diǎn)，共提取降軌模式下25 809個(gè)PS點(diǎn)，升軌模式下31 049個(gè)PS點(diǎn)。由于雷達(dá)入射角以及采樣方式的差異，導(dǎo)致升降軌解算結(jié)果中同名相干點(diǎn)數(shù)極為有限，為實(shí)現(xiàn)升降軌解算結(jié)果互檢，本研究通過(guò)分別獲取抬升區(qū)和沉降區(qū)各30對(duì)同名相干點(diǎn)用于對(duì)比分析實(shí)驗(yàn)。如圖3所示，沉降區(qū)內(nèi)升降軌垂直形變量相關(guān)系數(shù)為0.866 5，線性擬合方程為y=0.941 4x-45.264 3； 抬升區(qū)內(nèi)升降軌垂直形變量相關(guān)系數(shù)為0.802 6，線性擬合方程為y=1.056 7x-33.536 8，表明升降軌模式下同名點(diǎn)垂直形變量具有較高的相關(guān)性。

(a)降軌垂直向形變量

(b)升軌垂直向形變量圖2 升降軌垂直向形變量圖Fig.2 Vertical Deformation of Lift Rail

3.3 升降軌解算結(jié)果融合方法對(duì)比

通過(guò)對(duì)升軌解算結(jié)果進(jìn)行插值處理，得到與降軌解算結(jié)果同名相干點(diǎn)的偏差數(shù)組，利用平均值法和抗差最小二乘擬合法分別對(duì)偏差數(shù)組進(jìn)行處理，并對(duì)升軌解算結(jié)果進(jìn)行校正，使升降軌解算結(jié)果參考基準(zhǔn)統(tǒng)一，實(shí)現(xiàn)升降軌解算結(jié)果的數(shù)據(jù)融合。從圖4中可明顯看出，圖4(a)中1、2區(qū)域存在多個(gè)與周圍不同形變量級(jí)的PS點(diǎn)，而圖4(b)中1、2區(qū)域幾乎所有PS點(diǎn)形變量級(jí)保持一致，因此，抗差最小二乘擬合法升降軌融合效果明顯優(yōu)于平均值法。選取研究區(qū)域內(nèi)升降軌均勻分布的30對(duì)同名相干點(diǎn)，將降軌相干點(diǎn)垂直形變量設(shè)為真值，分別驗(yàn)證平均值法和抗差最小二乘擬合法對(duì)升軌相干點(diǎn)的矯正程度，如圖5所示。通過(guò)三組數(shù)據(jù)對(duì)比計(jì)算得到，平均值法條件下升降軌同名相干點(diǎn)均方根誤差為6.8 mm，而抗差最小二乘擬合法條件下升降軌同名相干點(diǎn)均方根誤差為3.4 mm，表明在校正升降軌參考基準(zhǔn)中抗差最小二乘擬合法較平均值法精度提高了3.4 mm。

(a)沉降區(qū)線性擬合方程圖(b)抬升區(qū)線性擬合方程圖(c)沉降區(qū)垂直形變量相關(guān)性圖(d)抬升區(qū)垂直形變量相關(guān)性圖圖3 升降軌垂直形變量相關(guān)性圖Fig.3 Correlation diagram of vertical deformation of lift rail

(a)平均值法升降軌融合圖(b)抗差最小二乘擬合法升降軌融合圖圖4 升降軌解算結(jié)果數(shù)據(jù)融合方法對(duì)比圖Fig.4 Comparison chart of data fusion methods for lifting track solution results

圖5 平均值法和抗差最小二乘擬合法對(duì)比圖Fig.5 Comparison chart of mean method and robust least squares fitting method

4 地面沉降監(jiān)測(cè)結(jié)果與成因分析

利用抗差最小二乘擬合法融合升降軌解算，對(duì)東營(yíng)市進(jìn)行地表沉降監(jiān)測(cè)，圖6為2018年1月至2019年7月東營(yíng)市地表垂直形變量圖。東營(yíng)市地面沉降監(jiān)測(cè)結(jié)果顯示，東營(yíng)市主要沉降區(qū)位于東部沿海區(qū)域，城鎮(zhèn)沉降區(qū)主要位于廣饒城區(qū)，主要存在以下4個(gè)較為明顯的沉降區(qū)：

圖6 融合升降軌數(shù)據(jù)東營(yíng)市地表垂直形變量圖Fig.6 Surface vertical deformation map of Dongying city from ascending and descending data fusion

1)A沉降區(qū)位于河口區(qū)仙河鎮(zhèn)東部及西北部，其最大沉降量為-458 mm；B沉降區(qū)位于東營(yíng)市勝利機(jī)場(chǎng)東部及東北部，其最大沉降量為-202.3 mm；C沉降區(qū)位于廣南水庫(kù)南部廣饒鹽場(chǎng)附近，其最大沉降量為-290.6 mm。A、B、C三個(gè)沉降區(qū)皆存在著大范圍的鹽場(chǎng)、鹽田，沿海鹽場(chǎng)的主要生產(chǎn)模式是通過(guò)抽取大量的地下鹵水來(lái)曬鹽。因此，地下鹵水的大量開采是造成三個(gè)區(qū)域沉降的主要因素。

2)D沉降區(qū)位于廣饒縣城區(qū)北部及石村鎮(zhèn)，其最大沉降量為-187.5 mm，沉降面積約為40 km2，為東營(yíng)市面積最大的沉降區(qū)域，其呈現(xiàn)整體均勻沉降的趨勢(shì)。該區(qū)域存在多個(gè)工業(yè)園區(qū)，對(duì)水資源需求量較大，導(dǎo)致該區(qū)域深層地下水大量開，。廣饒縣深層地下水過(guò)度開采情況一直存在，其年平均開采量在3000萬(wàn)m3以上 。因此，過(guò)度開采深層地下水是造成D區(qū)域沉降的主要因素。

根據(jù)A、B、C、D四個(gè)沉降區(qū)的分布特點(diǎn)，可以得知東營(yíng)市地面沉降呈現(xiàn)東部沉降大于西部沉降，沿海地區(qū)沉降大于城區(qū)內(nèi)陸沉降。除過(guò)度開采鹵水和深層地下水之外，影響東營(yíng)市地面沉降的因素還有很多，東營(yíng)市在地質(zhì)構(gòu)造上位于郯廬斷裂帶的西側(cè)，斷裂帶活動(dòng)較為頻繁，在斷裂帶的影響下，地殼發(fā)生持續(xù)緩慢的凹陷，導(dǎo)致東營(yíng)市地面沉降現(xiàn)象不斷發(fā)展。從地質(zhì)土層方面來(lái)看，由于東營(yíng)市地處黃河三角洲，受黃河沖積影響，沉積層多為黏性土層，其自然固結(jié)時(shí)間短，土質(zhì)松軟易壓縮，容易發(fā)生地面沉降現(xiàn)象，隨著東營(yíng)市經(jīng)濟(jì)社會(huì)的不斷發(fā)展，地面建筑的不斷增多，加速了土層的壓實(shí)作用，加劇了地面沉降的程度。

5 結(jié)論

采用PS-InSAR技術(shù)獲取東營(yíng)市地表形變量，并利用抗差最小二乘擬合法進(jìn)行升降軌解算結(jié)果融合，實(shí)現(xiàn)PS-InSAR解算結(jié)果的內(nèi)部互檢及形變結(jié)果分析。

1)分別利用升降軌數(shù)據(jù)PS-InSAR技術(shù)進(jìn)行地面沉降監(jiān)測(cè)，實(shí)現(xiàn)解算結(jié)果的內(nèi)部互檢，升降軌解算結(jié)果中同名相干點(diǎn)垂直形變量相關(guān)系數(shù)在0.8以上，保證了地面沉降監(jiān)測(cè)的可靠性。

2)利用平均值法和抗差最小二乘擬合法修正升軌解算結(jié)果參考基準(zhǔn)，通過(guò)與降軌同名相干點(diǎn)位比較，抗差最小二乘擬合法均方根誤差3.4 mm，優(yōu)于平均值法的6.8 mm，提高了融合升降軌PS-InSAR解算結(jié)果精度。

3)通過(guò)融合升降軌PS-InSAR獲取東營(yíng)市地面沉降信息，加密單一軌道的觀測(cè)序列，豐富了地面相干點(diǎn)的空間分布，監(jiān)測(cè)結(jié)果顯示：2018年1月至2019年7月期間東營(yíng)市存在4個(gè)明顯沉降區(qū)，最大沉降量為-458 mm，最大沉降面積為40 km2。東營(yíng)市沉降整體呈現(xiàn)空間分布不均勻，且東部沿海大于西部城區(qū)，其沉降成因主要與過(guò)度開采地下鹵水與深層地下水有關(guān)。