CN104699966B - 一种融合GNSS和InSAR数据获取高时空分辨率形变序列的方法 - Google Patents

Abstract

本发明公开了一种融合GNSS和InSAR数据获取高时空分辨率形变序列的方法，根据InSAR数据选择合适的趋势场；利用变差函数模型将去除趋势的残余量进行拟合，得到空间变差函数；依据GPS点位以及变差函数和趋势项构建克里金Kalman滤波中的空间场矩阵；使用EM估计估计克里金Kalman滤波模型中的模型参数；进行克里金Kalman进行滤波；在时空上进行插值得到高时空分辨率的形变序列。与现有的方法相比，本发明的变差函数以及趋势项的选取采用的是具有高空间分辨率的InSAR数据，而不是采用点位稀疏的GPS数据。因此本发明可以有效提高空间场的精度，并且得到高时空分辨率的形变序列。

Description

一种融合GNSS和InSAR数据获取高时空分辨率形变序列的 方法

技术领域

本发明涉及一种融合GNSS和InSAR数据获取高时空分辨率形变序列的方法。

背景技术

GNSS数据具有高时间分辨率低空间分辨率的特点，InSAR数据具有高空间分辨率低时间分辨率的特点。如何融合这两种技术，将其优缺点互补，得到高时空分辨率的形变序列，是亟待解决的问题。克里金Kalman滤波是一种时空Kalman滤波的方法，它不仅考虑了时间上的相关性，也考虑了空间点位的相关性，并且还可以在时空上进行插值得到高时空分辨率的形变序列。而传统的克里金Kalman滤波针对的是单一的GNSS数据。GNSS数据的时间分辨率很高，基本上可以达到1天的分辨率。但是GNSS点位分布稀疏，GNSS基线有的长达几十公里，这并不能较好的体现空间上的相关性。因此直接使用GNSS来构建克里金的空间场是不合理的，需要综合使用具有高空间分辨率的InSAR数据来提高克里金Kalman滤波的空间场精度。

克里金Kalman滤波模型：

其中Z为观测值矢量，H为空间场，α_t为待估状态量，Φ为状态转移矩阵。v为观测噪声，其期望为0，协方差矩阵为R。w为状态噪声，其期望为0，协方差矩阵为Q。

发明内容

本发明所要解决的技术问题是，针对上述现有技术的不足，提供一种融合GNSS和InSAR数据获取高时空分辨率形变序列的方法。

为解决上述技术问题，本发明所采用的技术方案是：一种融合GNSS和InSAR数据获取高时空分辨率形变序列的方法，包括以下步骤：

1)计算InSAR数据的趋势场F_insar，并去除InSAR数据的趋势项；

2)采用曲面拟合的方式，利用所述趋势场F_insar构造曲面拟合趋势项的误差方程，使用最小二乘求得趋势项曲面的系数矩阵，计算去除趋势项的平稳残余项，选择变差函数模型，并利用所述平稳残余项拟合变差函数γ(h)；

3)计算GNSS数据的空间趋势场矩阵F_gnss，利用上述变差函数γ(h)计算GNSS数据点位之间的空间协方差矩阵

4)根据GNSS数据的空间趋势场矩阵F_gnss以及空间协方差矩阵计算弯曲能量矩阵B；

5)利用所述弯曲能量矩阵B以及趋势场矩阵F_gnss构建克里金Kalman滤波的空间场H；

6)使用EM估计计算克里金Kalman滤波中的模型参数，进行克里金Kalman滤波；

7)进行时空插值，得到高时空分辨率的形变序列。

与现有技术相比，本发明所具有的有益效果为：本发明使用克里金Kalman融合GNSS和InSAR数据，其中克里金Kalman空间场中趋势项的选择和变差函数的拟合采用的是具有高空间分辨率的InSAR数据。现有的克里金Kalman采用的是单一的GNSS数据，但由于其点位稀疏的原因并不能得到精度高的空间趋势面以及平稳的残余项。本发明具体首先是采用InSAR数据使用曲面拟合的方式求得空间趋势曲面的系数矩阵，然后再计算平稳残余项，使用变异函数模型拟合残余项的变异函数，这能有效的提高计算空间场的精度，并可以在时空上插值得到高时空分辨率的形变序列。

附图说明

图1为本发明实施方法流程图；

图2(a)为GNSS数据实验变差函数值；图2(b)为InSAR数据实验变差函数值。

具体实施方式

如图1所示，本发明具体步骤为：

1)根据InSAR数据选择合适的趋势场

空间变量往往并不是平稳的，通常含有一个趋势。为了得到一个平稳的序列，需要先去除趋势项。常用的趋势场通常有常数场、线性场以及二次曲面场。若InSAR数据共有m_insar个观测值，即坐标x,y以及观测值Z_insar均为维数为m_insar*1的矢量。则：

(1)若趋势场为常数场，则趋势场F_insar为值全为1且维数为m_insar*1的矢量，即：

F_insar＝[1,1,...,1]′

(2)若趋势场为线性场，则趋势场F_insar为维数为m_insar*3的矩阵，如下：

F_insar＝[x′；y′；1]′

(3)若趋势场为二次曲面场，则趋势场F_insar为m_insar*6的矩阵，如下：

F_insar＝[(x²)′；(y²)′；(xy)′；x′；y′；1]′

2)去除趋势并拟合残余的空间协变差函数模型

(1)利用步骤1)中选择的趋势场拟合趋势项，即根据选择的趋势场，构造拟合趋势项的误差方程：

v＝F_insarβ-Z_insar

(2)采用最小二乘估计得到趋势项的系数矩阵β：

β＝(F_insar′F_insar)^-1F_insarZ_insar

(3)计算去除趋势项的平稳残余项矢量：

(4)计算残余项的变差函数值：

根据InSAR坐标计算InSAR点位每两点的距离，然后对距离进行分组，当距离落入[kh±ε(h)](k为常数，ε(h)为容差)范围内就分在同一组，认为这些点的距离为kh，统计落入[kh±ε(h)]范围内的距离数目，记为N(h)，然后分别求出距离为kh实验变差函数值γ(h)，其值的计算公式如下：

其中s为空间位置，h为空间距离，以及分别为平稳残余项矢量中位置为s及s+h的元素，N(h)为空间位置相距h的观测值数。

(5)变差函数拟合

常用的变差函数模型主要有球状模型、指数模型、高斯模型以及幂函数模型等。而球状模型是最能代表空间变化的模型，其函数模型为：

其中C₀为块金值，C₀+C为基台值，a为变程。

将步骤(4)中计算的变差函数值γ(h)选择球状模型进行函数拟合求得变差函数参数C₀、C和a。

3)计算空间场矩阵H

(1)计算GNSS点位的空间协方差矩阵

依据步骤2)得到的空间协变差函数计算GNSS点位的空间协方差矩阵 中第ij行元素Σ_ij是通过下式计算的：

σ_ij＝C₀+C-γ(h)

其中C₀+C为拟合出来的变差函数的基台值，γ(h)为GNSS点位相距h的变差函数值。

(2)计算趋势矩阵A以及弯曲能量矩阵B

其中F_gnss为GNSS数据的趋势场，维数为m*q。其构造方法和步骤1)中的一致，且选择的趋势场模型是相同的。

(3)计算空间场矩阵H

进行降维处理，将B进行谱分解：

B＝UDU′ Bu_i＝d_iu_i

其中u_i是第i个特征向量，d_i是第i个特征值。m_gnss为GNSS观测值数，diag(·)为求对角阵算子。将m_gnss个特征值进行非降序排列：

而空间场的维数p是由下式决定的：

当e大于某一比例(一般为90％-95％)时即可得到p的取值。

空间场的计算如下式：

其中q为趋势场F的列数。

4)使用EM估计计算克里金Kalman滤波中的模型参数

EM估计是一种在数据缺失的情况下的极大似然估计。在Kalman滤波中完全数据为：{α₀,α_t,Z_t}＝{α₀,α₁,...,α_n,Z₁,...,Z_n}，其中缺失数据是不可观测的待求状态量α_t，已知数据是观测值Z_t，Z_t为t时刻的观测值。而待估计的模型参数为θ＝{μ，∑，Φ，Q，R}，其中μ和Σ为初始状态参数α₀的期望和协方差矩阵，Φ为状态转移矩阵，R和Q分别为观测噪声协方差矩阵与状态噪声协方差矩阵。完全数据的对数似然函数为：

因为缺失的数据α_t，所以并不能直接极大化完全数据的对数似然函数来得到模型参数θ。需要采用EM估计估计模型参数θ。EM估计分为E步和M步。E步是在参数θ_k下计算完全数据的对数似然函数的期望，M步是指极大化完全数据的对数似然函数的期望求得新的参数θ_k+1。反复迭代E步和M步直到对数似然函数稳定且不再增加为止。具体步骤如下：

假设第k次迭代计算的模型参数值为θ_k＝{μ_k,Σ_k,Φ_k,Q_k,R_k}，则在第k+1次迭代中：

E步：计算完整数据的对数似然函数期望

其中tr(·)为求迹算子，并且

而是在参数为θ_k通过克里金Kalman向前向后滤波得到的，具体滤波流程如下：

定义Kalman向前向后滤波估计如下：

其方差-协方差矩阵为：

Kalman向前滤波：

计算t时刻的状态预报值：

计算预报值的协方差矩阵：

计算t时刻的增益矩阵：

计算t时刻的状态估值：

计算状态估值的协方差矩阵:

其中t＝1,2,…,n，α₀的期望为μ，协方差矩阵为Σ。

Kalman向后滤波：

计算向后滤波的t-1时刻的增益矩阵：

计算向后滤波的t-1时刻的状态估值

计算向后滤波的t-1时刻状态估值的协方差矩阵：

其中t＝n,n-1,…,1。而t时刻与t-1时刻的状态估值之间的协方差矩阵可以由下式计算：

其中t＝n,n-1,…,2，并且当t＝n时

M步：通过最大化完全数据的对数似然函数的期望来求得新的模型参数值θ_k+1。对完全数据的对数似然函数的期望求导可得：

Φ＝BA^-1

Q＝n^-1(C-BA^-1B′)

此时求得k+1次迭代的模型参数的新值θ_k+1＝{μ_k+1,Σ_k+1,Φ_k+1,Q_k+1,R_k+1}。

完全数据的对数似然函数值可以由下式代替：

重复迭代E步与M步直到上式稳定不再增加为止。若稳定，则第k+1次迭代的新值θ_k+1就为EM估计所求得的模型参数最终估值。

5)克里金Kalman滤波

使用EM估计估计出来的模型参数进行克里金Kalman滤波，使用具有高时间分辨率的GNSS数据求得高时间分辨率的状态估值其主要是由5个迭代公式组成的：

t时刻状态量的预报值：

预报值的协方差矩阵：

增益矩阵：

t时刻状态量的估值：

估值的协方差矩阵：

6)时空插值预测

为得到高时空分辨率的成果，需要在时空上进行预测。

设Z_t(s)表示在t时刻且位置为s未观测点的观测值，H^*表示未观测点的空间场，其构建方法和构造已知点的方法一致，v^*(t)表示Z_t(s)的噪声。Z_t为t时刻已知点的观测值，v(t)为量测噪声，Σ_v为其协方差矩阵。令Σ_V、和使用2)步骤中拟合的协变差函数来描述。则克里金预报值为：

增加s的空间分辨率，所得到的Z_t(s)就为高时空分辨率的成果。

预报值精度可由克里金方差表示：

实例分析：

图2(a)、图2(b)反映了GNSS数据和InSAR数据计算的空间实验变差值与距离分布情况。

由图2(a)图可知，GNSS监测点在空间上分布稀疏，根本不能体现形变数据在空间上的平稳性。所计算出来的实验变差函数值杂乱无章，这样变差函数值强制拟合出来的变差函数明显是不合理的。

由图2(b)图可知，InSAR点位在空间上分布密集，可以较好的体现出空间上的平稳性。所计算出来的实验变差函数值在小于变程时逐步升高，在大于变程时逐渐趋于稳定于基台值。这样拟合出来的空间变差函数具有一定的可靠性。

Claims (6)

1.一种融合GNSS和InSAR数据获取高时空分辨率形变序列的方法，其特征在于，包括以下步骤：

1)计算InSAR数据的趋势场F_insar，并去除InSAR数据的趋势项；

2)采用曲面拟合的方式，利用所述趋势场F_insar构造曲面拟合趋势项的误差方程，使用最小二乘求得趋势项曲面的系数矩阵，计算去除趋势项的平稳残余项，选择变差函数模型，并利用所述平稳残余项拟合变差函数γ(h)；

3)计算GNSS数据的空间趋势场矩阵F_gnss，利用上述变差函数γ(h)计算GNSS数据点位之间的空间协方差矩阵

4)根据GNSS数据的空间趋势场矩阵F_gnss以及空间协方差矩阵计算弯曲能量矩阵B；

5)利用所述弯曲能量矩阵B以及趋势场矩阵F_gnss构建克里金Kalman滤波的空间场H；

6)使用EM估计计算克里金Kalman滤波中的模型参数，进行克里金Kalman滤波；

7)进行时空插值，得到高时空分辨率的形变序列。

2.根据权利要求1所述的融合GNSS和InSAR数据获取高时空分辨率形变序列的方法，其特征在于，所述步骤1)中，InSAR数据的趋势场F_insar计算方法为：若趋势场为常数场，则趋势场F_insar为值全为1且维数为m_insar×1的矢量，即：F_insar＝[1,1,...,1]′；若趋势场为线性场，则趋势场F_insar为维数为m_insar×3的矩阵，即：F_insar＝[x′；y′；1]′；若趋势场为二次曲面场，则趋势场F_insar为m_insar×6的矩阵，即：F_insar＝[(x²)′；(y²)′；(xy)′；x′；y′；1]′；其中，m_insar为InSAR数据观测值的个数，InSAR数据观测值坐标x,y均为维数为m_insar×1的矢量。

3.根据权利要求2所述的融合GNSS和InSAR数据获取高时空分辨率形变序列的方法，其特征在于，所述步骤2)中，去除趋势项的平稳残余项矢量其中，Z_insar为观测值；趋势项曲面的系数矩阵β＝(F_insar′F_insar)^-1F_insarZ_insar。

4.根据权利要求3所述的融合GNSS和InSAR数据获取高时空分辨率形变序列的方法，其特征在于，所述步骤3)中，空间协方差矩阵中第ij行元素Σ_ij通过下式计算：Σ_ij＝C₀+C-γ(h)；其中，γ(h)为变差函数，而s为空间位置，h为空间距离，以及分别为平稳残余项矢量中位置为s及s+h的元素，N(h)为以h为中心且给定的容差ε(h)为半径的观测值数，C₀+C为变差函数的基台值。

5.根据权利要求4所述的融合GNSS和InSAR数据获取高时空分辨率形变序列的方法，其特征在于，所述步骤4)中，弯曲能量矩阵其中，

6.根据权利要求5所述的融合GNSS和InSAR数据获取高时空分辨率形变序列的方法，其特征在于，空间场其中，d_q+1、d_q+2、…d_p为弯曲能量矩阵B进行谱分解后的特征值；u_q+1、u_q+1、…u_p为弯曲能量矩阵B进行谱分解后的特征向量；p为空间场的维数。