鄧曉景，曲國慶，張建霞，席換，王暉

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

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

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

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

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

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

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

哨兵1號(Sentinel-1)衛(wèi)星是歐空局哥白尼計劃(GMES)中的地球觀測衛(wèi)星，發(fā)射于2014年，由兩顆衛(wèi)星組成，載有C波段合成孔徑雷達，本研究選擇覆蓋東營市的Sentinel_1A雷達衛(wèi)星數(shù)據(jù)，該數(shù)據(jù)采樣模式為IW寬幅干涉模式，在該模式下采用TOPS (terrain observation with progressive scans)掃描方式，可以同時獲取3個子條帶，其幅寬為250 km,重訪周期為12 d，空間分辨率為5 m×20 m[9]。實驗數(shù)據(jù)選取分別為降軌數(shù)據(jù)(Track-76)左側(cè)子條帶和升軌數(shù)據(jù)(Track-69)中間子條帶，數(shù)據(jù)參數(shù)見表1。參考DEM數(shù)據(jù)為美國宇航局提供的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ù)處理、差分干涉計算和形變量計算三個步驟。首先，選取覆蓋研究區(qū)的N+1幅SAR影像，根據(jù)基線參數(shù)和多普勒效應(yīng)，選取一幅SAR影像作為公共主影像，其余影像為輔影像，分別與公共主影像配準進行干涉處理，生成N幅干涉圖；通過精密軌道文件去除平地效應(yīng)，同時利用參考DEM數(shù)據(jù)去除地形相位，對N幅干涉圖進行差分處理，生成差分干涉圖[10]。其中，任意差分干涉圖像元的相位Δφ可以表示為：

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

(1)

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

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

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

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

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

(2)

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

dLOS≈dEsinθ+dVcosθ

(3)

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

dV≈dLOS/cosθ

(4)

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

(5)

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

2.3 抗差最小二乘曲面擬合

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

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

(6)

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

其誤差方程矩陣形式為:

V=KM-Z

(7)

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

K=(MTPM)-1MTPZ

(8)

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

(9)

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

(10)

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

(11)

(12)

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

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

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

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

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

通過獲得升降軌公共覆蓋區(qū)域，提取該區(qū)域內(nèi)升降軌模式下各自的PS點，共提取降軌模式下25 809個PS點，升軌模式下31 049個PS點。由于雷達入射角以及采樣方式的差異，導(dǎo)致升降軌解算結(jié)果中同名相干點數(shù)極為有限，為實現(xiàn)升降軌解算結(jié)果互檢，本研究通過分別獲取抬升區(qū)和沉降區(qū)各30對同名相干點用于對比分析實驗。如圖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，表明升降軌模式下同名點垂直形變量具有較高的相關(guān)性。

(a)降軌垂直向形變量

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

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

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

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

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

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

圖6 融合升降軌數(shù)據(jù)東營市地表垂直形變量圖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ū)位于東營市勝利機場東部及東北部，其最大沉降量為-202.3 mm；C沉降區(qū)位于廣南水庫南部廣饒鹽場附近，其最大沉降量為-290.6 mm。A、B、C三個沉降區(qū)皆存在著大范圍的鹽場、鹽田，沿海鹽場的主要生產(chǎn)模式是通過抽取大量的地下鹵水來曬鹽。因此，地下鹵水的大量開采是造成三個區(qū)域沉降的主要因素。

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

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

5 結(jié)論

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

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

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

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