融合升降軌PS-InSAR東營市地面沉降監測

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

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

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

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

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

1 研究區概況及數據選取

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

圖1 研究區概況及數據范圍Fig.1 Overview of study area and data range

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

表1 Sentinel-1數據參數Tab.1 Sentinel-1 data parameters

2 研究方法

2.1 PS-InSAR技術

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

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

(1)

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

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

2.2 升降軌數據融合簡化模型

PS-InSAR技術獲取的地表形變量可分為東西方向和垂直方向的地表形變量，對于同一研究區的升降軌數據，PS-InSAR技術獲取的觀測結果東西方向上水平形變量呈現相反的特點，而垂直方向上垂直形變量基本相同[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)

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

(5)

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

2.3 抗差最小二乘曲面擬合

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

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

(6)

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

其誤差方程矩陣形式為:

V=KM-Z

(7)

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

K=(MTPM)-1MTPZ

(8)

當偏差數據含有粗差時,如果使用初始權值,會使擬合曲面差生彎曲，影響擬合結果精度[15]。利用抗差估計理論，構造極值函數為：

(9)

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

(10)

單位權中誤差：

(11)

(12)

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

3 升降軌解算結果互檢與融合

3.1 升降軌垂直形變量解算結果

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

3.2 升降軌解算結果內部互檢

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

(a)降軌垂直向形變量

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

3.3 升降軌解算結果融合方法對比

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

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

(a)平均值法升降軌融合圖(b)抗差最小二乘擬合法升降軌融合圖圖4 升降軌解算結果數據融合方法對比圖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 地面沉降監測結果與成因分析

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

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

1)A沉降區位于河口區仙河鎮東部及西北部，其最大沉降量為-458 mm；B沉降區位于東營市勝利機場東部及東北部，其最大沉降量為-202.3 mm；C沉降區位于廣南水庫南部廣饒鹽場附近，其最大沉降量為-290.6 mm。A、B、C三個沉降區皆存在著大范圍的鹽場、鹽田，沿海鹽場的主要生產模式是通過抽取大量的地下鹵水來曬鹽。因此，地下鹵水的大量開采是造成三個區域沉降的主要因素。

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

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

5 結論

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

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

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

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