2. 江苏海洋大学海洋技术与测绘学院, 江苏 连云港 222005;
3. 湖南省精密工程测量与形变灾害监测重点实验室, 湖南 长沙 410083
2. School of Marine Technology and Geomatics, Jiangsu Ocean University, Lianyungang 222005, China;
3. Key Laboratory of Precise Engineering Surveying & Deformation Disaster Monitoring of Hunan Province, Changsha 410083, China
随着传感器网络、GNSS、InSAR等现代变形监测技术的发展，人们获取了大量的变形数据[1-3]，这些变形数据具有较为复杂的时空相关性，传统的空间和时间序列分析方法已无法满足时空变形分析的需求，因此一些学者开始将时空数据处理的理论与方法(例如时空地统计学模型、时空自回归模型、时空地理加权回归模型及时空Kalman滤波等)引入变形分析领域[4]。时空地统计学模型是在传统地统计学的基础上，运用时空插值及随机模拟估算和量化时空变化，动态反映地理空间的变化过程，其代表性的方法是时空克里金，在边坡[5-6]及隧道围岩[7]变形分析中得到了较好的应用。时空克里金模型虽然同时考虑了变形数据的时间和空间相关性，但是它计算量较大，不能满足变形监测大数据时空分析的需求。时空自回归模型是根据空间权矩阵及AR模型描述时空数据的时空相关性，已用于大坝[8-9]、桥梁[10]、地下管线[11]及建筑物[12]时空变形预测，且具有较高的预测精度，但是该模型分别考虑空间和时间上的相关性，并不是真正意义上的时空分析模型。时空地理加权回归模型是在普通地理加权回归模型的基础上进行时空扩展，构建了顾及时空相关性的时空回归模型，它已被用于大坝变形分析[13-14]，但是该模型需要对监测点的位移与影响因素建立多元变系数回归关系，计算量也较大。时空Kalman滤波是一种融入空间模型的Kalman滤波方法，在时空数据分析方面具有很大的优越性，能对时空数据进行时空动态滤波、时空插值及时空预测，已被广泛应用于众多时空数据分析领域，如传染病动态演变[15]、遥感影像融合[16]、空气污染物建模[17]、气候变化[18]等领域。文献[19]提出了第一个具有统计意义的时空Kalman滤波，发展形成了3种代表性的时空Kalman滤波模型：Kriged Kalman filter(KKF)模型[20]、space time Kalman filter(STKF)模型[21]、spatio-temporal mixed effects (STME)[22-23]模型。文献[24—27]将时空Kalman滤波引入变形监测领域，并在动态滤波插补GPS地面沉降数据、融合GNSS和InSAR沉降数据获取高时空分辨率地表沉降数据及大坝变形的时空预测方面取得了良好的应用效果。但上述研究仅对时空Kalman滤波用于变形分析进行了初步探索，未对已有典型的时空Kalman滤波模型的理论及其在变形分析中的适用性进行深入分析。
为了更好地利用时空Kalman滤波进行时空变形分析，本文从基本原理及试验两方面对比分析KKF、STKF及STME 3种典型时空Kalman滤波模型的异同点及其在变形分析应用中的适用性，并对比分析其在时空滤波去噪、数据插补及变形预测3方面的应用效果。
1 基本原理 1.1 KKF基本模型KKF模型是将泛Kriging插值法与Kalman滤波结合在一起的时空Kalman滤波模型，其基本模型如式(1)所示
(1)
式中，Zt(s)表示t时刻位置s处的观测值；αt为状态量；εt(s)~N(0, R1)，为观测噪声；Φ为状态转移矩阵；wt~N(0, Q1)，为状态噪声；H为描述空间相关性的空间场，由式(2)—式(5)确定[20, 26]，H(S)′为H(s)的转置，为便于描述，下文将H(s)简记为H
(2)
式中，Σ为协方差矩阵；X为趋势场；X′为X转置。
为降维运算，将式(2)中的B进行谱分解，得到
(3)
式中，U=(u1, u2, …, un)，为B的特征向量；U′为U的转置；D=diag(d1, d2, …, dn)，为对应的特征值。
则式(1)中的空间场H可表示为
(4)
式中，σ(s0)为协方差矢量；Hj表示空间场的前q个元素为趋势场，代表空间的整体变异；Hk表示空间场的p-q个元素为主要场，代表空间的局部变异。其中，当式(5)设定某个值时(如90%或95%)，可以得到p的值
(5)
KKF模型没有顾及空间小尺度细微变异，其建模结果过于光滑[28]。针对这个问题，文献[21]提出了顾及小尺度细微变异的STKF模型，其模型为
(6)
式中，Zt(s)表示t时刻位置s处的观测值；S(s)为描述空间相关性的正交空间基；αt为状态量；ξ(s, t)~N(m, U2)，表示具有空间相关性的小尺度细微变异；εt(s)~N(0, R2)，为观测噪声；Φ为状态转移矩阵；J=[S(s)′ S(s)]-1S(s)，为状态噪声的系数矩阵；ηt~N(0, Q2)，为状态噪声。
1.3 STME基本模型针对KKF及STKF模型存在海量数据计算困难的问题，文献[22—23]提出了固定秩的STME模型，其模型为
(7)
式中，Zt(s)表示t时刻位置s处的观测值；μt(s)为趋势项；S(s)为描述空间相关性的空间基；αt为状态量；ξt(s)~N(0, U3)，表示小尺度细微变异；εt(s)~N(0, R3)，为观测噪声；Φ为状态转移矩阵；ηt~N(0, Q3)，为状态噪声。
1.4 时空动态变形建模方法3种时空Kalman滤波模型的建模方法类似，均是先基于观测数据估计模型未知参数，然后基于估计的未知参数进行时空动态变形分析。本节以KKF模型为例，介绍利用时空Kalman滤波进行时空动态变形分析的具体建模方法，其建模方法主要包含构建空间场、估计模型未知参数及变形分析3部分。
1.4.1 构建KKF模型空间场H由式(4)可知，KKF模型空间场H的前q个元素为趋势场X，一般将测站的位置坐标(如经纬度)设为X的值。选定趋势类型(常数趋势、线性趋势或二次曲面趋势)后，利用最小二乘法估计趋势场的系数矩阵β。空间场H剩余的p-q个元素(主要场)可根据式(2)—式(5)计算，具体计算过程如下：
(1) 从原始观测数据中去除趋势项，得到剩余残差v=Zt(s)-Xβ；利用变差函数(常用的变差函数模型有高斯模型、球状模型、指数模型)拟合残差v的空间变异，基于拟合的变差函数求取所有已知测站之间的协方差矩阵Σ。
(2) 根据式(2)和式(3)获得矩阵B的特征值及特征向量，将式(3)得到的特征值D=diag(d1, d2, …, dn) 进行非降序排列，即0=d1=d2=…= dq < dq+1 < …dn。由式(2)可知BX =0，即前q个特征值为0，若式(5)的ratio取值为90%，可求取p的值。
(3) 将空间位置s0的位置坐标代入拟合的变差函数即可求得σ(s0)′，并将p的值代入式(4)可求解空间场H的主要场Hk。将求得的趋势场Hj和主要场Hk组合即可得到KKF模型的空间场H=[Hj，Hk]。
1.4.2 估计未知参数假定式(1)中状态量αt的初始状态为α0~N(a0, P0)，由KKF模型的基本原理可知，KKF模型的未知参数为Ψ=[Φ, R1, Q1]。通过经验设定未知参数的初始值，利用EM(expectation-maximum)[22]算法迭代估计未知参数的最佳估值。
对于t=1, 2, …, T，向前滤波形式如下
(8)
(9)
(10)
(11)
(12)
当获得未知参数最优估计值后，可通过式(13)获得监测站观测数据的滤波值、缺失数据的插补值及未进行实际观测站点(简称非观测站点)当前时刻的空间预测值
(13)
实际应用中，当式(13)中H(s)表示有观测数据的监测站空间场时，
监测站点与非观测站点未来时刻变形数据的预测值可由式(14)计算得到
(14)
式中，
由式(1)—式(7)可以看出，KKF、STKF及STME模型虽然是3种不同的时空Kalman滤波模型，其时空数据相关性描述方式基本相同，具体包括：①均是基于一阶矢量自回归的动力学模型描述时空数据的时间相关性；②均是利用空间基函数描述时空数据的空间相关性；③均是利用空间基函数及动力学模型组合形式描述时空数据的时空相关性。
同时也可以看出，3种模型在趋势项、小尺度细微变异、观测噪声、空间基函数、状态转移矩阵计算及降维方法6方面存在差异，具体有：①KKF和STME模型包含趋势项μt(s)，STKF模型不包含趋势项，建立STKF模型时需要先去除数据的趋势。②KKF模型不包含小尺度细微变异，将小尺度细微变异建模在观测噪声中。STKF和STME模型在观测方程中单独建立小尺度细微变异ξt(s)。STKF模型的小尺度细微变异为空间相关的有色噪声，STME模型的小尺度细微变异为空间不相关的高斯白噪声。③KKF模型的观测噪声包含小尺度细微变异及观测误差，STKF和STME模型的观测噪声是高斯白噪声。④KKF模型的空间基函数H是根据式(2)—式(5)构建的。STKF和STME模型的空间基函数S(s) 的选择具有多样性，需要根据经验选择合适的空间基函数。其中，STKF模型的空间基函数必须是正交的，而STME模型的空间基函数无正交性限制。⑤STKF模型通过时空交互函数来构建状态转移矩阵Φ，KKF和STME模型的状态转移矩阵可根据具体的应用情况灵活选择，一般情况下，将KKF和STME模型状态转移矩阵的初始值设为单位阵，然后通过EM估计方法求解其最佳估计值。⑥KKF模型基于式(5)选取特定的特征值实现降维功能，STKF及STME模型基于选取的空间基函数对时空数据进行降维。
2.2 试验分析由以上理论分析可知，KKF、STKF和STME这3种模型的主要差异在于空间变异的描述形式。为进一步分析3种模型对空间变异描述的精细程度，笔者模拟了一组空间数据，并分别用KKF、STKF和STME这3种模型对模拟数据进行空间变异建模拟合，以拟合残差作为指标对比分析3种模型对空间变异的表达能力。
2.2.1 模拟数据模拟的空间变形数据由趋势项、局部变异项、细微变异项及噪声项组成，具体如下。
(1) 趋势项：采用一个与空间位置线性相关的趋势，即μt(s)=0.03·s，其中0 m≤s≤500 m。
(2) 局部变异项：构造满足ν= Sη，ν~N(0, K)的局部变异，其中K建模为exp(-||si-sj||/45)的指数分布；S为在一维区域D={s|s=1, 2, …, 500}内从3个尺度布设的bisquare空间基函数。其中，第1尺度空间基的中心位置为cx1={x1: x1=50+60k, k=0, 1, …, 7}；第2尺度空间基的中心位置为cx2={x2: x2=30+30k, k=0, 1, …, 15}；第3尺度空间基的中心位置为cx3={x3: x3=20+15k, k=0, 1, …, 32}，共布设57个空间基。
(3) 细微变异项：通过细微变异与总变异的比例关系获得细微变异的方差。本文设定FVP=0.05，利用式(15)求得细微变异方差σξ2，进而构造满足ξ~N(0, σξ2)的细微变异
(15)
式中，n为观测点的个数；σ2ξ为细微变异方差；S为选取的空间基函数；K为局部变异方差；S′为S的转置。
(4) 噪声项：基于信噪比公式获得观测噪声的方差。本文设定SNR=10，利用式(16)求得观测噪声方差σε2，进而构造满足ε~N(0, σε2)的观测噪声
(16)
利用以上方法模拟的空间数据如图 1所示。
|
| 图 1 模拟的空间变形数据 Fig. 1 The simulated spatial deformation data |
2.2.2 KKF模型空间变异拟合效果分析
由图 1可知，模拟数据包含与空间位置相关的线性趋势项，因此式(4)中的趋势场设为：Hj= Xj，j=1, 2, …, 500，其中Xj为位置坐标。根据式(2)和式(3)可计算式(4)的主要场Hk。为了分析特征值与主要场的关系，本文对比分析了第1个主要场、第200个主要场及第500个主要场描述的空间变异与局部变异之间的关系，如图 2所示。
|
| 图 2 模拟数据的局部变异及KKF主要场 Fig. 2 Local variations of simulated data and principal field of KKF model |
对比图 2中的局部变异及第1、200、500个主要场描述的空间变异，可以看出，KKF建模的第1个主要场能描述较多的空间变异；第200个主要场描述的空间变异较第1个主要场少；第500个主要场仅描述部分区域(100~150 m、250~300 m、350~400 m及450~470 m)的空间变异。本文试验的第1个主要场对应最小的特征值，第500个主要场对应最大的特征值。通过上述分析可知，较小的特征值能描述大尺度空间变异，较大的特征值描述小尺度空间变异。因此，KKF建模过程中保留的特征值越多，越能有效地描述空间变异。但是当数据量较大时，保留的特征值越多，其计算量越大。为了减轻计算压力，通常根据式(5)选取部分特征值参与建模。本文试验中，若仅建模90%的空间变异(即ratio=0.9)，由式(5)可得至少需要保留451个特征值，空间拟合及残差如图 3所示。由图 3可知，KKF能有效地描述空间变异，其拟合残差的RMS值为0.25 mm。
|
| 图 3 KKF空间变异拟合结果及残差 Fig. 3 Spatial modeling results and residuals of KKF model |
2.2.3 STKF模型空间变异拟合效果分析
由1.2节可知，STKF模型利用正交的空间基函数描述时空数据的空间相关性，本文试验基于经验正交函数构建STKF模型的空间基函数。式(6)的STKF模型不包含时空数据的趋势项，因此在空间建模拟合之前需要去除观测数据的趋势项。采用最小二乘方法去除模拟数据的趋势项，得到模拟数据的局部变异ν。基于矩估计求取局部变异ν的经验协方差矩阵Q=E(νν′)-R，其中R表示观测噪声方差。将Q进行特征分解：Q=ΨΛΨ′，其中Ψ=[ψ1, ψ2, …, ψn]为特征向量，Λ=diag(λ1, λ2, …, λn)为特征值。本文试验将特征向量ψi按降序排列，并选取前n个特征向量作为STKF模型的空间基。类似于2.2.2节，本文对比分析了前3个空间基描述的空间变异与局部变异之间的关系，如图 4所示。
|
| 图 4 局部变异及STKF空间基 Fig. 4 Local variations of simulated data and spatial basis function of STKF model |
由图 4可以看出，STKF模型第1个空间基的变化趋势与模拟数据的局部变异趋势非常相似，表明第1个空间基就可以描述局部变异的整体变化趋势。第2个空间基和第3个空间基仅描述了部分区域的空间变异，表明少量的空间基即可描述观测数据的空间变异，因此STKF模型的降维效果明显优于KKF模型。为了分析STKF模型对空间变异的建模效果，本文试验选取了前200个STKF空间基对模拟数据进行空间拟合，其拟合效果如图 5所示。
|
| 图 5 STKF模型空间变异拟合结果及残差 Fig. 5 Spatial modeling results and residuals of STKF model |
由图 5可知，前200个STKF空间基能有效地描述本文试验模拟数据的空间变异，其拟合残差的RMS值为0.24 mm。对比KKF模型的空间拟合效果可以看出，STKF模型基于少量空间基函数即可获得优于较多空间基函数的KKF空间变异拟合效果，表明STKF模型的降维效果明显优于KKF模型。
STKF模型的空间基在描述空间变异方面存在一个问题：最主要的空间基描述了固定的空间变异，然而在实际应用中，空间变异是不断变化的，导致最主要的空间基无法描述新的空间变异。例如在地表变形过程中，前期数据仅有一处区域存在沉降漏斗(记为区域A)，而后期数据中区域A的沉降漏斗已无变形，区域A外的其他区域又出现一个新的沉降漏斗(记为区域B)。在建立STKF模型时，根据前期的变形数据选取最主要的空间基，由于前期数据中区域B内无变形，因此选取的空间基将不会包含区域B的变形信息，选定的空间基也无法描述区域B后期发生的变形。
2.2.4 STME模型空间变异拟合效果分析STME模型同样需要先去除趋势项，然后对去除趋势后的大尺度局部变异及小尺度细微变异进行建模。本文试验选取bisquare空间基[25]作为STME模型的空间基，bisquare空间基形式如下
(17)
式中，ck为第k层尺度下的空间基位置；gk为第k层尺度下变异函数的变程(一般将第k层尺度下最短距离的1.5倍作为变程)。为了准确描述空间变异，一般从多个尺度布设空间基函数，直至式(18)的空间基拟合残差趋于平稳
(18)
本文试验共设置3个尺度的空间基，其中第1个尺度下的空间基位置为：c1=(1+60i)，i=0, 1, …, 8；第2个尺度下的空间基位置为：c2=(2+20i)，i=0, 1, …, 24；第3个尺度下的空间基位置为：c3=(3+5i)，i=0, 1, …, 99；共布设134个空间基。STME模型空间基拟合效果如图 6所示。
|
| 图 6 空间基拟合结果及残差 Fig. 6 Spatial modeling results and residuals of spatial basis function |
由图 6中可以看出，空间基的拟合效果与原始观测数据趋势一致，能较好地描述原始数据的空间变异，但是空间基拟合后的残差中仍存在较大的空间变异。笔者将残差中的空间变异建模为式(7)中的小尺度细微变异ξt(s)，利用球状模型的变异函数描述小尺度细微变异ξt(s)的协方差矩阵。将求取的小尺度细微变异ξt(s)加入空间基拟合结果中，得到STME模型的结果，如图 7所示。
|
| 图 7 STME模型空间变异拟合结果及残差 Fig. 7 Spatial modeling results and residuals of STME model |
对比图 6与图 7可知，加入小尺度细微变异ξt(s)后，STME模型的空间拟合结果与原始观测数据更吻合，其拟合残差更稳定，残差RMS值为0.21 mm。对比本文试验3种时空Kalman滤波模型的拟合残差可以看出，3种时空Kalman滤波模型对空间变异的拟合效果基本相当。
3 3种模型用于时空变形分析的试验分析 3.1 模拟试验 3.1.1 时空变形模拟数据模拟的时空变形数据由趋势项、局部变异项、细微变异项及噪声项组成，其中趋势项与2.2.1节模拟的空间变形监测数据的趋势一致，在2.2.1节细微变异和噪声项的基础上分别模拟600期的细微变异和噪声项。在2.2.1节局部变异的基础上，采用随机游走模型模拟时空变形监测数据的局部变异，即ηt=Hηt-1+δt，t=1, 2, …, 600 d。空间基S乘以每一期的状态量得到每一期变形监测数据的真值。为了展示模拟的时空变形监测数据，以100 d为时间间隔均匀选取7期的模拟数据，抽取的每期变形数据在前一期变形数据基础上加10 mm，如图 8所示。
|
| 图 8 模拟的时空变形数据 Fig. 8 The simulated spatio-temporal deformation data |
3.1.2 滤波去噪效果分析
从模拟的时空变形数据中随机选取20个测站作为已知测站，利用选取的20个已知测站数据分别建立KKF、STKF及STME模型。从20个已知测站中随机选取9个测站分析3种时空Kalman滤波模型的滤波去噪效果，为对比分析3种时空Kalman滤波模型的滤波效果，本文利用普通Kalman滤波对选取的变形数据进行滤波去噪。时空Kalman滤波模型的滤波结果及普通Kalman滤波结果如图 9所示，其滤波结果残差的RMS统计见表 1。
|
| 图 9 时空Kalman滤波模型及普通Kalman滤波模型的滤波结果 Fig. 9 Filtering results of spatio-temporal Kalman filter and Kalman filter |
| 测站 | 不同滤波模型的RMS/mm | 改善率/(%) | ||||||
| Kalman | KKF | STKF | STME | 改善率1 | 改善率2 | 改善率3 | ||
| 1 | 0.41 | 0.37 | 0.31 | 0.29 | 9.6 | 24.8 | 30.6 | |
| 2 | 0.44 | 0.36 | 0.37 | 0.27 | 19.6 | 17.1 | 39.4 | |
| 3 | 0.40 | 0.34 | 0.32 | 0.27 | 16.3 | 20.1 | 33.2 | |
| 4 | 0.41 | 0.35 | 0.34 | 0.28 | 15.1 | 17.9 | 30.6 | |
| 5 | 0.37 | 0.32 | 0.36 | 0.28 | 13.8 | 2.5 | 25.0 | |
| 6 | 0.34 | 0.30 | 0.28 | 0.28 | 10.2 | 18.0 | 16.6 | |
| 7 | 0.36 | 0.31 | 0.34 | 0.26 | 13.3 | 4.1 | 25.9 | |
| 8 | 0.35 | 0.33 | 0.26 | 0.26 | 4.7 | 25.4 | 24.5 | |
| 9 | 0.36 | 0.32 | 0.15 | 0.29 | 11.1 | 58.0 | 21.6 | |
| 平均值 | 0.38 | 0.33 | 0.30 | 0.27 | 3.2 | 21.1 | 28.9 | |
表 1中改善率1、改善率2、改善率3分别表示KKF、STKF及STME模型滤波结果相对于普通Kalman滤波结果的改善率。由表 1可知，3种时空Kalman滤波模型均能有效地去除观测数据的噪声，其时空滤波结果相对于普通Kalman滤波结果的平均改善率分别为13.2%、21.1%、28.9%。
3.1.3 缺失数据插补效果分析将上一节选取的9个已知测站数据连续缺失100 d，然后利用20个已知测站数据分别建立KKF、STKF及STME模型，并对选取的9个已知测站缺失的100 d数据进行时空插值。为对比分析3种时空Kalman滤波模型的时空插值效果，笔者利用分段三次Hermite插值多项式(PCHIP)插值出缺失的变形数据。时空插值结果及PCHIP插值结果如图 10所示，其插值结果残差的RMS统计见表 2。
|
| 图 10 时空Kalman滤波模型时空插值及PCHIP时间插值结果 Fig. 10 The interpolation results of spatio-temporal Kalman and PCHIP |
| 测站 | 不同滤波模型的RMS/mm | 改善率/(%) | ||||||
| PCHIP | KKF | STKF | STME | 改善率1 | 改善率2 | 改善率3 | ||
| 1 | 1.17 | 0.36 | 0.88 | 0.63 | 69.5 | 24.6 | 46.6 | |
| 2 | 1.40 | 0.29 | 1.17 | 1.08 | 79.0 | 16.5 | 23.0 | |
| 3 | 0.89 | 0.53 | 0.64 | 0.80 | 41.0 | 28.1 | 10.6 | |
| 4 | 0.99 | 0.54 | 0.73 | 0.69 | 45.4 | 26.5 | 30.5 | |
| 5 | 0.91 | 0.59 | 0.76 | 0.63 | 35.5 | 16.5 | 30.8 | |
| 6 | 1.95 | 0.36 | 0.74 | 0.77 | 81.4 | 62.3 | 60.4 | |
| 7 | 1.19 | 0.75 | 0.72 | 0.70 | 37.1 | 39.7 | 41.0 | |
| 8 | 1.10 | 0.57 | 0.86 | 0.95 | 48.3 | 21.9 | 13.4 | |
| 9 | 0.96 | 0.15 | 0.70 | 0.81 | 84.2 | 26.8 | 15.7 | |
| 平均值 | 1.18 | 0.46 | 0.80 | 0.78 | 61.0 | 32.2 | 33.9 | |
图 10中每幅图内灰色竖虚线表示缺失数据的时间段，表 2的改善率1、改善率2、改善率3分别表示KKF、STKF及STME模型时空插值结果相对于PCHIP插值结果的改善率。由图 10及表 2可知，3种时空Kalman滤波模型均能有效地插补缺失的观测数据，其时空插值结果相对于PCHIP插值结果的平均改善率分别为61.0%、32.2%、33.9%。
3.1.4 空间预测效果分析时空Kalman滤波既可以基于已知测站变形数据的时空相关性预测非观测站点当前时刻的变形数据，也可以预测变形体的未来发展趋势。为了便于区分，本文将非观测站点当前时刻变形数据的时空预测记为空间预测，将未来时刻变形数据的时空预测记为时空预测。本节基于模拟的时空变形数据检验3种时空Kalman滤波模型的空间预测性能。
利用20个已知测站DOY 1—DOY 500的数据分别建立KKF、STKF及STME模型，从未参与建模的480个测站中随机选取9个测站作为验证测站，基于建立的时空Kalman滤波模型预测9个验证测站DOY 1—DOY 500的变形数据，为对比分析3种时空Kalman滤波模型的空间预测效果，笔者利用普通Kriging模型插值出9个验证测站的变形数据。3种时空Kalman滤波模型空间预测结果及普通Kriging模型空间插值结果如图 11所示，其预测结果残差的RMS统计见表 3。
|
| 图 11 时空Kalman滤波模型空间预测及Kriging空间插值结果 Fig. 11 Results of spatio-temporal Kalman and Kriging |
| 测站 | 不同滤波模型的RMS/mm | 改善率/(%) | ||||||
| Kriging | KKF | STKF | STME | 改善率1 | 改善率2 | 改善率3 | ||
| 1 | 1.62 | 0.55 | 0.56 | 0.55 | 66.2 | 65.3 | 66.3 | |
| 2 | 1.37 | 0.54 | 0.52 | 0.49 | 60.7 | 61.6 | 63.8 | |
| 3 | 1.56 | 0.54 | 0.55 | 0.51 | 65.6 | 65.0 | 67.3 | |
| 4 | 1.33 | 0.49 | 0.52 | 0.59 | 62.8 | 60.6 | 55.3 | |
| 5 | 1.39 | 0.49 | 0.51 | 0.44 | 64.5 | 63.6 | 68.2 | |
| 6 | 1.49 | 0.54 | 0.57 | 0.56 | 63.6 | 61.6 | 62.2 | |
| 7 | 1.52 | 0.52 | 0.55 | 0.52 | 65.8 | 63.5 | 65.5 | |
| 8 | 1.58 | 0.50 | 0.54 | 0.49 | 68.2 | 65.8 | 68.9 | |
| 9 | 1.68 | 0.45 | 0.50 | 0.42 | 73.2 | 70.2 | 74.9 | |
| 平均值 | 1.50 | 0.50 | 0.54 | 0.51 | 66.7 | 64.0 | 66.0 | |
表 3的改善率1、改善率2、改善率3分别表示KKF、STKF、STME模型空间预测结果相对于Kriging空间插值结果的改善率。由图 11及表 3可知，3种时空Kalman滤波模型均能有效地预测非观测站点当前时刻的变形数据，其空间预测结果相对于Kriging空间插值结果的平均改善率分别为66.7%、64.0%、66.0%。
3.1.5 时空预测效果分析时空Kalman滤波既可以预测已知测站未来时刻的变形数据，也可以预测非观测站点未来时刻的变形数据。本节分别利用时空Kalman滤波预测已知测站和非观测站点未来时刻的变形数据，并对比分析其时空预测精度。本文仅对比分析了单步预测结果，即利用前期变形数据仅预测下一天的变形，然后基于新观测的变形数据更新模型参数，继续预测下一天的变形。
利用20个已知测站DOY 1—DOY 500的变形数据分别建立KKF、STKF及STME模型，从未参与建模的480个测站中随机选取5个测站作为验证测站(非观测站点)，基于建立的时空Kalman滤波模型对20个已知测站及5个验证测站未来100 d(DOY 501—DOY 600)的变形数据进行时空预测。为对比分析3种时空Kalman滤波模型对未来变形数据的时空预测效果，笔者利用普通Kalman滤波对20个已知测站未来100 d的变形数据进行时间预测，利用Kalman滤波+Kriging组合模型对5个验证测站(非观测站点)未来100 d的变形数据进行时间和空间预测。从20个已知测站中随机选取4个测站的预测结果进行对比分析，4个已知测站及5个验证测站(非观测站点)未来100 d的预测结果如图 12所示。4个已知测站预测结果残差的RMS统计见表 4，5个验证测站预测结果残差的RMS统计见表 5。
|
| 图 12 已知测站及非观测站点未来100天变形数据的预测结果 Fig. 12 Prediction results of known and unknown stations for the next 100 days |
| 测站 | 不同滤波模型的RMS/mm | 改善率/(%) | ||||||
| Kalman | KKF | STKF | STME | 改善率1 | 改善率2 | 改善率3 | ||
| 1 | 0.64 | 0.48 | 0.56 | 0.47 | 23.8 | 12.4 | 26.2 | |
| 2 | 0.66 | 0.49 | 0.54 | 0.46 | 25.8 | 18.3 | 30.8 | |
| 3 | 0.58 | 0.45 | 0.55 | 0.43 | 22.1 | 4.0 | 25.3 | |
| 4 | 0.64 | 0.48 | 0.54 | 0.46 | 25.2 | 16.4 | 27.8 | |
| 平均值 | 0.63 | 0.48 | 0.55 | 0.46 | 23.8 | 12.7 | 27.0 | |
| 测站 | 不同滤波模型的RMS/mm | 改善率/(%) | ||||||
| Kalman+ Kriging | KKF | STKF | STME | 改善率1 | 改善率2 | 改善率3 | ||
| 5 | 0.88 | 0.62 | 0.64 | 0.59 | 29.5 | 27.3 | 33.0 | |
| 6 | 0.68 | 0.52 | 0.61 | 0.55 | 23.5 | 10.3 | 19.1 | |
| 7 | 0.69 | 0.55 | 0.64 | 0.59 | 20.3 | 7.2 | 14.5 | |
| 8 | 0.73 | 0.54 | 0.65 | 0.55 | 26.0 | 11.0 | 24.7 | |
| 9 | 0.68 | 0.48 | 0.67 | 0.49 | 29.4 | 1.5 | 27.9 | |
| 平均值 | 0.73 | 0.54 | 0.64 | 0.55 | 26.0 | 12.3 | 24.7 | |
图 12(a)—(d)表示4个已知测站的时空Kalman滤波模型时空预测结果和普通Kalman滤波时间预测结果，图 12(e)—(i)表示5个验证测站的时空Kalman滤波模型时空预测结果和Kalman滤波+Kriging组合模型预测结果。表 4的改善率1、改善率2、改善率3分别表示KKF、STKF及STME模型时空预测结果相对于普通Kalman滤波时间预测结果的改善率。由图 12可以看出，3种时空Kalman滤波模型预测的已知测站未来变形与真实变形的变化趋势一致。由表 4的RMS统计可知，3种时空Kalman滤波模型均能有效地基于已知测站时空变形数据预测已知测站未来时刻的变形数据，其时空预测结果相对于普通Kalman滤波时间预测结果的平均改善率分别为23.8%、12.7%、27.0%。
表 5的改善率1、改善率2、改善率3分别表示KKF、STKF及STME模型时空预测结果相对于Kalman滤波+Kriging组合模型预测结果的改善率。由表 5及图 12可知，KKF、STKF及STME模型对非观测站点未来变形的时空预测结果与真实变形吻合较好，其时空预测结果相对于Kalman滤波+Kriging组合模型预测结果的平均改善率分别为26.0%、12.3%、24.7%。对比表 4和表 5时空预测结果残差的平均RMS值可以看出，非观测站点时空预测结果的平均RMS值略高于已知测站时空预测结果的平均RMS值。其主要原因是已知测站通过观测数据更新建立的时空Kalman滤波模型，使其与实测数据更吻合，而非观测站点的时空预测结果是基于建立的时空Kalman滤波模型进行时空预测的，无实际观测数据对模型预测数据进行修正。
3.2 南加州地区GPS地面沉降数据试验为进一步对比分析3种时空Kalman滤波模型在变形分析中的应用效果，本节选取了南加州地区GPS地面沉降数据进行处理与分析。从SCIGN网站(http://www.scign.org/)下载了2005年1月1日—2007年1月1日共731 d的南加州综合GPS监测网中209个测站的沉降数据，选取数据缺失较少的4个GPS站点(7ODM、BKMS、BSRY和HOGS)作为验证测站，其测站位置分布如图 13所示。
|
| 图 13 南加州GPS测站位置分布 Fig. 13 Location distribution of GPS stations in Southern California |
由于实际的变形数据含有观测噪声，无法定量分析3种时空Kalman滤波模型的滤波去噪精度，因此本节基于南加州地区GPS地面沉降数据定量地对比分析3种时空Kalman滤波模型在缺失数据插补、空间预测及时空预测3方面的应用效果。
将4个验证测站连续删除200 d的沉降数据模拟大量缺失的GPS变形序列，利用PCHIP、KKF、STKF及STME模型插补缺失的数据，其插值结果残差的RMS统计见表 6。利用205个GPS测站DOY 1—DOY 731的沉降数据分别建立Kriging、KKF、STKF及STME模型，并预测4个验证测站DOY 1—DOY 731的沉降数据，其预测结果残差的RMS统计见表 7。利用205个GPS测站DOY 1—DOY 631的沉降数据分别建立Kalman滤波+Kriging组合模型、KKF、STKF及STME模型，并预测4个验证测站未来100 d(DOY 632—DOY 731)的沉降数据，其预测结果残差的RMS统计见表 8。
| 测站 | RMS/mm | 改善率/(%) | ||||||
| PCHIP | KKF | STKF | STME | 改善率1 | 改善率2 | 改善率3 | ||
| 7ODM-U | 12.5 | 10.5 | 11.2 | 10.9 | 16.0 | 10.4 | 12.8 | |
| BKMS-U | 5.4 | 4.2 | 4.5 | 3.9 | 22.2 | 16.7 | 27.8 | |
| BSRY-U | 6.6 | 2.7 | 4.2 | 3.9 | 59.1 | 36.4 | 40.9 | |
| HOGS-U | 8.4 | 3.3 | 4.0 | 2.9 | 60.7 | 52.4 | 65.5 | |
| 平均值 | 8.2 | 5.2 | 6.0 | 5.4 | 36.6 | 26.8 | 34.1 | |
| 测站 | RMS/mm | 改善率/(%) | ||||||
| Kriging | KKF | STKF | STME | 改善率1 | 改善率2 | 改善率3 | ||
| 7ODM-U | 14.1 | 11.9 | 13.2 | 12.0 | 15.6 | 6.4 | 14.9 | |
| BKMS-U | 19.2 | 11.6 | 7.1 | 5.5 | 39.6 | 63.0 | 71.4 | |
| BSRY-U | 5.4 | 3.1 | 4.3 | 3.3 | 42.6 | 20.4 | 38.9 | |
| HOGS-U | 6.0 | 4.1 | 4.9 | 3.6 | 31.7 | 18.3 | 40.0 | |
| 平均值 | 11.2 | 7.7 | 7.4 | 6.1 | 31.3 | 33.9 | 45.5 | |
| 测站 | RMS/mm | 改善率/(%) | ||||||
| Kalman+ Kriging | KKF | STKF | STME | 改善率1 | 改善率2 | 改善率3 | ||
| 7ODM-U | 13.2 | 10.1 | 11.2 | 10.0 | 23.5 | 15.2 | 24.2 | |
| BKMS-U | 7.5 | 4.5 | 5.7 | 4.6 | 40.0 | 24.0 | 38.7 | |
| BSRY-U | 6.4 | 4.6 | 5.1 | 4.8 | 28.1 | 20.3 | 25.0 | |
| HOGS-U | 7.2 | 5.0 | 6.3 | 5.3 | 30.6 | 12.5 | 26.4 | |
| 平均值 | 8.6 | 6.1 | 7.1 | 6.2 | 29.1 | 17.4 | 27.9 | |
由表 6、表 7及表 8可知，3种时空Kalman滤波模型对于实测变形数据的缺失数据插补精度、空间预测精度及时空预测精度均优于普通模型的插补及预测精度，进一步检验了3种时空Kalman滤波在变形分析应用中的有效性。
4 结论本文从基本原理和试验两方面对比分析了3种典型时空Kalman滤波模型的异同点及其在变形分析中的适用性和应用效果。理论上，3种时空Kalman滤波模型均基于空间基函数及动力学模型组合形式描述时空数据的时空相关性，其主要差异在于空间变异的描述形式，如空间数据的趋势项、小尺度细微变异、观测噪声及空间基函数。在时空变形分析应用效果方面，3种时空Kalman滤波模型均具有较高精度的时空滤波、时空插值及时空预测性能：KKF模型降维效果不明显，更适合于监测点稀疏(如GNSS监测)的时空变形分析应用；STKF模型和STME模型降维效果明显，更适合于海量测站(如InSAR监测或GNSS与InSAR集成监测等)的时空变形分析应用。但STKF模型的空间基描述了固定的空间变异，无法适用于复杂变形规律的变形分析，因此优先选用STME模型用于海量测站的时空变形分析。
| [1] |
张勤, 黄观文, 杨成生. 地质灾害监测预警中的精密空间对地观测技术[J]. 测绘学报, 2017, 46(10): 1300-1307. ZHANG Qin, HUANG Guanwen, YANG Chengsheng. Precision space observation technique for geological hazard monitoring and early warning[J]. Acta Geodaetica et Cartographica Sinica, 2017, 46(10): 1300-1307. DOI:10.11947/j.AGCS.2017.20170453 |
| [2] |
许强, 朱星, 李为乐, 等. "天-空-地"协同滑坡监测技术进展[J]. 测绘学报, 2022, 51(7): 1416-1436. XU Qiang, ZHU Xing, LI Weile, et al. Technical progress of space-air-ground collaborative monitoring of landslide[J]. Acta Geodaetica et Cartographica Sinica, 2022, 51(7): 1416-1436. DOI:10.11947/j.AGCS.2022.20220320 |
| [3] |
ZHOU Lü, ZHAO Yizhan, ZHU Zilin, et al. Spatial and temporal evolution of surface subsidence in Tianjin from 2015 to 2020 based on SBAS-InSAR technology[J]. Journal of Geodesy and Geoinformation Science, 2022, 5(1): 60-72. DOI:10.11947/j.JGGS.2022.0107 |
| [4] |
尹晖. 时空变形分析与预报的理论和方法[M]. 北京: 测绘出版社, 2002. YIN Hui. Theory and method of spatio-temporal deformation analysis and prediction[M]. Beijing: Surveying and Mapping Press, 2002. |
| [5] |
王建民, 张锦, 邓增兵, 等. 时空Kriging插值在边坡变形监测中的应用[J]. 煤炭学报, 2014, 39(5): 874-879. WANG Jianmin, ZHANG Jin, DENG Zengbing, et al. Slope deformation analyses with space-time Kriging interpolation method[J]. Journal of China Coal Society, 2014, 39(5): 874-879. |
| [6] |
王建民. 矿山边坡变形监测数据的高斯过程智能分析与预测[D]. 太原: 太原理工大学, 2016. WANG Jianmin. Intelligent analysis and prediction of mine slope deformation monitoring data based on Gaussian process[D]. Taiyuan: Taiyuan University of Technology, 2016. |
| [7] |
张可能, 胡达, 何杰, 等. 基于Kriging时空统一模型的隧道动态施工位移预测[J]. 中南大学学报(自然科学版), 2017, 48(12): 3328-3334. ZHANG Keneng, HU Da, HE Jie, et al. Tunnel construction of dynamic displacement prediction based on unified space-time Kriging model[J]. Journal of Central South University (Science and Technology), 2017, 48(12): 3328-3334. |
| [8] |
李广春, 戴吾蛟, 杨国祥, 等. 时空自回归模型在大坝变形分析中的应用[J]. 武汉大学学报(信息科学版), 2015, 40(7): 877-881. LI Guangchun, DAI Wujiao, YANG Guoxiang, et al. Application of space-time auto-regressive model in dam deformation analysis[J]. Geomatics and Information Science of Wuhan University, 2015, 40(7): 877-881. |
| [9] |
杨志佳, 戴吾蛟, 陈必焰, 等. 克里金时空自回归模型在变形建模中的应用[J]. 测绘科学, 2019, 44(7): 40-45. YANG Zhijia, DAI Wujiao, CHEN Biyan, et al. The application of Kriging's space-time auto-regressive model in deformation modeling[J]. Science of Surveying and Mapping, 2019, 44(7): 40-45. |
| [10] |
孙志鹏, 郭玉平. 基于时空自回归模型的大型桥梁变形监测分析与预报[J]. 全球定位系统, 2015, 40(6): 83-85. SUN Zhipeng, GUO Yuping. Large-scale bridge deformation monitoring analysis and forecasting based on space-time auto-regressive model[J]. GNSS World of China, 2015, 40(6): 83-85. |
| [11] |
柳新强, 王涛. 时空序列模型在地下管线沉降监测中的应用[J]. 北京测绘, 2018, 32(7): 809-813. LIU Xinqiang, WANG Tao. Application of space-time series model in subsidence monitoring of underground pipeline[J]. Beijing Surveying and Mapping, 2018, 32(7): 809-813. |
| [12] |
柳新强, 王涛, 胡泊. 时空序列模型在沉降监测中的应用[J]. 测绘与空间地理信息, 2019, 42(2): 86-89, 93. LIU Xinqiang, WANG Tao, HU Bo. Application of space-time series model in subsidence monitoring[J]. Geomatics & Spatial Information Technology, 2019, 42(2): 86-89, 93. |
| [13] |
YANG Zhijia, DAI Wujiao, SANTERRE R, et al. A spatiotemporal deformation modelling method based on geographically and temporally weighted regression[J]. Mathematical Problems in Engineering, 2019, 4352396. |
| [14] |
YANG Zhijia, DAI Wujiao, YU Wenkun, et al. Mixed geographically and temporally weighted regression for spatio-temporal deformation modelling[J]. Survey Review, 2022, 54(385): 290-300. DOI:10.1080/00396265.2021.1935578 |
| [15] |
ZHUANG Lili, CRESSIE N. Bayesian hierarchical statistical SIRS models[J]. Statistical Methods & Applications, 2014, 23(4): 601-646. |
| [16] |
NGUYEN H, KATZFUSS M, CRESSIE N, et al. Spatio-temporal data fusion for very large remote sensing datasets[J]. Technometrics, 2014, 56(2): 174-185. |
| [17] |
ZAMMIT-MANGION A, CRESSIE N, GANESAN A L, et al. Spatio-temporal bivariate statistical models for atmospheric trace-gas inversion[J]. Chemometrics and Intelligent Laboratory Systems, 2015, 149: 227-241. |
| [18] |
ZHANG Bohai, CRESSIE N. Bayesian inference of spatio-temporal changes of Arctic sea ice[J]. Bayesian Analysis, 2020, 15(2): 605-631. |
| [19] |
HUANG H C, CRESSIE N. Spatio-temporal prediction of snow water equivalent using the Kalman filter[J]. Computational Statistics & Data Analysis, 1996, 22(2): 159-175. |
| [20] |
MARDIA K V, GOODALL C, REDFERN E J, et al. The Kriged Kalman filter[J]. Test, 1998, 7(2): 217-282. |
| [21] |
WIKLE C K, CRESSIE N. A dimension-reduced approach to space-time Kalman filtering[J]. Biometrika, 1999, 86(4): 815-829. |
| [22] |
CRESSIE N, SHI Tao, KANG E L. Fixed rank filtering for spatio-temporal data[J]. Journal of Computational and Graphical Statistics, 2010, 19(3): 724-745. |
| [23] |
KANG E L, CRESSIE N, SHI Tao. Using temporal variability to improve spatial mapping with application to satellite data[J]. Canadian Journal of Statistics, 2010, 38(2): 271-289. |
| [24] |
LIU Ning, DAI Wujiao, SANTERRE R, et al. A MATLAB-based Kriged Kalman filter software for interpolating missing data in GNSS coordinate time series[J]. GPS Solutions, 2018, 22(1): 25. |
| [25] |
LIU Ning, DAI Wujiao, SANTERRE R, et al. High spatio-temporal resolution deformation time series with the fusion of InSAR and GNSS data using spatio-temporal random effect model[J]. IEEE Transactions on Geoscience and Remote Sensing, 2019, 57(1): 364-380. |
| [26] |
DAI Wujiao, LIU Ning, SANTERRE R, et al. Dam deformation monitoring data analysis using space-time Kalman filter[J]. ISPRS International Journal of Geo-Information, 2016, 5(12): 236. |
| [27] |
刘宁, 戴吾蛟, 刘斌. 一种抗差的形变数据插补方法[J]. 测绘科学, 2017, 42(9): 126-131, 190. LIU Ning, DAI Wujiao, LIU Bin. A interpolation method of deformation monitoring data series[J]. Science of Surveying and Mapping, 2017, 42(9): 126-131, 190. |
| [28] |
CRESSIE N A, WIKLE C. Strategies for dynamic space-time statistical modeling: discussion of "the Kriged Kalman filter" by Mardia et al[J]. Test, 1998, 7(2): 257-264. |