應(yīng)用時空復(fù)合趨勢面Kriging 方法插值地表沉降

黃丙湖，李欣芮，范芷睿，潘海燕，廖一蘭

（1. 中國石油大學(xué)（華東） 海洋與空間信息學(xué)院，山東 青島 266580；2. 新疆大學(xué)建筑工程學(xué)院，新疆 烏魯木齊 830046；3. 中國科學(xué)院地理科學(xué)與資源研究所資源與環(huán)境國家重點(diǎn)實(shí)驗(yàn)室，北京 100101）

地鐵施工中土質(zhì)、荷載等因素將導(dǎo)致地面沉降，易引發(fā)坍塌、開裂等事故，因此必須建立相應(yīng)的觀測體系進(jìn)行風(fēng)險(xiǎn)評估，以便采取相應(yīng)措施防止事故災(zāi)害的發(fā)生。目前的常用方法為在地鐵施工線上設(shè)立地表沉降監(jiān)測點(diǎn)，并定期觀測監(jiān)測點(diǎn)沉降數(shù)據(jù)，但所得的采樣數(shù)據(jù)較離散，無法直觀體現(xiàn)整體沉降趨勢，必須對采樣數(shù)據(jù)進(jìn)行插值[1-2]。空間插值考慮了沉降監(jiān)測點(diǎn)在空間上的相關(guān)性，能對整體沉降區(qū)域的監(jiān)測點(diǎn)進(jìn)行預(yù)測。常用的空間插值方法包括幾何方法、統(tǒng)計(jì)方法、函數(shù)方法和綜合方法[3]，其中以Kriging 為代表的地學(xué)統(tǒng)計(jì)方法可定量描述空間隨機(jī)分布變量的相關(guān)性，已應(yīng)用于測繪、地質(zhì)學(xué)、GIS 等學(xué)科，具有高度可行性和有效性[4-8]。陳雷[9]等利用Kriging 插值方法對地面沉降數(shù)據(jù)進(jìn)行內(nèi)插，但選取樣本量小，只能反映局部趨勢；劉志平[10]等利用Kriging插值方法研究了邊坡變形速率，但其將變形速率演化為與時間相關(guān)的函數(shù)，弱化了時空關(guān)聯(lián)性。在對地表沉降量等時空變量進(jìn)行插值前，除分析變量空間相關(guān)性外，還應(yīng)考慮時間相關(guān)性，王建民[11]等利用時空Kriging插值方法，將時間域擴(kuò)展到空間域，建立了時空半變異函數(shù)，該方法計(jì)算量大，選擇變異函數(shù)模型時需人工干預(yù)，會造成一定誤差；Kyriakidis P C[12]等利用時空Kriging 插值方法研究區(qū)域變化量的時空分布特征，但未與傳統(tǒng)地統(tǒng)計(jì)插值方法進(jìn)行比較，無法體現(xiàn)其優(yōu)勢。時空復(fù)合趨勢面（CSTM）兼顧了沉降點(diǎn)的時間相關(guān)性和空間相關(guān)性，利用已知沉降點(diǎn)監(jiān)測值進(jìn)行趨勢面擬合得到區(qū)域內(nèi)每個時空變量的全局平均趨勢。因此，本文在傳統(tǒng)地統(tǒng)計(jì)插值模型的基礎(chǔ)上引入了CSTM，使模型能更準(zhǔn)確地反映時空變量特征，從而提高模型預(yù)測精度；并與傳統(tǒng)方法進(jìn)行了精度對比，驗(yàn)證了時空插值的優(yōu)勢。

1 基于CSTM的Kriging插值方法

1.1 CSTM

CSTM 首先通過區(qū)域內(nèi)已知時空變量的屬性值進(jìn)行趨勢面擬合，再利用指數(shù)平滑法得到每個時空變量的屬性值。該方法將區(qū)域內(nèi)每個屬性值分為趨勢值和殘差值兩個部分，以此反映時空變量的變異情況。Lee[13]等提出CSTM 是基于時空隨機(jī)場等理論定義的。時空隨機(jī)場[14]即z(p)=z(s't)，包括時間尺度和空間尺度兩個維度，s為二維空間位置（本文中為沉降監(jiān)測點(diǎn)空間位置），t為時間（本文中為沉降監(jiān)測點(diǎn)觀測時間），將Y(p)=log[z(p)]作為z(p)的對數(shù)變換。全局平均趨勢mY(p)是由Y(p) 構(gòu)建的確定性函數(shù)，反映Y(p)的全局平均趨勢。CSTM是基于復(fù)合時間/空間指標(biāo)獲得mY(p)的有效方法[13]。傳統(tǒng)理論獲得mY(p)（即mY(s't)）的方法是將其分為時間和空間兩部分，而CSTM的優(yōu)勢在于兼顧了時間分量與空間分量[15]。

對于一個周期性變量，CSTM 可利用基于指數(shù)平滑函數(shù)的權(quán)重來計(jì)算擬合函數(shù)，以獲得每個時空坐標(biāo)平滑的平均趨勢值，在任意時空坐標(biāo)pj=(sj'tj)下，平均趨勢值的計(jì)算公式為[13]：

式中，Y(sj'tj)為沉降量在可測量監(jiān)測點(diǎn)pj處的對數(shù)變換；wi為權(quán)重；d(si'sj)為沉降監(jiān)測點(diǎn)pi、pj之間的歐氏距離，；as、at分別為指數(shù)平滑函數(shù)的空間范圍和時間范圍。

1.2 基于CSTM的Kriging插值

CSTM 將區(qū)域內(nèi)每個屬性值分成趨勢值和殘差值，基于CSTM的Kriging插值是對殘差進(jìn)行插值，再與趨勢值求和，最后得到時空變量的預(yù)測值。

1）數(shù)據(jù)空間分布的判斷。使用Kriging 插值方法的前提是待插值數(shù)據(jù)滿足正態(tài)分布且具有空間自相關(guān)性，因此可利用Moran’s I等空間自相關(guān)性統(tǒng)計(jì)量指標(biāo)判定待插值數(shù)據(jù)的空間分布特征，滿足該前提才能使用該方法。

2）時空尺度的確定。式（2）中as的選取原則為小于等于兩個測量沉降監(jiān)測點(diǎn)之間的最大歐氏距離，at的設(shè)置原則為小于等于沉降監(jiān)測點(diǎn)多次觀測的最大時間差。根據(jù)式（2）計(jì)算得到指數(shù)平滑函數(shù)的權(quán)重。

3）mY(p)的計(jì)算。沉降過程中根據(jù)已知監(jiān)測點(diǎn)沉降量來構(gòu)建CSTM，利用式（1）計(jì)算得到每個時空監(jiān)測點(diǎn)的mY(p)。

4）Kriging插值。Y(p)由確定性變換X(p)=Y(p)-mY(p) 產(chǎn)生，由局部穩(wěn)定平均值mx(p)=E[X(p)]=E[Y(p)-mY(p)]獲得固有平穩(wěn)時空隨機(jī)場，其中E[.]表示期望算子[13]，X(p) 為殘差。由于沉降監(jiān)測點(diǎn)p=(s't)和p'=(s''t')之間空間滯后γ=‖s-s'‖ 和時間滯后τ=|t-t' |的影響產(chǎn)生誤差，二者之間的總體誤差用協(xié)方差計(jì)算，能衡量一對沉降監(jiān)測點(diǎn)間的時空相關(guān)性，即

Kriging 插值方法考慮了樣本之間的空間相關(guān)性，通過最小化均方誤差權(quán)重來獲得未監(jiān)測沉降點(diǎn)殘差的最佳估計(jì)值。X(p)在點(diǎn)k處的估計(jì)值⌒是測量值xd的線性組合，計(jì)算公式為：

式中，λ為采用傳統(tǒng)Kriging 插值估算權(quán)重的列向量，沉降點(diǎn)pk、pd的時空復(fù)合距離越近，權(quán)重越大；mx(pk)、mx(pd)為X(p)在預(yù)測點(diǎn)pk、pd的平均趨勢。

式中，cx(pk'pd)為估計(jì)點(diǎn)與數(shù)據(jù)點(diǎn)之間X(p)的協(xié)方差行向量；cx(pd'pd)為數(shù)據(jù)點(diǎn)之間X(p)的協(xié)方差矩陣[16]。

5）監(jiān)測量插值結(jié)果的估計(jì)。將區(qū)域內(nèi)每個沉降監(jiān)測點(diǎn)的mY pj=(sj'tj) 與Kriging 插值得到的殘差對應(yīng)求和，最后進(jìn)行指數(shù)變換得到每個監(jiān)測點(diǎn)的沉降量。

2 實(shí)驗(yàn)與結(jié)果分析

2.1 數(shù)據(jù)來源

本文的實(shí)驗(yàn)數(shù)據(jù)來源于烏魯木齊市軌道交通3 號線一期工程05標(biāo)明園站—友好站區(qū)間，該工程處于施工階段，以2 號豎井和橫通道開挖。考慮到地層自穩(wěn)能力差，受開挖、支護(hù)工程和地下水等因素影響，地表會出現(xiàn)不同程度沉降，因此豎井施工過程中針對可能沉降變化較大的區(qū)域布設(shè)41個監(jiān)測點(diǎn)，構(gòu)成了6條監(jiān)測線（圖1），在井口頂部圈梁短邊周圍布設(shè)6個監(jiān)測點(diǎn)，其余監(jiān)測點(diǎn)分別布設(shè)在豎井橫通道附近和隧道線路附近。采用徠卡電子水準(zhǔn)儀（±0.3 mm/km），按照二等水準(zhǔn)測量方法，對地表監(jiān)測點(diǎn)進(jìn)行觀測，經(jīng)平差合格后，獲得監(jiān)測點(diǎn)處的高程；再統(tǒng)計(jì)觀測期間監(jiān)測點(diǎn)的高程變化情況，獲得沉降變化量。實(shí)驗(yàn)獲取了2017-06-24—2017-09-08 共15 期沉降數(shù)據(jù)，監(jiān)測頻率為5 d/次，監(jiān)測過程中個別監(jiān)測點(diǎn)數(shù)據(jù)缺失，導(dǎo)致總數(shù)據(jù)量小于n=41×15。

圖1 地鐵地表沉降監(jiān)測點(diǎn)布設(shè)圖

2.2 插值過程

首先，對獲取的數(shù)據(jù)進(jìn)行空間自相關(guān)性檢驗(yàn)，利用Moran’s I 進(jìn)行判定，p值小于0.05 且z 得分超過1.65，認(rèn)為通過95%置信度檢驗(yàn)，結(jié)果見表1，可以看出，4個時段的數(shù)據(jù)未通過顯著性檢驗(yàn)，且它們的Moran’s I 均接近于0，空間自相關(guān)性很弱，將其剔除。本文分別利用Kriging插值和基于CSTM的Kriging插值方法對研究區(qū)內(nèi)符合條件的11個時段的監(jiān)測點(diǎn)沉降量進(jìn)行插值。實(shí)驗(yàn)中選取每期監(jiān)測數(shù)據(jù)利用Kriging方法進(jìn)行獨(dú)立插值發(fā)現(xiàn)，預(yù)測沉降監(jiān)測點(diǎn)與半徑為2 m 區(qū)域內(nèi)的已知沉降點(diǎn)空間相關(guān)性最強(qiáng)，因此確定CSTM基于指數(shù)平滑函數(shù)權(quán)重的空間尺度為2 m，時間尺度選取待插值數(shù)據(jù)的前后1 個時間段。值得注意的是，由于剔除了7 月30 日—8 月4 日的數(shù)據(jù)，導(dǎo)致該時段周期變?yōu)?0 d，即7月30日—8月9日，本文將這種時段相距天數(shù)不相等的間隔稱為非等時間間隔，反之稱為等時間間隔。然后，計(jì)算區(qū)域內(nèi)每個沉降監(jiān)測點(diǎn)的mY(p)。最后，對每個地表沉降監(jiān)測點(diǎn)的殘差進(jìn)行Kriging插值。

表1 實(shí)驗(yàn)數(shù)據(jù)空間自相關(guān)性檢驗(yàn)

兩種方法的插值結(jié)果見圖2，可以看出，圈梁、豎井橫通道附近的沉降程度相對明顯，因此施工過程中需加強(qiáng)監(jiān)測這兩個部分；基于CSTM的Kriging插值方法考慮了時空相關(guān)性，7月15—25日圈梁附近地表緩慢下降，而7月25日—8月14日則開始抬升，后期逐步趨于穩(wěn)定狀態(tài)；豎井和橫通道附近的點(diǎn)在不同時段內(nèi)沉降波動較大，這可能是由施工推進(jìn)造成的。

2.3 精度評價(jià)

實(shí)驗(yàn)采用留一交叉驗(yàn)證的方法，即對于N個樣本，每次選擇N-1 個樣本訓(xùn)練數(shù)據(jù)，一個樣本驗(yàn)證模型預(yù)測的好壞。該方法的優(yōu)勢在于，每個回合中幾乎所有的沉降數(shù)據(jù)皆用于訓(xùn)練模型，因此最接近原始沉降數(shù)據(jù)的分布，評估結(jié)果比較可靠。首先假定某一監(jiān)測點(diǎn)沉降量為未知，然后利用插值方法計(jì)算該監(jiān)測點(diǎn)沉降量，最后計(jì)算該監(jiān)測點(diǎn)實(shí)際觀測值與插值得到的沉降值之間的誤差。這種誤差主要是由插值算法的準(zhǔn)確性和監(jiān)測點(diǎn)的分布情況造成的[17]。實(shí)驗(yàn)選取均方根誤差（RMSE）和R2作為精度驗(yàn)證指標(biāo)（表2），RMSE 可反映監(jiān)測點(diǎn)數(shù)據(jù)插值結(jié)果的離散程度，計(jì)算公式為：

表2 兩種方法的RMSE和R 2 比較

式中，vi為第i個監(jiān)測點(diǎn)的實(shí)際沉降量；v'i為第i個監(jiān)測點(diǎn)的差值估算值；n為參與交叉驗(yàn)證的監(jiān)測點(diǎn)數(shù)目。

R2可體現(xiàn)模型的擬合效果，由殘差平方和總體平方組成，其值越接近1，表明線性相關(guān)程度越高，計(jì)算公式為：

由表2 可知，空間插值和時空插值都表現(xiàn)出較強(qiáng)的線性相關(guān)性，本文方法比傳統(tǒng)方法的RMSE 降低約35%，說明本文方法的總體離散程度明顯小于傳統(tǒng)方法；傳統(tǒng)方法的R2大部分在0.5以下，而本文方法的擬合程度更好；綜合來看，本文方法提高了模型精度，更適用于地表沉降監(jiān)測。

然而，CSTM 未給出標(biāo)準(zhǔn)的時間尺度確定方法，因此本文在前后一個以非等間隔時間段為指數(shù)平滑的時間尺度實(shí)驗(yàn)之后，又將時間尺度設(shè)為等間隔時間段進(jìn)行探索性實(shí)驗(yàn)，結(jié)果見表3，可以看出，采用等間隔時間段的沉降數(shù)據(jù)構(gòu)建CSTM 得到的精度更好。插值結(jié)果見圖3。

表3 不同時間尺度的本文方法的精度評價(jià)

圖3 同時期非等間隔與等間隔本文方法插值結(jié)果比較

CSTM 同時考慮了沉降監(jiān)測點(diǎn)的時間相關(guān)性與空間相關(guān)性，因此選取研究區(qū)域內(nèi)均勻時段進(jìn)行預(yù)測將得到更好的效果，由于本次實(shí)驗(yàn)數(shù)據(jù)量的局限性，等間隔沉降量差值的時間序列較短，未來可選取更長時間序列進(jìn)行實(shí)驗(yàn)。

3 結(jié) 語

本文將CSTM 引入Kriging 插值，同時考慮了空間相關(guān)性和時間相關(guān)性，使該插值方法能更真實(shí)地反映時空變量變化趨勢，提高了地表沉降監(jiān)測精度。實(shí)驗(yàn)表明，與傳統(tǒng)方法相比，本文方法的RMSE 降低約35%，具有更高的擬合精度；從插值結(jié)果來看，能更精細(xì)地反映沉降量的時空差異，以實(shí)現(xiàn)對地鐵地表沉降的有效監(jiān)測。受限于數(shù)據(jù)觀測天數(shù)較少，本文具有一定局限性，未來可考慮對更長時間的沉降監(jiān)測進(jìn)行實(shí)驗(yàn)，同時可引入周邊建筑物、土質(zhì)、施工方式、人口密度、車流量等輔助變量來優(yōu)化插值模型，實(shí)現(xiàn)長時間、大范圍的地表沉降趨勢分析。除此之外，本文方法還可與三維激光掃描數(shù)據(jù)的點(diǎn)云數(shù)據(jù)相結(jié)合，對格網(wǎng)化后的點(diǎn)云數(shù)據(jù)進(jìn)行插值，實(shí)現(xiàn)更高效快捷的整體監(jiān)測分析。