石 強(qiáng)，戴吾蛟，晏慧能，劉 寧

1. 中南大學(xué)測(cè)繪與遙感科學(xué)系，湖南 長(zhǎng)沙 410083； 2. 江蘇海洋大學(xué)海洋技術(shù)與測(cè)繪學(xué)院，江蘇 連云港 222005； 3. 湖南省精密工程測(cè)量與形變?yōu)暮ΡO(jiān)測(cè)重點(diǎn)實(shí)驗(yàn)室，湖南 長(zhǎng)沙 410083

隨著傳感器網(wǎng)絡(luò)、GNSS、InSAR等現(xiàn)代變形監(jiān)測(cè)技術(shù)的發(fā)展，人們獲取了大量的變形數(shù)據(jù)[1-3]，這些變形數(shù)據(jù)具有較為復(fù)雜的時(shí)空相關(guān)性，傳統(tǒng)的空間和時(shí)間序列分析方法已無法滿足時(shí)空變形分析的需求，因此一些學(xué)者開始將時(shí)空數(shù)據(jù)處理的理論與方法(例如時(shí)空地統(tǒng)計(jì)學(xué)模型、時(shí)空自回歸模型、時(shí)空地理加權(quán)回歸模型及時(shí)空Kalman濾波等)引入變形分析領(lǐng)域[4]。時(shí)空地統(tǒng)計(jì)學(xué)模型是在傳統(tǒng)地統(tǒng)計(jì)學(xué)的基礎(chǔ)上，運(yùn)用時(shí)空插值及隨機(jī)模擬估算和量化時(shí)空變化，動(dòng)態(tài)反映地理空間的變化過程，其代表性的方法是時(shí)空克里金，在邊坡[5-6]及隧道圍巖[7]變形分析中得到了較好的應(yīng)用。時(shí)空克里金模型雖然同時(shí)考慮了變形數(shù)據(jù)的時(shí)間和空間相關(guān)性，但是它計(jì)算量較大，不能滿足變形監(jiān)測(cè)大數(shù)據(jù)時(shí)空分析的需求。時(shí)空自回歸模型是根據(jù)空間權(quán)矩陣及AR模型描述時(shí)空數(shù)據(jù)的時(shí)空相關(guān)性，已用于大壩[8-9]、橋梁[10]、地下管線[11]及建筑物[12]時(shí)空變形預(yù)測(cè)，且具有較高的預(yù)測(cè)精度，但是該模型分別考慮空間和時(shí)間上的相關(guān)性，并不是真正意義上的時(shí)空分析模型。時(shí)空地理加權(quán)回歸模型是在普通地理加權(quán)回歸模型的基礎(chǔ)上進(jìn)行時(shí)空擴(kuò)展，構(gòu)建了顧及時(shí)空相關(guān)性的時(shí)空回歸模型，它已被用于大壩變形分析[13-14]，但是該模型需要對(duì)監(jiān)測(cè)點(diǎn)的位移與影響因素建立多元變系數(shù)回歸關(guān)系，計(jì)算量也較大。時(shí)空Kalman濾波是一種融入空間模型的Kalman濾波方法，在時(shí)空數(shù)據(jù)分析方面具有很大的優(yōu)越性，能對(duì)時(shí)空數(shù)據(jù)進(jìn)行時(shí)空動(dòng)態(tài)濾波、時(shí)空插值及時(shí)空預(yù)測(cè)，已被廣泛應(yīng)用于眾多時(shí)空數(shù)據(jù)分析領(lǐng)域，如傳染病動(dòng)態(tài)演變[15]、遙感影像融合[16]、空氣污染物建模[17]、氣候變化[18]等領(lǐng)域。文獻(xiàn)[19]提出了第一個(gè)具有統(tǒng)計(jì)意義的時(shí)空Kalman濾波，發(fā)展形成了3種代表性的時(shí)空Kalman濾波模型：Kriged Kalman filter(KKF)模型[20]、space time Kalman filter(STKF)模型[21]、spatio-temporal mixed effects (STME)[22-23]模型。文獻(xiàn)[24—27]將時(shí)空Kalman濾波引入變形監(jiān)測(cè)領(lǐng)域，并在動(dòng)態(tài)濾波插補(bǔ)GPS地面沉降數(shù)據(jù)、融合GNSS和InSAR沉降數(shù)據(jù)獲取高時(shí)空分辨率地表沉降數(shù)據(jù)及大壩變形的時(shí)空預(yù)測(cè)方面取得了良好的應(yīng)用效果。但上述研究?jī)H對(duì)時(shí)空Kalman濾波用于變形分析進(jìn)行了初步探索，未對(duì)已有典型的時(shí)空Kalman濾波模型的理論及其在變形分析中的適用性進(jìn)行深入分析。

為了更好地利用時(shí)空Kalman濾波進(jìn)行時(shí)空變形分析，本文從基本原理及試驗(yàn)兩方面對(duì)比分析KKF、STKF及STME 3種典型時(shí)空Kalman濾波模型的異同點(diǎn)及其在變形分析應(yīng)用中的適用性，并對(duì)比分析其在時(shí)空濾波去噪、數(shù)據(jù)插補(bǔ)及變形預(yù)測(cè)3方面的應(yīng)用效果。

1 基本原理

1.1 KKF基本模型

KKF模型是將泛Kriging插值法與Kalman濾波結(jié)合在一起的時(shí)空Kalman濾波模型，其基本模型如式(1)所示

(1)

式中，Zt(s)表示t時(shí)刻位置s處的觀測(cè)值；αt為狀態(tài)量；εt(s)～N(0,R1)，為觀測(cè)噪聲；Φ為狀態(tài)轉(zhuǎn)移矩陣；wt～N(0,Q1)，為狀態(tài)噪聲；H為描述空間相關(guān)性的空間場(chǎng)，由式(2)—式(5)確定[20,26]，H(S)′為H(s)的轉(zhuǎn)置，為便于描述，下文將H(s)簡(jiǎn)記為H

B=Σ-1-Σ-1X(X′Σ-1X)-1X′Σ-1

(2)

式中，Σ為協(xié)方差矩陣；X為趨勢(shì)場(chǎng)；X′為X轉(zhuǎn)置。

為降維運(yùn)算，將式(2)中的B進(jìn)行譜分解，得到

(3)

式中，U=(u1,u2,…,un)，為B的特征向量；U′為U的轉(zhuǎn)置；D=diag(d1,d2,…,dn)，為對(duì)應(yīng)的特征值。

則式(1)中的空間場(chǎng)H可表示為

(4)

式中，σ(s0)為協(xié)方差矢量；Hj表示空間場(chǎng)的前q個(gè)元素為趨勢(shì)場(chǎng)，代表空間的整體變異；Hk表示空間場(chǎng)的p-q個(gè)元素為主要場(chǎng)，代表空間的局部變異。其中，當(dāng)式(5)設(shè)定某個(gè)值時(shí)(如90%或95%)，可以得到p的值

(5)

1.2 STKF基本模型

KKF模型沒有顧及空間小尺度細(xì)微變異，其建模結(jié)果過于光滑[28]。針對(duì)這個(gè)問題，文獻(xiàn)[21]提出了顧及小尺度細(xì)微變異的STKF模型，其模型為

(6)

式中，Zt(s)表示t時(shí)刻位置s處的觀測(cè)值；S(s)為描述空間相關(guān)性的正交空間基；αt為狀態(tài)量；ξ(s,t)～N(m,U2)，表示具有空間相關(guān)性的小尺度細(xì)微變異；εt(s)～N(0,R2)，為觀測(cè)噪聲；Φ為狀態(tài)轉(zhuǎn)移矩陣；J=[S(s)′S(s)]-1S(s)，為狀態(tài)噪聲的系數(shù)矩陣；ηt～N(0,Q2)，為狀態(tài)噪聲。

1.3 STME基本模型

針對(duì)KKF及STKF模型存在海量數(shù)據(jù)計(jì)算困難的問題，文獻(xiàn)[22—23]提出了固定秩的STME模型，其模型為

(7)

式中，Zt(s)表示t時(shí)刻位置s處的觀測(cè)值；μt(s)為趨勢(shì)項(xiàng)；S(s)為描述空間相關(guān)性的空間基；αt為狀態(tài)量；ξt(s)～N(0,U3)，表示小尺度細(xì)微變異；εt(s)～N(0,R3)，為觀測(cè)噪聲；Φ為狀態(tài)轉(zhuǎn)移矩陣；ηt～N(0,Q3)，為狀態(tài)噪聲。

1.4 時(shí)空動(dòng)態(tài)變形建模方法

3種時(shí)空Kalman濾波模型的建模方法類似，均是先基于觀測(cè)數(shù)據(jù)估計(jì)模型未知參數(shù)，然后基于估計(jì)的未知參數(shù)進(jìn)行時(shí)空動(dòng)態(tài)變形分析。本節(jié)以KKF模型為例，介紹利用時(shí)空Kalman濾波進(jìn)行時(shí)空動(dòng)態(tài)變形分析的具體建模方法，其建模方法主要包含構(gòu)建空間場(chǎng)、估計(jì)模型未知參數(shù)及變形分析3部分。

1.4.1 構(gòu)建KKF模型空間場(chǎng)H

由式(4)可知，KKF模型空間場(chǎng)H的前q個(gè)元素為趨勢(shì)場(chǎng)X，一般將測(cè)站的位置坐標(biāo)(如經(jīng)緯度)設(shè)為X的值。選定趨勢(shì)類型(常數(shù)趨勢(shì)、線性趨勢(shì)或二次曲面趨勢(shì))后，利用最小二乘法估計(jì)趨勢(shì)場(chǎng)的系數(shù)矩陣β?？臻g場(chǎng)H剩余的p-q個(gè)元素(主要場(chǎng))可根據(jù)式(2)—式(5)計(jì)算，具體計(jì)算過程如下：

(1) 從原始觀測(cè)數(shù)據(jù)中去除趨勢(shì)項(xiàng)，得到剩余殘差v=Zt(s)-Xβ；利用變差函數(shù)(常用的變差函數(shù)模型有高斯模型、球狀模型、指數(shù)模型)擬合殘差v的空間變異，基于擬合的變差函數(shù)求取所有已知測(cè)站之間的協(xié)方差矩陣Σ。

(2) 根據(jù)式(2)和式(3)獲得矩陣B的特征值及特征向量，將式(3)得到的特征值D=diag(d1,d2,…,dn)進(jìn)行非降序排列，即0=d1=d2=…=dq

(3) 將空間位置s0的位置坐標(biāo)代入擬合的變差函數(shù)即可求得σ(s0)′，并將p的值代入式(4)可求解空間場(chǎng)H的主要場(chǎng)Hk。將求得的趨勢(shì)場(chǎng)Hj和主要場(chǎng)Hk組合即可得到KKF模型的空間場(chǎng)H=[Hj，Hk]。

1.4.2 估計(jì)未知參數(shù)

假定式(1)中狀態(tài)量αt的初始狀態(tài)為α0～N(a0,P0)，由KKF模型的基本原理可知，KKF模型的未知參數(shù)為Ψ=[Φ,R1,Q1]。通過經(jīng)驗(yàn)設(shè)定未知參數(shù)的初始值，利用EM(expectation-maximum)[22]算法迭代估計(jì)未知參數(shù)的最佳估值。

對(duì)于t=1,2,…,T，向前濾波形式如下

αt|t-1=Φαt-1|t-1

(8)

Pt|t-1=ΦPt-1|t-1Φ′+Q1

(9)

K=Pt|t-1H′(HPt|t-1H′+R1)-1

(10)

αt|t=αt|t-1+K(Zt(s)-Hαt|t-1)

(11)

Pt|t=Pt|t-1-KHPt|t-1

(12)

1.4.3 時(shí)空變形分析應(yīng)用

當(dāng)獲得未知參數(shù)最優(yōu)估計(jì)值后，可通過式(13)獲得監(jiān)測(cè)站觀測(cè)數(shù)據(jù)的濾波值、缺失數(shù)據(jù)的插補(bǔ)值及未進(jìn)行實(shí)際觀測(cè)站點(diǎn)(簡(jiǎn)稱非觀測(cè)站點(diǎn))當(dāng)前時(shí)刻的空間預(yù)測(cè)值

(13)

監(jiān)測(cè)站點(diǎn)與非觀測(cè)站點(diǎn)未來時(shí)刻變形數(shù)據(jù)的預(yù)測(cè)值可由式(14)計(jì)算得到

(14)

2 3種時(shí)空Kalman濾波模型比較分析

2.1 理論分析

由式(1)—式(7)可以看出，KKF、STKF及STME模型雖然是3種不同的時(shí)空Kalman濾波模型，其時(shí)空數(shù)據(jù)相關(guān)性描述方式基本相同，具體包括：①均是基于一階矢量自回歸的動(dòng)力學(xué)模型描述時(shí)空數(shù)據(jù)的時(shí)間相關(guān)性；②均是利用空間基函數(shù)描述時(shí)空數(shù)據(jù)的空間相關(guān)性；③均是利用空間基函數(shù)及動(dòng)力學(xué)模型組合形式描述時(shí)空數(shù)據(jù)的時(shí)空相關(guān)性。

同時(shí)也可以看出，3種模型在趨勢(shì)項(xiàng)、小尺度細(xì)微變異、觀測(cè)噪聲、空間基函數(shù)、狀態(tài)轉(zhuǎn)移矩陣計(jì)算及降維方法6方面存在差異，具體有：①KKF和STME模型包含趨勢(shì)項(xiàng)μt(s)，STKF模型不包含趨勢(shì)項(xiàng)，建立STKF模型時(shí)需要先去除數(shù)據(jù)的趨勢(shì)。②KKF模型不包含小尺度細(xì)微變異，將小尺度細(xì)微變異建模在觀測(cè)噪聲中。STKF和STME模型在觀測(cè)方程中單獨(dú)建立小尺度細(xì)微變異ξt(s)。STKF模型的小尺度細(xì)微變異為空間相關(guān)的有色噪聲，STME模型的小尺度細(xì)微變異為空間不相關(guān)的高斯白噪聲。③KKF模型的觀測(cè)噪聲包含小尺度細(xì)微變異及觀測(cè)誤差，STKF和STME模型的觀測(cè)噪聲是高斯白噪聲。④KKF模型的空間基函數(shù)H是根據(jù)式(2)—式(5)構(gòu)建的。STKF和STME模型的空間基函數(shù)S(s)的選擇具有多樣性，需要根據(jù)經(jīng)驗(yàn)選擇合適的空間基函數(shù)。其中，STKF模型的空間基函數(shù)必須是正交的，而STME模型的空間基函數(shù)無正交性限制。⑤STKF模型通過時(shí)空交互函數(shù)來構(gòu)建狀態(tài)轉(zhuǎn)移矩陣Φ，KKF和STME模型的狀態(tài)轉(zhuǎn)移矩陣可根據(jù)具體的應(yīng)用情況靈活選擇，一般情況下，將KKF和STME模型狀態(tài)轉(zhuǎn)移矩陣的初始值設(shè)為單位陣，然后通過EM估計(jì)方法求解其最佳估計(jì)值。⑥KKF模型基于式(5)選取特定的特征值實(shí)現(xiàn)降維功能，STKF及STME模型基于選取的空間基函數(shù)對(duì)時(shí)空數(shù)據(jù)進(jìn)行降維。

2.2 試驗(yàn)分析

由以上理論分析可知，KKF、STKF和STME這3種模型的主要差異在于空間變異的描述形式。為進(jìn)一步分析3種模型對(duì)空間變異描述的精細(xì)程度，筆者模擬了一組空間數(shù)據(jù)，并分別用KKF、STKF和STME這3種模型對(duì)模擬數(shù)據(jù)進(jìn)行空間變異建模擬合，以擬合殘差作為指標(biāo)對(duì)比分析3種模型對(duì)空間變異的表達(dá)能力。

2.2.1 模擬數(shù)據(jù)

模擬的空間變形數(shù)據(jù)由趨勢(shì)項(xiàng)、局部變異項(xiàng)、細(xì)微變異項(xiàng)及噪聲項(xiàng)組成，具體如下。

(1) 趨勢(shì)項(xiàng)：采用一個(gè)與空間位置線性相關(guān)的趨勢(shì)，即μt(s)=0.03·s，其中0 m≤s≤500 m。

(15)

(16)

利用以上方法模擬的空間數(shù)據(jù)如圖1所示。

圖1 模擬的空間變形數(shù)據(jù)Fig.1 The simulated spatial deformation data

2.2.2 KKF模型空間變異擬合效果分析

由圖1可知，模擬數(shù)據(jù)包含與空間位置相關(guān)的線性趨勢(shì)項(xiàng)，因此式(4)中的趨勢(shì)場(chǎng)設(shè)為：Hj=Xj，j=1,2,…,500，其中Xj為位置坐標(biāo)。根據(jù)式(2)和式(3)可計(jì)算式(4)的主要場(chǎng)Hk。為了分析特征值與主要場(chǎng)的關(guān)系，本文對(duì)比分析了第1個(gè)主要場(chǎng)、第200個(gè)主要場(chǎng)及第500個(gè)主要場(chǎng)描述的空間變異與局部變異之間的關(guān)系，如圖2所示。

圖2 模擬數(shù)據(jù)的局部變異及KKF主要場(chǎng)Fig.2 Local variations of simulated data and principal field of KKF model

對(duì)比圖2中的局部變異及第1、200、500個(gè)主要場(chǎng)描述的空間變異，可以看出，KKF建模的第1個(gè)主要場(chǎng)能描述較多的空間變異；第200個(gè)主要場(chǎng)描述的空間變異較第1個(gè)主要場(chǎng)少；第500個(gè)主要場(chǎng)僅描述部分區(qū)域(100～150 m、250～300 m、350～400 m及450～470 m)的空間變異。本文試驗(yàn)的第1個(gè)主要場(chǎng)對(duì)應(yīng)最小的特征值，第500個(gè)主要場(chǎng)對(duì)應(yīng)最大的特征值。通過上述分析可知，較小的特征值能描述大尺度空間變異，較大的特征值描述小尺度空間變異。因此，KKF建模過程中保留的特征值越多，越能有效地描述空間變異。但是當(dāng)數(shù)據(jù)量較大時(shí)，保留的特征值越多，其計(jì)算量越大。為了減輕計(jì)算壓力，通常根據(jù)式(5)選取部分特征值參與建模。本文試驗(yàn)中，若僅建模90%的空間變異(即ratio=0.9)，由式(5)可得至少需要保留451個(gè)特征值，空間擬合及殘差如圖3所示。由圖3可知，KKF能有效地描述空間變異，其擬合殘差的RMS值為0.25 mm。

圖3 KKF空間變異擬合結(jié)果及殘差Fig.3 Spatial modeling results and residuals of KKF model

2.2.3 STKF模型空間變異擬合效果分析

由1.2節(jié)可知，STKF模型利用正交的空間基函數(shù)描述時(shí)空數(shù)據(jù)的空間相關(guān)性，本文試驗(yàn)基于經(jīng)驗(yàn)正交函數(shù)構(gòu)建STKF模型的空間基函數(shù)。式(6)的STKF模型不包含時(shí)空數(shù)據(jù)的趨勢(shì)項(xiàng)，因此在空間建模擬合之前需要去除觀測(cè)數(shù)據(jù)的趨勢(shì)項(xiàng)。采用最小二乘方法去除模擬數(shù)據(jù)的趨勢(shì)項(xiàng)，得到模擬數(shù)據(jù)的局部變異ν?；诰毓烙?jì)求取局部變異ν的經(jīng)驗(yàn)協(xié)方差矩陣Q=E(νν′)-R，其中R表示觀測(cè)噪聲方差。將Q進(jìn)行特征分解：Q=ΨΛΨ′，其中Ψ=[ψ1,ψ2,…,ψn]為特征向量，Λ=diag(λ1,λ2,…,λn)為特征值。本文試驗(yàn)將特征向量ψi按降序排列，并選取前n個(gè)特征向量作為STKF模型的空間基。類似于2.2.2節(jié)，本文對(duì)比分析了前3個(gè)空間基描述的空間變異與局部變異之間的關(guān)系，如圖4所示。

由圖4可以看出，STKF模型第1個(gè)空間基的變化趨勢(shì)與模擬數(shù)據(jù)的局部變異趨勢(shì)非常相似，表明第1個(gè)空間基就可以描述局部變異的整體變化趨勢(shì)。第2個(gè)空間基和第3個(gè)空間基僅描述了部分區(qū)域的空間變異，表明少量的空間基即可描述觀測(cè)數(shù)據(jù)的空間變異，因此STKF模型的降維效果明顯優(yōu)于KKF模型。為了分析STKF模型對(duì)空間變異的建模效果，本文試驗(yàn)選取了前200個(gè)STKF空間基對(duì)模擬數(shù)據(jù)進(jìn)行空間擬合，其擬合效果如圖5所示。

圖4 局部變異及STKF空間基Fig.4 Local variations of simulated data and spatial basis function of STKF model

圖5 STKF模型空間變異擬合結(jié)果及殘差Fig.5 Spatial modeling results and residuals of STKF model

由圖5可知，前200個(gè)STKF空間基能有效地描述本文試驗(yàn)?zāi)M數(shù)據(jù)的空間變異，其擬合殘差的RMS值為0.24 mm。對(duì)比KKF模型的空間擬合效果可以看出，STKF模型基于少量空間基函數(shù)即可獲得優(yōu)于較多空間基函數(shù)的KKF空間變異擬合效果，表明STKF模型的降維效果明顯優(yōu)于KKF模型。

STKF模型的空間基在描述空間變異方面存在一個(gè)問題：最主要的空間基描述了固定的空間變異，然而在實(shí)際應(yīng)用中，空間變異是不斷變化的，導(dǎo)致最主要的空間基無法描述新的空間變異。例如在地表變形過程中，前期數(shù)據(jù)僅有一處區(qū)域存在沉降漏斗(記為區(qū)域A)，而后期數(shù)據(jù)中區(qū)域A的沉降漏斗已無變形，區(qū)域A外的其他區(qū)域又出現(xiàn)一個(gè)新的沉降漏斗(記為區(qū)域B)。在建立STKF模型時(shí)，根據(jù)前期的變形數(shù)據(jù)選取最主要的空間基，由于前期數(shù)據(jù)中區(qū)域B內(nèi)無變形，因此選取的空間基將不會(huì)包含區(qū)域B的變形信息，選定的空間基也無法描述區(qū)域B后期發(fā)生的變形。

2.2.4 STME模型空間變異擬合效果分析

STME模型同樣需要先去除趨勢(shì)項(xiàng)，然后對(duì)去除趨勢(shì)后的大尺度局部變異及小尺度細(xì)微變異進(jìn)行建模。本文試驗(yàn)選取bisquare空間基[25]作為STME模型的空間基，bisquare空間基形式如下

(17)

式中，ck為第k層尺度下的空間基位置；gk為第k層尺度下變異函數(shù)的變程(一般將第k層尺度下最短距離的1.5倍作為變程)。為了準(zhǔn)確描述空間變異，一般從多個(gè)尺度布設(shè)空間基函數(shù)，直至式(18)的空間基擬合殘差趨于平穩(wěn)

V=Zt(s)-μt(s)-S((S′S)-1S′(Zt(s)-

μt(s)))

(18)

本文試驗(yàn)共設(shè)置3個(gè)尺度的空間基，其中第1個(gè)尺度下的空間基位置為：c1=(1+60i)，i=0,1,…,8；第2個(gè)尺度下的空間基位置為：c2=(2+20i)，i=0,1,…,24；第3個(gè)尺度下的空間基位置為：c3=(3+5i)，i=0,1,…,99；共布設(shè)134個(gè)空間基。STME模型空間基擬合效果如圖6所示。

由圖6中可以看出，空間基的擬合效果與原始觀測(cè)數(shù)據(jù)趨勢(shì)一致，能較好地描述原始數(shù)據(jù)的空間變異，但是空間基擬合后的殘差中仍存在較大的空間變異。筆者將殘差中的空間變異建模為式(7)中的小尺度細(xì)微變異ξt(s)，利用球狀模型的變異函數(shù)描述小尺度細(xì)微變異ξt(s)的協(xié)方差矩陣。將求取的小尺度細(xì)微變異ξt(s)加入空間基擬合結(jié)果中，得到STME模型的結(jié)果，如圖7所示。

對(duì)比圖6與圖7可知，加入小尺度細(xì)微變異ξt(s)后，STME模型的空間擬合結(jié)果與原始觀測(cè)數(shù)據(jù)更吻合，其擬合殘差更穩(wěn)定，殘差RMS值為0.21 mm。對(duì)比本文試驗(yàn)3種時(shí)空Kalman濾波模型的擬合殘差可以看出，3種時(shí)空Kalman濾波模型對(duì)空間變異的擬合效果基本相當(dāng)。

圖6 空間基擬合結(jié)果及殘差Fig.6 Spatial modeling results and residuals of spatial basis function

圖7 STME模型空間變異擬合結(jié)果及殘差Fig.7 Spatial modeling results and residuals of STME model

3 3種模型用于時(shí)空變形分析的試驗(yàn)分析

3.1 模擬試驗(yàn)

3.1.1 時(shí)空變形模擬數(shù)據(jù)

模擬的時(shí)空變形數(shù)據(jù)由趨勢(shì)項(xiàng)、局部變異項(xiàng)、細(xì)微變異項(xiàng)及噪聲項(xiàng)組成，其中趨勢(shì)項(xiàng)與2.2.1節(jié)模擬的空間變形監(jiān)測(cè)數(shù)據(jù)的趨勢(shì)一致，在2.2.1節(jié)細(xì)微變異和噪聲項(xiàng)的基礎(chǔ)上分別模擬600期的細(xì)微變異和噪聲項(xiàng)。在2.2.1節(jié)局部變異的基礎(chǔ)上，采用隨機(jī)游走模型模擬時(shí)空變形監(jiān)測(cè)數(shù)據(jù)的局部變異，即ηt=Hηt-1+δt，t=1,2,…,600 d。空間基S乘以每一期的狀態(tài)量得到每一期變形監(jiān)測(cè)數(shù)據(jù)的真值。為了展示模擬的時(shí)空變形監(jiān)測(cè)數(shù)據(jù)，以100 d為時(shí)間間隔均勻選取7期的模擬數(shù)據(jù)，抽取的每期變形數(shù)據(jù)在前一期變形數(shù)據(jù)基礎(chǔ)上加10 mm，如圖8所示。

圖8 模擬的時(shí)空變形數(shù)據(jù)Fig.8 The simulated spatio-temporal deformation data

3.1.2 濾波去噪效果分析

從模擬的時(shí)空變形數(shù)據(jù)中隨機(jī)選取20個(gè)測(cè)站作為已知測(cè)站，利用選取的20個(gè)已知測(cè)站數(shù)據(jù)分別建立KKF、 STKF及STME模型。從20個(gè)已知測(cè)站中隨機(jī)選取9個(gè)測(cè)站分析3種時(shí)空Kalman濾波模型的濾波去噪效果，為對(duì)比分析3種時(shí)空Kalman濾波模型的濾波效果，本文利用普通Kalman濾波對(duì)選取的變形數(shù)據(jù)進(jìn)行濾波去噪。時(shí)空Kalman濾波模型的濾波結(jié)果及普通Kalman濾波結(jié)果如圖9所示，其濾波結(jié)果殘差的RMS統(tǒng)計(jì)見表1。

圖9 時(shí)空Kalman濾波模型及普通Kalman濾波模型的濾波結(jié)果Fig.9 Filtering results of spatio-temporal Kalman filter and Kalman filter

表1中改善率1、改善率2、改善率3分別表示KKF、STKF及STME模型濾波結(jié)果相對(duì)于普通Kalman濾波結(jié)果的改善率。由表1可知，3種時(shí)空Kalman濾波模型均能有效地去除觀測(cè)數(shù)據(jù)的噪聲，其時(shí)空濾波結(jié)果相對(duì)于普通Kalman濾波結(jié)果的平均改善率分別為13.2%、21.1%、28.9%。

表1 濾波結(jié)果殘差RMS統(tǒng)計(jì)

3.1.3 缺失數(shù)據(jù)插補(bǔ)效果分析

將上一節(jié)選取的9個(gè)已知測(cè)站數(shù)據(jù)連續(xù)缺失100 d，然后利用20個(gè)已知測(cè)站數(shù)據(jù)分別建立KKF、STKF及STME模型，并對(duì)選取的9個(gè)已知測(cè)站缺失的100 d數(shù)據(jù)進(jìn)行時(shí)空插值。為對(duì)比分析3種時(shí)空Kalman濾波模型的時(shí)空插值效果，筆者利用分段三次Hermite插值多項(xiàng)式(PCHIP)插值出缺失的變形數(shù)據(jù)。時(shí)空插值結(jié)果及PCHIP插值結(jié)果如圖10所示，其插值結(jié)果殘差的RMS統(tǒng)計(jì)見表2。

圖10 時(shí)空Kalman濾波模型時(shí)空插值及PCHIP時(shí)間插值結(jié)果Fig.10 The interpolation results of spatio-temporal Kalman and PCHIP

圖10中每幅圖內(nèi)灰色豎虛線表示缺失數(shù)據(jù)的時(shí)間段，表2的改善率1、改善率2、改善率3分別表示KKF、STKF及STME模型時(shí)空插值結(jié)果相對(duì)于PCHIP插值結(jié)果的改善率。由圖10及表2可知，3種時(shí)空Kalman濾波模型均能有效地插補(bǔ)缺失的觀測(cè)數(shù)據(jù)，其時(shí)空插值結(jié)果相對(duì)于PCHIP插值結(jié)果的平均改善率分別為61.0%、32.2%、33.9%。

表2 時(shí)空Kalman濾波模型時(shí)空插值及PCHIP時(shí)間插值結(jié)果殘差RMS統(tǒng)計(jì)

3.1.4 空間預(yù)測(cè)效果分析

時(shí)空Kalman濾波既可以基于已知測(cè)站變形數(shù)據(jù)的時(shí)空相關(guān)性預(yù)測(cè)非觀測(cè)站點(diǎn)當(dāng)前時(shí)刻的變形數(shù)據(jù)，也可以預(yù)測(cè)變形體的未來發(fā)展趨勢(shì)。為了便于區(qū)分，本文將非觀測(cè)站點(diǎn)當(dāng)前時(shí)刻變形數(shù)據(jù)的時(shí)空預(yù)測(cè)記為空間預(yù)測(cè)，將未來時(shí)刻變形數(shù)據(jù)的時(shí)空預(yù)測(cè)記為時(shí)空預(yù)測(cè)。本節(jié)基于模擬的時(shí)空變形數(shù)據(jù)檢驗(yàn)3種時(shí)空Kalman濾波模型的空間預(yù)測(cè)性能。

利用20個(gè)已知測(cè)站DOY 1—DOY 500的數(shù)據(jù)分別建立KKF、STKF及STME模型，從未參與建模的480個(gè)測(cè)站中隨機(jī)選取9個(gè)測(cè)站作為驗(yàn)證測(cè)站，基于建立的時(shí)空Kalman濾波模型預(yù)測(cè)9個(gè)驗(yàn)證測(cè)站DOY 1—DOY 500的變形數(shù)據(jù)，為對(duì)比分析3種時(shí)空Kalman濾波模型的空間預(yù)測(cè)效果，筆者利用普通Kriging模型插值出9個(gè)驗(yàn)證測(cè)站的變形數(shù)據(jù)。3種時(shí)空Kalman濾波模型空間預(yù)測(cè)結(jié)果及普通Kriging模型空間插值結(jié)果如圖11所示，其預(yù)測(cè)結(jié)果殘差的RMS統(tǒng)計(jì)見表3。

圖11 時(shí)空Kalman濾波模型空間預(yù)測(cè)及Kriging空間插值結(jié)果Fig.11 Results of spatio-temporal Kalman and Kriging

表3 時(shí)空Kalman濾波模型空間預(yù)測(cè)及Kriging空間插值結(jié)果殘差RMS統(tǒng)計(jì)

表3的改善率1、改善率2、改善率3分別表示KKF、STKF、STME模型空間預(yù)測(cè)結(jié)果相對(duì)于Kriging空間插值結(jié)果的改善率。由圖11及表3可知，3種時(shí)空Kalman濾波模型均能有效地預(yù)測(cè)非觀測(cè)站點(diǎn)當(dāng)前時(shí)刻的變形數(shù)據(jù)，其空間預(yù)測(cè)結(jié)果相對(duì)于Kriging空間插值結(jié)果的平均改善率分別為66.7%、64.0%、66.0%。

3.1.5 時(shí)空預(yù)測(cè)效果分析

時(shí)空Kalman濾波既可以預(yù)測(cè)已知測(cè)站未來時(shí)刻的變形數(shù)據(jù)，也可以預(yù)測(cè)非觀測(cè)站點(diǎn)未來時(shí)刻的變形數(shù)據(jù)。本節(jié)分別利用時(shí)空Kalman濾波預(yù)測(cè)已知測(cè)站和非觀測(cè)站點(diǎn)未來時(shí)刻的變形數(shù)據(jù)，并對(duì)比分析其時(shí)空預(yù)測(cè)精度。本文僅對(duì)比分析了單步預(yù)測(cè)結(jié)果，即利用前期變形數(shù)據(jù)僅預(yù)測(cè)下一天的變形，然后基于新觀測(cè)的變形數(shù)據(jù)更新模型參數(shù)，繼續(xù)預(yù)測(cè)下一天的變形。

利用20個(gè)已知測(cè)站DOY 1—DOY 500的變形數(shù)據(jù)分別建立KKF、STKF及STME模型，從未參與建模的480個(gè)測(cè)站中隨機(jī)選取5個(gè)測(cè)站作為驗(yàn)證測(cè)站(非觀測(cè)站點(diǎn))，基于建立的時(shí)空Kalman濾波模型對(duì)20個(gè)已知測(cè)站及5個(gè)驗(yàn)證測(cè)站未來100 d(DOY 501—DOY 600)的變形數(shù)據(jù)進(jìn)行時(shí)空預(yù)測(cè)。為對(duì)比分析3種時(shí)空Kalman濾波模型對(duì)未來變形數(shù)據(jù)的時(shí)空預(yù)測(cè)效果，筆者利用普通Kalman濾波對(duì)20個(gè)已知測(cè)站未來100 d的變形數(shù)據(jù)進(jìn)行時(shí)間預(yù)測(cè)，利用Kalman濾波+Kriging組合模型對(duì)5個(gè)驗(yàn)證測(cè)站(非觀測(cè)站點(diǎn))未來100 d的變形數(shù)據(jù)進(jìn)行時(shí)間和空間預(yù)測(cè)。從20個(gè)已知測(cè)站中隨機(jī)選取4個(gè)測(cè)站的預(yù)測(cè)結(jié)果進(jìn)行對(duì)比分析，4個(gè)已知測(cè)站及5個(gè)驗(yàn)證測(cè)站(非觀測(cè)站點(diǎn))未來100 d的預(yù)測(cè)結(jié)果如圖12所示。4個(gè)已知測(cè)站預(yù)測(cè)結(jié)果殘差的RMS統(tǒng)計(jì)見表4，5個(gè)驗(yàn)證測(cè)站預(yù)測(cè)結(jié)果殘差的RMS統(tǒng)計(jì)見表5。

圖12 已知測(cè)站及非觀測(cè)站點(diǎn)未來100天變形數(shù)據(jù)的預(yù)測(cè)結(jié)果Fig.12 Prediction results of known and unknown stations for the next 100 days

表4 已知測(cè)站時(shí)空Kalman濾波模型時(shí)空預(yù)測(cè)及Kalman濾波模型時(shí)間預(yù)測(cè)結(jié)果殘差RMS統(tǒng)計(jì)

表5 非觀測(cè)站點(diǎn)時(shí)空Kalman濾波模型時(shí)空預(yù)測(cè)及Kalman+Kriging組合模型預(yù)測(cè)結(jié)果殘差RMS統(tǒng)計(jì)

圖12(a)—(d)表示4個(gè)已知測(cè)站的時(shí)空Kalman濾波模型時(shí)空預(yù)測(cè)結(jié)果和普通Kalman濾波時(shí)間預(yù)測(cè)結(jié)果，圖12(e)—(i)表示5個(gè)驗(yàn)證測(cè)站的時(shí)空Kalman濾波模型時(shí)空預(yù)測(cè)結(jié)果和Kalman濾波+Kriging組合模型預(yù)測(cè)結(jié)果。表4的改善率1、改善率2、改善率3分別表示KKF、STKF及STME模型時(shí)空預(yù)測(cè)結(jié)果相對(duì)于普通Kalman濾波時(shí)間預(yù)測(cè)結(jié)果的改善率。由圖12可以看出，3種時(shí)空Kalman濾波模型預(yù)測(cè)的已知測(cè)站未來變形與真實(shí)變形的變化趨勢(shì)一致。由表4的RMS統(tǒng)計(jì)可知，3種時(shí)空Kalman濾波模型均能有效地基于已知測(cè)站時(shí)空變形數(shù)據(jù)預(yù)測(cè)已知測(cè)站未來時(shí)刻的變形數(shù)據(jù)，其時(shí)空預(yù)測(cè)結(jié)果相對(duì)于普通Kalman濾波時(shí)間預(yù)測(cè)結(jié)果的平均改善率分別為23.8%、12.7%、27.0%。

表5的改善率1、改善率2、改善率3分別表示KKF、STKF及STME模型時(shí)空預(yù)測(cè)結(jié)果相對(duì)于Kalman濾波+Kriging組合模型預(yù)測(cè)結(jié)果的改善率。由表5及圖12可知，KKF、STKF及STME模型對(duì)非觀測(cè)站點(diǎn)未來變形的時(shí)空預(yù)測(cè)結(jié)果與真實(shí)變形吻合較好，其時(shí)空預(yù)測(cè)結(jié)果相對(duì)于Kalman濾波+Kriging組合模型預(yù)測(cè)結(jié)果的平均改善率分別為26.0%、12.3%、24.7%。對(duì)比表4和表5時(shí)空預(yù)測(cè)結(jié)果殘差的平均RMS值可以看出，非觀測(cè)站點(diǎn)時(shí)空預(yù)測(cè)結(jié)果的平均RMS值略高于已知測(cè)站時(shí)空預(yù)測(cè)結(jié)果的平均RMS值。其主要原因是已知測(cè)站通過觀測(cè)數(shù)據(jù)更新建立的時(shí)空Kalman濾波模型，使其與實(shí)測(cè)數(shù)據(jù)更吻合，而非觀測(cè)站點(diǎn)的時(shí)空預(yù)測(cè)結(jié)果是基于建立的時(shí)空Kalman濾波模型進(jìn)行時(shí)空預(yù)測(cè)的，無實(shí)際觀測(cè)數(shù)據(jù)對(duì)模型預(yù)測(cè)數(shù)據(jù)進(jìn)行修正。

3.2 南加州地區(qū)GPS地面沉降數(shù)據(jù)試驗(yàn)

為進(jìn)一步對(duì)比分析3種時(shí)空Kalman濾波模型在變形分析中的應(yīng)用效果，本節(jié)選取了南加州地區(qū)GPS地面沉降數(shù)據(jù)進(jìn)行處理與分析。從SCIGN網(wǎng)站(http:∥www.scign.org/)下載了2005年1月1日—2007年1月1日共731 d的南加州綜合GPS監(jiān)測(cè)網(wǎng)中209個(gè)測(cè)站的沉降數(shù)據(jù)，選取數(shù)據(jù)缺失較少的4個(gè)GPS站點(diǎn)(7ODM、BKMS、BSRY和HOGS)作為驗(yàn)證測(cè)站，其測(cè)站位置分布如圖13所示。

圖13 南加州GPS測(cè)站位置分布Fig.13 Location distribution of GPS stations in Southern California

由于實(shí)際的變形數(shù)據(jù)含有觀測(cè)噪聲，無法定量分析3種時(shí)空Kalman濾波模型的濾波去噪精度，因此本節(jié)基于南加州地區(qū)GPS地面沉降數(shù)據(jù)定量地對(duì)比分析3種時(shí)空Kalman濾波模型在缺失數(shù)據(jù)插補(bǔ)、空間預(yù)測(cè)及時(shí)空預(yù)測(cè)3方面的應(yīng)用效果。

將4個(gè)驗(yàn)證測(cè)站連續(xù)刪除200 d的沉降數(shù)據(jù)模擬大量缺失的GPS變形序列，利用PCHIP、KKF、STKF及STME模型插補(bǔ)缺失的數(shù)據(jù)，其插值結(jié)果殘差的RMS統(tǒng)計(jì)見表6。利用205個(gè)GPS測(cè)站DOY 1—DOY 731的沉降數(shù)據(jù)分別建立Kriging、KKF、STKF及STME模型，并預(yù)測(cè)4個(gè)驗(yàn)證測(cè)站DOY 1—DOY 731的沉降數(shù)據(jù)，其預(yù)測(cè)結(jié)果殘差的RMS統(tǒng)計(jì)見表7。利用205個(gè)GPS測(cè)站DOY 1—DOY 631的沉降數(shù)據(jù)分別建立Kalman濾波+Kriging組合模型、KKF、STKF及STME模型，并預(yù)測(cè)4個(gè)驗(yàn)證測(cè)站未來100 d(DOY 632—DOY 731)的沉降數(shù)據(jù)，其預(yù)測(cè)結(jié)果殘差的RMS統(tǒng)計(jì)見表8。

由表6、表7及表8可知，3種時(shí)空Kalman濾波模型對(duì)于實(shí)測(cè)變形數(shù)據(jù)的缺失數(shù)據(jù)插補(bǔ)精度、空間預(yù)測(cè)精度及時(shí)空預(yù)測(cè)精度均優(yōu)于普通模型的插補(bǔ)及預(yù)測(cè)精度，進(jìn)一步檢驗(yàn)了3種時(shí)空Kalman濾波在變形分析應(yīng)用中的有效性。

表6 時(shí)空Kalman濾波模型時(shí)空插值及PCHIP時(shí)間插值結(jié)果殘差RMS統(tǒng)計(jì)

表7 時(shí)空Kalman濾波模型空間預(yù)測(cè)及Kriging空間插值結(jié)果殘差RMS統(tǒng)計(jì)

表8 時(shí)空Kalman濾波模型空間預(yù)測(cè)及Kalman+Kriging組合模型預(yù)測(cè)結(jié)果殘差RMS統(tǒng)計(jì)

4 結(jié) 論

本文從基本原理和試驗(yàn)兩方面對(duì)比分析了3種典型時(shí)空Kalman濾波模型的異同點(diǎn)及其在變形分析中的適用性和應(yīng)用效果。理論上，3種時(shí)空Kalman濾波模型均基于空間基函數(shù)及動(dòng)力學(xué)模型組合形式描述時(shí)空數(shù)據(jù)的時(shí)空相關(guān)性，其主要差異在于空間變異的描述形式，如空間數(shù)據(jù)的趨勢(shì)項(xiàng)、小尺度細(xì)微變異、觀測(cè)噪聲及空間基函數(shù)。在時(shí)空變形分析應(yīng)用效果方面，3種時(shí)空Kalman濾波模型均具有較高精度的時(shí)空濾波、時(shí)空插值及時(shí)空預(yù)測(cè)性能：KKF模型降維效果不明顯，更適合于監(jiān)測(cè)點(diǎn)稀疏(如GNSS監(jiān)測(cè))的時(shí)空變形分析應(yīng)用；STKF模型和STME模型降維效果明顯，更適合于海量測(cè)站(如InSAR監(jiān)測(cè)或GNSS與InSAR集成監(jiān)測(cè)等)的時(shí)空變形分析應(yīng)用。但STKF模型的空間基描述了固定的空間變異，無法適用于復(fù)雜變形規(guī)律的變形分析，因此優(yōu)先選用STME模型用于海量測(cè)站的時(shí)空變形分析。