利用水準數據的天津市地面沉降時空特征分析

常 明，尹海權，蘇廣利，田 曉，徐玉健

(中國地震局第一監測中心，天津 300180)

地面沉降是區域性地面高程降低的環境地質災害，具有沉降速度相對較慢、難以察覺、持續時間長、影響面積廣、防治難度大、幾乎不可恢復等特點。地面沉降會導致樓房墻壁開裂、地基下沉、鐵路路基下沉的問題，給人們的生產、生活造成巨大損失[1-3]。針對地面沉降的監測和分析已經越來越引起人們的重視。

造成地面沉降的因素分為自然因素和人為因素。其中，自然因素主要包括地殼運動、地震、地質災害等，具有地質構造發生變化使得地面出現相應的形變特征[4-7]。人為因素則分為外部因素和內部因素。其中，內部因素主要包括地下水、石油、天然氣及礦物質等地層物質的長期超量開采；外部因素則主要是大規模的施工建設導致的地面高程變化。目前我國遭受地面沉降災害的城市超過50個，地下水的過量開采是導致我國城市地面沉降的主要原因[8-13]。

國內許多學者借助多種觀測手段對地面沉降情況進行觀測，并借助趨勢面分析等算法模型對地面沉降情況進行了分析處理。文獻[14]利用InSAR手段對北京地面沉降的時空分布特征進行了分析。文獻[15]利用GPS對天津市沉降情況進行處理并與水準觀測結果進行比對，證明了GPS在平原區域對沉降情況進行監測的可行性。文獻[16]利用目標層次分析法構建了滄州地區地面沉降風險評估指標體系，并對指標體系進行定性分析，借助GIS將滄州分成了高中低3級地面沉降危險區。文獻[17]以溫黃平原為例，借助多種模型算法對離散沉降觀測數據進行擬合，構建地面沉降量與時間序列的關系，為地面沉降預測模型選擇和沉降監測防控提供了依據。文獻[18]綜合考慮建筑物荷載對地面沉降的影響，對Kriging差值方法進行改進，準確反映了沉降形變量大小和趨勢。文獻[19]在趨勢面擬合的基礎上添加時間序列，構建三維動態趨勢面模型，獲得了良好的時空插值結果，并使用GPS沉降監測網數據對地面沉降情況進行了分析。

雖然目前有多種觀測手段用于地面沉降監測，但水準觀測是最準確、最能直觀反映地面高程變化的觀測手段之一。通過測量地面水準點的高程變化，可以反映地面的沉降變化。自1923年天津市經歷了地面沉降的萌芽期、沉降中心的形成、地面沉降快速擴展幾大階段，直至20世紀70年代才開始針對地面沉降進行研究、治理。目前天津市已經成為我國北方地區地面沉降最嚴重的城市之一，持續的地面沉降導致經濟的發展受到損失、人們的生活受到威脅[20-22]。本文采用2005—2012年、2016—2017年水準觀測數據，利用擴展卡爾曼濾波模型對觀測數據進行濾波，求解沉降量、沉降速率及沉降加速率，構建多項式加權內插模型，對天津市地面沉降的時空特征進行分析。

1 水準監測情況概述

天津市地處華北平原東北部，渤海西岸，除薊縣外均為低山丘陵。為了避免降水等其他自然因素造成地下水位變化從而導致地面沉降的情況，天津市的地面沉降動態監測一、二等水準測量工作都會選擇在每年8—10月進行，其中一等水準路線長度約1300 km。考慮地形的影響，一等水準點的布設主要集中在東麗、大港、塘沽等除薊縣以外區域。具體分布如圖1所示。

圖1 天津市水準點分布情況

隨著城市的快速發展，礙于施工影響，部分水準點遭到破壞，需要重新埋設。為了保證數據的連續性，本文選擇收集水準觀測成果中未被破壞的水準點進行地面沉降分析。

2 卡爾曼濾波模型

地面沉降過程可以看作是一個動態系統，因此水準數據是動態測量平差的過程。選擇一個合適的動態數學模型對其進行數據處理能夠更加準確地獲取地面的沉降信息。水準測量的觀測方程模型為

L=BX+Δ

(1)

式中，L為觀測值；B為已知系數矩陣；Δ為觀測誤差向量，即觀測噪聲向量。因此求矩陣X的過程即為測量平差。在測量平差中，含有誤差的觀測值求解參數的最佳估值方法分為兩種：最小二乘平差法和濾波。由于濾波考慮了參數的隨機性質，因此當已知參數的先驗特性時，所得到的估值比經典最小二乘的方法具有更高的精度。地面沉降可以理解為處于運動中的地面控制點與時間t的函數。為了能夠更好地對天津市沉降信息進行描述，本文以沉降量、沉降速率、沉降加速率作為狀態向量，構建狀態方程式如下

Xk+1=Φk+1,kXk+Γk+1ωk

(2)

式中，Xk為k時刻的狀態向量；Xk+1為k+1時刻的狀態向量；ωk為隨機干擾；Φk+1,k為狀態轉換參數,則

(3)

式中，h為某一水準點位的高程；v為沉降速率；a為沉降加速率。

Φk+1,k為狀態轉移矩陣

(4)

卡爾曼濾波的觀測方程式可表述為

Lk+1=Bk+1Xk+1+Δk+1

(5)

式中，Lk+1為k+1時刻的觀測值；Bk+1為k+1時刻的已知系數矩陣；Δk+1為k+1時刻的觀測誤差。

隨機模型為

(6)

在水準觀測測量中，假設某一個控制點i在t時刻的沉降信息為ξi(t),其沉降瞬時速率為λi(t)，將地面沉降的瞬時加速率Ωi(t)，則

(7)

(8)

針對一個控制點的狀態向量的描述可認為是沉降信息、瞬時速率及瞬時加速率的函數，記為

f(ξi(t),λi(t),Ωi(t))

(9)

令i點的狀態向量為Xi(t)，所有水準的狀態向量記為X(t)，則天津市的沉降水準觀測網的狀態模型可以表述為

(10)

水準測量中，以觀測點的高程和沉降速率為狀態向量，則狀態方程和觀測方程為

(11)

hi(k+1)=Bk+1Xk+1+Δi(k+1)

(12)

通過以上兩個方程可實現對天津市水準數據的濾波，得到水準觀測點的沉降量、沉降速率、沉降加速率。

3 地面沉降空間分布

水準觀測獲取固定點相對高程的變化，為了更加直觀地表現沉降的空間分布特征，本文針對求解的參數對沉降信息進行內插。常見的內插數據模型有很多種，主要包括反距離加權插值、線性內插等函數模型。

3.1 反距離加權插值

反距離加權插值是對于當前插值點i，選擇一定大小l的窗口。對窗口內的觀測點集Pn的沉降量按照距離倒數k次方的方法進行加權計算，k值可以根據實際大小進行選擇。具體公式為

(13)

式中，hi為內插點沉降量；hj為點集Pn中第j個點的沉降量；dij為觀測點與目標點之間的距離；σ為平滑因子，一般取為0.002；k為權的方次，一般可取為1。

3.2 線性內插

線性內插是假設試驗區域內的沉降量為線性變化，內插點i的取值認為是關于經度、緯度的函數，公式為

hi=a0+a1Bi+a2Li

(14)

式中，Bi、Li為內插點的經緯度。

3.3 多項式加權內插

地面沉降并不完全呈線性變化，且內插點沉降量的大小并不僅僅受距離的影響。距離、沉降速率、沉降加速率都是影響內插點沉降的因素。因此本文針對目標內插點i，選擇一定大小的窗口，并提取窗口內觀測點濾波后的狀態向量，構建多項式加權內插模型，對目標點的沉降量進行求解，公式如下

(15)

式中，hi為目標內插點i的沉降量；hj為窗口內點集Pn中觀測點j的沉降量；δij為觀測點j的權重；θij為由距離解算的權重；Δvj為觀測點j的沉降速率；Δaj為觀測點j的沉降加速率；dij為觀測點與目標點之間的距離；σ為平滑因子；k1、k2、k3分別為距離、沉降速率、沉降加速率的方次。

通過窗口的不斷移動，由以上公式實現對試驗區域內每個點的沉降量進行求解。

3.4 精度評定

本文利用觀測數據中80%作為試驗數據進行內插，利用20%作為驗證數據進行精度評定。通過比較驗證數據與內插數據的偏移量的絕對值，驗證多項式加權內插的準確性。統計結果見表1。

表1 試驗數據內插結果統計 mm

通過試驗數據的比較分析，多項式加權內插相比于其他方法具有更高的精度。

4 試驗分析

對天津市水準觀測數據進行濾波，并利用本文提出的內插方式進行處理分析發現，2005—2017年天津市整體區域累計形變量分布情況如圖2所示，最大形變量為1 143.5 mm；累計平均形變量為394.477 5 mm。以北辰區西部、漢沽區南部、靜海中部、東麗南部、大港中部、西青南部為中心，呈現漏斗式區域下沉結果，其中以北辰區地面沉降情況最為嚴重。

圖2 天津市2005—2017年總體沉降量

為了能夠更加直觀地查看沉降主要區域的時空變化情況，對天津市每年的沉降速率進行求解，如圖3所示。

由圖3(a)可知，2005—2017年天津市每年都在發生不均勻沉降，呈現漏斗式沉降，大體以北辰、西青、大港、塘沽等為主要區域，且北辰區沉降分布范圍逐漸擴大，沉降量逐漸增加。由圖3(e)—(g)可知，2009—2011年漢沽地區沉降量逐步增加，2011—2012年漢沽地區沉降量達到最大，速率達到220 mm/a；2012—2016年漢沽地區地面沉降呈現減弱趨勢；截至2017年，天津市地面沉降空間分布大體趨于穩定，僅北辰西部與大港中西部仍存在沉降情況。

圖3 平均形變速率

5 結 語

本文利用天津市2005—2012年、2016—2017年的水準觀測數據，以水準點的沉降量、沉降速率、沉降加速率為狀態向量，利用卡爾曼濾波算法對水準觀測數據進行了處理；并根據得到的結果構建了多項式加權內插模型，對內插結果進行了可視化展示。分析2005年以來天津市地面沉降情況的時空演化特征可對天津市的城市發展及建設提供一定的參考意義。