[发明专利]一种顾及多路径和同态误差的GNSS地震地表位移监测方法

权利要求书

1.一种顾及多路径和同态误差的GNSS地震地表位移监测方法，其特征在于，包括以下步骤：

步骤1，对测区内N个GNSS测站逐站进行PPP处理，获得N个GNSS测站的东向位移序列、北向位移序列、垂向位移序列，通过构造平均能量函数计算GNSS测站的位移序列中所携带的地震波信号的平均能量，进一步探测每个GNSS测站地震波的到达时刻并探测每个GNSS测站地震波的结束时刻；

步骤2，结合GNSS测站地震波的到达时刻计算GNSS测站的地震波到达时刻参考值，结合GNSS测站地震波的结束时刻计算GNSS测站的地震波结束时刻参考值，扣除震后GNSS测站东向、北向、垂向位移序列中的同震形变；

步骤3，利用主成分分析法分别计算非震时期GNSS测站东向、北向、垂向位移序列的空间响应和主成分；

步骤4，结合GNSS测站的地震波到达时刻参考值与GNSS测站的地震波结束时刻参考值，确定出多个不受地震波影响的GNSS测站，用不受地震波影响的GNSS测站位移分别计算出受地震波影响时段内测站东向、北向和垂向位移序列的主成分；

步骤5，结合步骤3计算的非震时期内GNSS测站东向、北向和垂向位移序列的空间响应，以及步骤4所述受地震波影响时段内GNSS测站东向、北向和垂向位移序列的主成分，计算出受地震波影响时段包括受震测站在内的GNSS测站在东向、北向和垂向位移序列中的同态误差；结合步骤3计算的非震时期内GNSS测站东向、北向和垂向位移序列的空间响应和主成分，计算非震时期内GNSS测站在东向、北向和垂向位移序列中的同态误差；

步骤6，从GNSS测站的东向、北向和垂向位移序列中扣除步骤5所得非震时期内的同态误差和受地震波影响时段内GNSS测站的同态误差，再进行恒星日滤波消除多路径误差，最终获得东向、北向和垂向高精度的GNSS地表位移序列；

步骤3中所述利用主成分分析法分别计算非震时期GNSS测站东向、北向、垂向位移序列的空间响应和主成分为：

所述计算非震时期GNSS测站第i个方向的位移分量矩阵可表示为：

其中，表示震前第i个方向GNSS测站的原始位移矩阵，表示震后第i个方向GNSS测站扣除同震形变后的位移矩阵，C_pre为震前观测时段内的观测历元总数，C_after为震后观测时段内的观测历元总数，C_pa＝C_pre+C_after为非震时期内的观测历元总数，N为GNSS测站数，矩阵的每一列代表各测站第i个方向的位移向量，i＝1表示东向，i＝2表示北向，i＝3表示垂向；

第i个方向位移分量构成的协方差矩阵可表示为其又可分解为：

其中，为特征值构成的矩阵，特征值从大到小顺序排列，表示非震期间内第i个方向特征值对应的特征向量构成的矩阵，特征向量之间相互正交，由于任意n维矩阵，均可由n组正交基表示，n∈[1,N]，N为GNSS测站数，i＝1表示东向，i＝2表示北向，i＝3表示垂向；

选取中的前n组基来重新表示第k个GNSS测站第i个方向在第t时刻的位移分量：

其中，表示非震期间内第i个方向在第t时刻构成的位移矩阵，表示非震期间内第i个方向第j个主成分在第t时刻的值，表示非震期间内第i个方向特征向量矩阵中第j个特征向量中的第k个值，t∈T_pa，j∈[1,n]，k∈[1,N]，i＝1表示东向，i＝2表示北向，i＝3表示垂向；

可表示为：

其中，表示非震期间内第k个测站第i个方向在t时刻的位移,表示非震期间内第i个方向第j个主成分，又称非震期间内第i个方向第j个主成分对应的空间响应，t_tep∈T_pa，j∈[1,n]，i＝1表示东向，i＝2表示北向，i＝3表示垂向；

第j个模态可由第j个主成分和第j个空间响应表示；

前n个模态代表区域中的同态误差；

步骤4所述不受地震波影响的GNSS测站为：

结合步骤1所述，认为在任意时刻GNSS测站位移序列中所携带的地震波信号的平均能量F_k均小于3倍F₀时，该GNSS测站既为不受地震波影响的测站；

步骤4所述用不受地震波影响的GNSS测站位移分别计算出受地震波影响时段内测站东向、北向和垂向位移序列的主成分：

用Ns个受地震波影响时段内不受地震波影响的GNSS测站位移计算出受地震波影响时段内GNSS测站第i个方向第j个主成分在第t_nom时刻的值：

其中，表示受地震波影响时段内第i个方向第j个主成分在第t_nom时刻的值，表示受地震波影响时段内第i个方向第j个主成分，为受地震波影响时段内第k个不受地震波影响的测站在第i个方向第t_nom时刻的位移，为受地震波影响时段内第i个方向第j个模态对应的空间响应中的第k个值，t_nom∈T_move，j∈[1,Ns]，k∈[1,Ns]，Ns为受地震波影响时段内不受地震波影响的测站个数，i＝1表示东向，i＝2表示北向，i＝3表示垂向；

结合步骤2所述，将地震波到达时刻参考值和地震波结束时刻参考值之间的时段定义为受地震波影响的时段:

T_move＝[t_arrival,t_stop]

步骤5所述的结合步骤3计算的非震时期内GNSS测站东向、北向和垂向位移序列的空间响应，以及步骤4所述受地震波影响时段内GNSS测站东向、北向和垂向位移序列的主成分，计算出受地震波影响时段包括受震测站在内的GNSS测站在东向、北向和垂向位移序列中的同态误差：

首先，通过判断标准化空间响应和特征值占比，确定参与同态误差计算的模态：

对非震时期内第i个方向第j个模态对应的空间响应为：

进行标准化为：

其中，表示为非震时期内第i个方向第j个模态对应的空间响应中的第s个值，表示非震时期内第i个方向第j个模态对应的空间响应中所有元素的最大值，s∈[1,N]，j∈[1,N]，N为GNSS测站数，i＝1表示东向，i＝2表示北向，i＝3表示垂向；

计算第i个方向第j个空间响应对应的统计量为：

其中，表示非震时期内第i个方向第j个标准化空间响应向量中大于阈值thre_sr的测站数，j∈[1,N]，N为GNSS测站数，i＝1表示东向，i＝2表示北向，i＝3表示垂向；

计算与非震时期内第i个方向第j个空间响应向量相关的特征值所对应的统计量为：

若统计量满足以下条件时，则认为第i个方向第j个模态可用于计算第i个方向位移的同态误差：

其中，thre1表示空间响应统计量对应的阈值，thre1∈[50％,100％]，thre2表示特征值统计量对应的阈值thre2∈[1％,100％]；

第j个模态可由第j个主成分和第j个空间响应表示；

前n个模态代表区域中的同态误差；

根据判断条件筛选出n0个模态，计算出受地震波影响时段内包含受震测站在内的所有测站的同态误差：

其中，表示受地震波影响时段内包括受震测站在内的第k个测站在第i个方向上位移序列中的同态误差，为受地震波影响时段内第i个方向第j个模态对应的主成分，为步骤2计算的第i个方向上第j个模态对应的空间响应中的第k个值，j∈[1,n0]，k∈[1,N]，N为GNSS测站数，i＝1表示东向，i＝2表示北向，i＝3表示垂向；

步骤5所述的结合步骤3计算的非震时期内GNSS测站东向、北向和垂向位移序列的空间响应和主成分，计算非震时期内GNSS测站东向、北向和垂向位移序列的同态误差：

根据筛选出的n0个模态，计算出非震时期内GNSS测站的同态误差：

其中，表示非震时期内第k个测站在第i个方向上位移序列中的同态误差，为步骤3计算的非震时期内第i个方向第j个模态对应的主成分，为步骤3计算的非震时期内第i个方向上第j个模态对应的空间响应中的第k个值j∈[1,n0]，k∈[1,N]，N为GNSS测站数，i＝1表示东向，i＝2表示北向，i＝3表示垂向；

步骤6从N个测站位移中扣除步骤5所得受地震波影响时段内包括受震测站在内的GNSS测站东向、北向和垂向位移序列的同态误差，以及非震时期内GNSS测站东向、北向和垂向位移序列的同态误差：

其中，表示所有观测时间内第k个测站在第i个方向上位移序列中的同态误差，表示受地震波影响时段内包括受震测站在内的第k个测站在第i个方向上位移序列中的同态误差，表示非震时期内第k个测站在第i个方向上位移序列中的同态误差，表示第k个测站在第i个方向上的原始位移序列，k∈[1,N]，N为GNSS测站数，i＝1表示东向，i＝2表示北向，i＝3表示垂向；

步骤6所述，将扣除同态误差后的位移序列进行恒星日滤波，最终获得高精度位置序列：

其中，表示第k个GNSS测站在第i个方向上经过恒星日滤波计算得到的多路径误差：

其中，表示第d天第k个测站在第i个方向上扣除同态误差后的位移序列，k∈[1,N]，d∈[1,D]，N为GNSS测站数，D表示参与计算多路径误差的天数，i＝1表示东向，i＝2表示北向，i＝3表示垂向。