基于站间相关性空间结构函数构建的GNSS成像方法

摘要

本发明提供一种基于站间相关性空间结构函数构建的GNSS成像方法，包括输入GNSS测站坐标时间序列观测值，以及每个GNSS测站的速度及不确定性；结合研究区域的地质学与大地测量学结果，对GNSS网内的测站进行聚类划分，得到聚类区域；计算各聚类区域内任意两测站所构成测站对之间的相关系数，根据测站间相关系数进行数据池划分，获取每个聚类区域中若干数据池及在各数据池内的GNSS测站对；在每个聚类区域内，计算每个数据池中所有测站对相关系数的中值和绝对中位差，构建各聚类区域的空间结构函数，标准化形成整个GNSS网最终的空间结构函数；依据速度不确定性和空间结构函数，确定研究范围内所有测站的权，利用空间插值方法进行空间插值，形成影像。

说明书

技术领域

本发明属于GNSS数据精密处理技术领域，涉及一种基于GNSS多测站时间序列构建空间结构函数的GNSS成像方法。

背景技术

近年来，GPS/GNSS监测网络的投入运行，在覆盖范围上获得了巨大的扩展，测站数目日益增多，包括中国地壳运动观测网络和中国大陆构造环境监测网络，美国的PBO网络，欧洲的EPN，日本的GEONET等网络，产生了大量的呈级数增长的观测数据。

利用这些密集的GNSS监测网络，可以合理利用GNSS时间序列确定测站速度及其不确定性，从而进行地壳形变成像；然而，利用GNSS成像技术，遇到的瓶颈问题包括如何构建反映地壳形变空间特征的空间结构函数。

发明内容

本发明提供一种GNSS成像方法，实现基于各向同性假设的空间结构函数的构建，该方法利用研究区域内GNSS测站速度及不确定性所反映的统计特征进行构建，有效提高GNSS成像的空间分辨率。

为解决上述技术问题，本发明采用的技术方案为一种基于站间相关性空间结构函数构建的GNSS成像方法，包括以下步骤，

步骤1，输入GNSS测站坐标时间序列观测值，以及每个GNSS测站的速度v_n及其不确定性σ_n，n＝1,2,…,N；其中，N为整个GNSS网中GNSS测站的总数，n为整个GNSS网中GNSS测站的标号；

步骤2，结合研究区域的地质学与大地测量学结果，对GNSS网内的测站进行聚类划分，设得到K个聚类区域，第k个聚类区域记作I_k，k＝1,2,...,K；

步骤3，计算各聚类区域内任意两测站所构成测站对之间的相关系数，设第k个聚类区域I_k内有I个测站，计算第k个聚类区域I_k内的GNSS网中对应的任意两测站i,j所构成测站对之间的相关系数γ_i,j；其中，k＝1,2,...,K，i＝1,2,...,I；j＝1,2,...,I；i≠j；

步骤4，在每个聚类区域I_k内，根据测站间相关系数γ_i,j进行数据池划分，获取每个聚类区域中若干数据池及在各数据池内的GNSS测站对，设聚类区域I_k有Q个数据池，某个数据池内测站数目记作m，该数据池内可组成的GNSS测站对为m(m-1)/2；

步骤5，在每个聚类区域I_k内，计算每个数据池内任意测站对之间的速度差值以及数据池中所有测站对相关系数的中值，实现方式如下：

依据步骤1中获得的速度，计算每个数据池内任意两测站a,b所构成测站对之间的速度差值dv_a,b以及数据池中所有测站对相关系数的中值M(γ_q)，q＝1,...,Q；

步骤6，在每个聚类区域I_k内，对每个数据池将各测站对的速度差值dv_a,b按从小到大排序，求取测站对的速度差值序列的中位数M(dv)；

步骤7，在每个聚类区域I_k内，对每个数据池分别计算数据池中每个测站相对于M(dv)的相对量，按从小到大排序，求取中位数，记为绝对中位差M(δv)；

步骤8，构建各聚类区域I_k的空间结构函数ssf_k如下，

其中，M(γ_q)表示聚类区域I_k内第q个数据池中的相关系数中位数，M(δv)_q-1表示聚类区域I_k内第q-1个数据池的绝对中位差，M(δv)_q表示聚类区域I_k内第q个数据池的绝对中位差；

步骤9，对各聚类区域I_k的空间结构函数ssf_k标准化，形成整个GNSS网最终的空间结构函数ssf(γ)；

步骤10，依据步骤1中获得的速度不确定性σ_n和步骤9中确定的空间结构函数ssf(γ)值，确定研究范围内所有测站的权利用空间插值方法进行空间插值，形成影像。

而且，步骤4中，在每个聚类区域I_k内，根据测站间相关系数γ_i,j进行数据池划分，采用以下方式实现，

根据聚类区域内所有的相关系数，获取相关系数的最大值γ_max和最小值γ_min；

输入聚类区域I_k内测站对之间的相关系数，并且按照相关系数从小到大的顺序排列；

根据相关系数的大小，取作为数据池的大小，如果某个测站对的相关系数cor大小满足则把该测站对放入第num个数据池。

和现有技术相比，本发明具有如下特点和有益效果：

本发明提供了一种基于GNSS时间序列进行成像中的关键问题的解决方案，即提出一种基于测站对相关系数的空间结构函数的构建方案，实现GNSS成像。本发明针对GNSS成像中的关键问题——如何构建空间结构函数，提出创新性的解决方法，得到令人满意的解决效果：利用GNSS时间序列估计的测站速率具有较高的时间分辨率和可靠性，有效克服了常规卫星影像成像方法在时间分辨率上的不足，测站对之间的相关性较好地描述站网内多个测站空间关系和数据统计特征，成为构建密集的GNSS站网间的空间结构的关键。基于站间相关性的空间结构函数构建，对异常点进行降权处理，充分利用空间可用点，从而实现GNSS成像，提高了对地壳形变空间模式的认知。

附图说明

图1为本发明实施例的流程示意图。

具体实施方式

为使本发明目的、技术方案及有益效果更加清楚明白，下面将结合附图及实施例，进一步说明本发明。

参见图1，本发明实施例提供的基于站间相关性空间结构函数构建的GNSS成像方法，包括步骤：

步骤1，输入GNSS测站坐标时间序列观测值，以及每个GNSS测站的速度v_n及其不确定性σ_n(n＝1，2，…，N)；其中，N为整个GNSS网中GNSS测站的总数，n为整个GNSS网中GNSS测站的标号；

步骤2，结合研究区域的地质学与大地测量学结果(如地质断层位置信息、基于长期GPS观测的活动块体划分模型等)，对GNSS网内的测站进行聚类划分，设得到K个聚类区域，第k个聚类区域记作I_k,(k＝1,2,...,K)；

地质学通常会给出断层和活动块体的空间分布，大地测量学结果会提供测站的位置信息；根据大地测量的位置信息，判断对应的测站点落在哪个活动块体内，或者落在哪几个断层之间的空间区域；例如，落在某个活动块体围限区域内的所有测站点属于一类，从而实现初步的聚类划分。

步骤3，计算各聚类区域内任意两测站所构成测站对之间的相关系数，设第k个聚类区域I_k内有I个测站，计算第k个聚类区域I_k内的GNSS网中对应的任意两测站i,j所构成测站对之间的相关系数γ_i,j；其中，k＝1,2,...,K，i＝1,2,...,I；j＝1,2,...,I；i≠j；

相关系数的计算可采用现有技术，本发明不予赘述；

步骤4，在每个聚类区域I_k,(k＝1,2,...,K)内，根据测站间相关系数γ_i,j进行数据池划分，获取每个聚类区域中若干数据池及在各数据池内的GNSS测站对，每个数据池内测站数目不等，设聚类区域I_k有Q个数据池，设某个数据池内测站数目记作m，该数据池内可组成的GNSS测站对为m(m-1)/2；

依据相关系数对每个聚类区域进行测站对分类，并为后续速度差计算提供前提，是本发明的关键创新之一。

当两测站重合时，两者之间的相关系数为1；而当两站之间相距较远时，两站之间的相关系数减小；一般而言，任意一个测站对(两个测站)之间的相关系数取值范围为[-1,1]。需要指出一点，随着研究区域范围不同，GNSS测站网的疏密程度不同，相关系数的取值范围也存在较大差异。因此，本发明进一步提出依据相关系数γ_l,j数据池划分的具体实现方式为，步骤3已经计算聚类区域中的任意两个测站对应的相关系数γ_l,j，根据聚类区域内所有的相关系数，获取相关系数的最大值γ_max和最小值γ_min。依据步骤2中的初步聚类划分结果，在每个聚类区域内开展以下工作：

(1)根据步骤3所得相关系数，输入聚类区域I_k内测站对之间的相关系数，并且按照相关系数从小到大的顺序排列；

(2)根据相关系数的大小，取作为数据池的大小，将聚类区域I_k内每个测站对放入对应的数据池，保证每个数据池内有足够的数据，一般不得少于5个数据对；如果数据对少于5个，将相邻的数据池进行合并；

具体实施时，将聚类区域I_k内测站对放入对应的数据池，可采用的方式为：

如果某个测站对的相关系数cor大小满足则把该测站对放入第num个数据池。具体实施时，数据池的数目Num可根据需要设置，num的取值为1,2,…Num。

步骤5，在每个聚类区域I_k,(k＝1,2,...,K)内，计算计算每个数据池内任意测站对之间的速度差值以及数据池中所有测站对相关系数的中值，实现方式如下：

依据步骤1中获得的速度，计算每个数据池内任意两测站a,b所构成测站对之间的速度差值dv_a,b以及数据池中所有测站对相关系数的中值，设第q个数据池中所有测站对相关系数的中值M(γ_q)(q＝1,...,Q)，其中，a＝1,2,...,m；b＝1,2,...,m；a≠b：

dv_a,b＝|v_a-v_b|；>

其中，v_a是数据池内第a个GNSS测站的速度，v_b是数据池内第b个GNSS测站的速度。

步骤6，在每个聚类区域I_k,(k＝1,2,...,K)内，对每个数据池求取测站对的速度差值序列的中位数，实现如下：

将某一个数据池中测站对的速度差值dv_a,b按从小到大排序，求取速度差值序列的中位数M(dv)：

其中，mod()表示整除，上式提供了当p为奇数和偶数时的不同取值，p表示数据池内的GNSS测站对的数目，p＝m(m-1)/2；dv_(p+1)/2表示p为奇数时，速度差值序列的中位数取第(p+1)/2测站对的速度差值，公式(2)中第二行表示当p为偶数时，速度差值序列的中位数为序列中段的两个速度差值的数值平均，dv_p/2表示第p/2测站对的速度差值，dv_(1+p/2)表示第(1+p/2)测站对的速度差值。

步骤7，在每个聚类区域I_k,(k＝1,2,...,K)内，对每个数据池计算数据池中每个测站相对于速度差值中位数的相对量，按从小到大排序，求取中位数，记为绝对中位差M(δv)；实现方式如下，

对于数据池中某个测站a，相对于速度差值中位数的相对量记为δv_a：

δv_a＝v_a-M(dv)，>

其中，v_a是数据池内第a个GNSS测站的速度，a＝1,2,...,m；

并将速度差值中位数的相对量δv_a按从小到大排序，按照步骤6的方法，求取中位数M(δv)；

步骤8，构建各聚类区域I_k的空间结构函数ssf_k，实现如下，

比较某个聚类区域I_k内Q个数据池中所获得的绝对中位差M(δv)的大小，构建整个聚类区域I_k的空间结构函数，

其中，M(γ_q)(q＝1,...,Q)表示聚类区域I_k内第q个数据池中的相关系数中位数，M(δv)_q-1表示聚类区域I_k内第q-1个数据池的绝对中位差，M(δv)_q表示聚类区域I_k内第q个数据池的绝对中位差，通过比较相邻数据池的绝对中位差，可以快速找到所有数据池的绝对中位差中最大者，具体实施时每个具体区域的数据池可能达到数万个，这样在数据池数量较大时可以显著提高效率；

步骤9，对各聚类区域I_k的空间结构函数ssf_k标准化，形成整个GNSS网最终的空间结构函数ssf(γ)；

标准化过程的具体实现为：

(1)令整个GNSS网内最大的空间结构函数取值为1，记作max(ssf_k)(k＝1,...,K)；

(2)每个聚类区域内的空间结构函数相对于最大空间结构函数的归一化结果可以表示为：

ssf_{k_new}＝ssf_k/max(ssf_t)

其中，ssf_k表示第k个聚类区域的空间结构函数,ssf_{k_new}表示第k个聚类区域归一化后的空间结构函数。

步骤10，依据步骤1中获得的速度不确定性σ_n和步骤9中确定的空间结构函数ssf(γ)值，确定研究范围内所有测站的权w_n，利用常规的空间插值方法进行空间插值，形成影像：

具体实施时，以上流程可采用计算机软件方式实现自动运行流程。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。