热膨胀效应对GPS坐标时间序列影响的精确量化方法

摘要

本发明公开了一种热膨胀效应对GPS坐标时间序列影响的精确量化方法，包括：收集IGS基准站的地表温度时间序列，采集测站信息；计算热膨胀效应对IGS基准站坐标时间序列的总影响，即序列A；采用地表质量负载重分布模型计算环境负载引起的IGS基准站的垂直方向位移时间序列，即序列B；采用GPS精密数据处理模型得到IGS基准站的坐标时间序列，即序列C；从序列C中扣除序列B，得序列C‑B；从序列C‑B中扣除序列A，得序列C‑B‑A；基于序列C‑B和序列C‑B‑A提出数值指标量化热膨胀效应对GPS坐标时间序列的影响。本发明可获得温度变化引起的IGS基准站垂直方向实际运动，并采用数值指标准确量化热膨胀效应的影响。

说明书

技术领域

本发明属于非线性时间序列分析技术领域，并具体涉及一种热膨胀效应对GPS坐标时间序列影响的精确量化方法。

背景技术

GNSS坐标时间序列的变化受多种因素的影响，除了固体潮、海洋潮汐等已模型化的因素以外，还会受到环境负载(非潮汐海洋负载、大气压负载、水文负载等)及热膨胀效应等其他地球物理效应的影响。然而，目前已发布的坐标参考框架(如ITRF2008)提供的速度场模型未考虑这部分非线性运动，不能精确描述基准站的实际变化特征。因此，精化热膨胀影响计算模型，定量确定热膨胀造成的测站坐标时间序列季节性信号的影响，并探究相关的地球物理过程及其物理机制，不仅能得到更为精确、真实的测站位置和速度信息，反映出测站真实的运动变化，而且有助于研究冰后回弹及海平面变化、反演冰雪质量变迁等地球动力学过程，具有非常重要的理论意义和应用价值。

当前的研究成果表明，GNSS测站附近温度变化引起的测站天线观测墩和基岩热膨胀效应是坐标时间序列季节性信号的潜在贡献源之一，尤其是垂直方向上，在中高纬度地区测站的周年影响振幅可达3～4mm(Yan等，Contributions of thermal expansion of monuments and nearby bedrock to observed GPS height changes，Geophysical Research Letters，2009)。但是受限于全球IGS基准站测站观测墩和基岩信息缺失，难以精确计算并量化全球IGS基准站热膨胀效应在垂直方向上的影响及其对季节性信号的贡献。此外，使用不同的数学模型和评价指标也会造成不同学者计算出的结果差异较大。

发明内容

针对目前全球GNSS基准站垂直方向热膨胀效应影响难以精确量化的问题，本发明提供了一种热膨胀效应对GPS坐标时间序列影响的精确量化方法，该方法可精确计算由热膨胀引起的GNSS基准站垂直方向的位移时间序列，即垂直方向热膨胀影响时间序列。

为解决上述技术问题，本发明采用如下技术方案：

热膨胀效应对GPS坐标时间序列影响的精确量化方法，包括：

S1收集待分析区域IGS基准站的地表温度时间序列，采集测站信息，所述的测站信息包括观测墩高度、观测墩热膨胀系数、基岩泊松比、基岩热膨胀系数和基岩热扩散系数；

S2利用热膨胀效应影响模型，分别计算各IGS基准站由基岩热膨胀引起的垂直方向位移Δh(t)和由观测墩热膨胀引起的垂直方向位移ΔL(t)，Δh(t)和ΔL(t)之和即热膨胀效应的总影响，记为序列A，其中，t表示时刻；

S3采用地表质量负载重分布模型计算环境负载引起的IGS基准站的垂直方向位移时间序列，记为序列B；

S4采用GPS精密数据处理模型得到IGS基准站的原始坐标时间序列，记为序列C；

S5从序列C中扣除序列B，得序列C-B；从序列C-B中扣除序列A，得序列C-B-A；

S6采用序列C-B和序列C-B-A间的相关系数、加权均方根误差的改变程度和周期振幅的变化程度三个指标，来量化热膨胀效应对GPS坐标时间序列的影响。

所述的地表温度时间序列来源于由NCEP/NCAR再分析数据集提供的全球地表温度格网或EMCWF提供的高分辨率全球/区域地表温度格网。

步骤S2进一步包括：

采用热膨胀效应影响模型计算IGS基准站由基岩热膨胀引起的垂直方向位移Δh(t)；

采用热膨胀效应影响模型计算IGS基准站由观测墩热膨胀引起的垂直方向位移ΔL(t)；

Δh(t)和ΔL(t)之和即热膨胀效应的总影响；

上述，μ为基岩泊松比；α₁为基岩热膨胀系数；κ为基岩热扩散系数；N表示谐波数；a_i、b_i、c_i分别表示地表温度时间序列第i个谐波的振幅、角频率和相位，可通过最小二乘拟合法对地表温度时间序列估计获得；α₂为观测墩热膨胀系数；H为观测墩高度；T(t)为t时刻观测墩附近的地表温度，为观测墩附近的地表温度年平均值，T(t)和值从地表温度时间序列获得。

步骤S3中，环境负载包括大气负载、非潮汐海洋负载和大陆水储量负载。

步骤S6中，序列C-B和序列C-B-A间的加权均方根误差的改变程度PWRMS为：

其中，WRMS(U)和WRMS(U-H)分别表示序列U和U-H的加权均方根误差，序列U和H分别表示序列C-B和序列C-B-A。

步骤S6中，序列C-B和序列C-B-A间的周期振幅的变化程度PE为：

其中，k为周期数；PE(k)表示第k个周期下的周期振幅变化程度；A_a(k)和A_b(k)分别表示序列C-B和序列C-B-A的第k个周期的振幅；

序列C-B和序列C-B-A的第k个周期的振幅采用如下方法获得：

采用频谱分析法确定序列C-B和序列C-B-A的周期，采用最小二乘拟合法估计各周期的振幅。

本发明采用高空间分辨率的地表温度数据并顾及观测墩高度等测站和基岩信息，计算热膨胀效应引起的垂直方向位移时间序列，并通过计算序列间相关系数、扣除热膨胀影响前后均方根误差是否改变及改变大小、以及扣除热膨胀影响前后序列显著周期振幅的变化程度，量化并评价热膨胀对GPS坐标时间序列的影响和贡献。

与现有技术相比，本发明具有如下特点：

(1)顾及不同测站的观测墩高度、泊松比等信息，精化热膨胀效应影响计算模型，得到的时间序列更能够表征由温度变化引起的IGS基准站垂直方向实际运动；

(2)采用相关系数、序列加权均方根值(WRMS)和季节性信号贡献量等定量指标量化热膨胀对GPS坐标时间序列的影响和贡献，较已有的评价标准更为准确，有助于进一步揭示GPS坐标时间序列季节性信号的物理机制，为建立估计基准站非线性变化的地球参考框架提供支持与借鉴。

附图说明

图1为实施例计算得到的各IGS基准站由热膨胀引起的垂直方向位移周年振幅；

图2为实施例计算得到的各IGS基准站由热膨胀引起的垂直方向位移半周年振幅；

图3为计算得到的全球各IGS基准站热膨胀效应影响评价指标结果，其中，图(a)、图(b)、图(c)分别为相关系数、PWRMS值和PE值的评价结果；

图4为扣除热膨胀效应和环境负载影响前后的垂直方向坐标时间序列堆栈功率谱密度；

图5为加入环境负载和热膨胀效应改正后部分测站周年振幅；

图6为部分测站扣除环境负载和热膨胀效应影响前后的垂直方向坐标时间序列；

图7为本发明方法的具体流程示意图。

具体实施方式

下面结合附图和实施例对本发明技术方案做进一步的叙述。

本实施例是计算560个全球分布的IGS(International GNSS Service)基准站的热膨胀效应影响，参见图7，具体步骤如下：

步骤1，获取由美国国家环境预报中心(National Centers for Environmental Prediction，NCEP)和美国国家大气研究中心(National Center for Atmospheric Research，NCAR)联合推出的NCEP/NCAR再分析数据集，空间分辨率为2.5°×2.5°，并利用双线性内插法得到560个IGS基准站从1998年1月1日到2013年12月31日的地表温度时间序列。

步骤2，采集各IGS基准站的测站信息，所述的测站信息包括观测墩高度、观测墩热膨胀系数、基岩泊松比、基岩热膨胀系数和基岩热扩散系数。基于测站信息，利用热膨胀效应影响模型计算560个IGS基准站由热膨胀引起的垂直方向位移时间序列，记为序列A。图1和图2分别为采用公式(9)计算的IGS基准站由热膨胀引起的垂直方向位移周年和半周年振幅。

下面将对垂直方向位移时间序列的计算过程进行详细说明。

热膨胀效应引起的垂直方向位移分为地上和地下两个部分，即基岩热膨胀和观测墩热膨胀引起的垂直方向位移，两部分垂直方向位移相加即热膨胀效应的总影响。

根据弹性半空间模型的假设，在半无限空间体内有热传导方程：

式(1)中，t表示时刻；T表示温度；κ表示热扩散系数，d为基岩到地壳表面的垂直距离，d≥0，地壳表面为0；本发明中，κ为基岩的热扩散系数。

同时考虑到边际条件：

式(2)中，T′表示年平均温度；N表示谐波数；a_i、b_i、c_i分别表示地表温度时间序列第i个谐波的振幅、角频率和相位，可通过最小二乘拟合法对地表温度时间序列估计获得。

通过求解微分方程得到不同时刻t下不同深度d处基岩的温度解析解T(d,t)：

考虑广义胡克定律，作用在单元体上的热应力产生的线变形为：

式(4)中，E表示杨氏模量；μ表示弹性体的泊松比，本发明中μ表示基岩的泊松比；α表示线性膨胀系数(/℃)，ΔT表示温度变化量，σ_x、σ_y、σ_z分别表示x、y、z方向的应力元，ε_x、ε_y、ε_z分别表示x、y、z方向的单元体线形变。

假设平均温度对应热应力为0的状态，同时半无限空间体为统一的弹性介质，那么可以定义水平方向上有：

ε_x＝0,ε_y＝0,σ_z＝0>

得到单位体的垂直方向线形变ε_z为：

式(6)左右两边同时对深度d由0至∞积分，得到基岩热膨胀造成的垂直方向位移Δh(t)：

式(7)中，α₁表示基岩的线性热膨胀系数。

天线观测墩热膨胀造成的垂直方向位移ΔL(t)采用简单的线性膨胀模型近似描述：

式(8)中，α₂表示天线观测墩的线性热膨胀系数，与观测墩材质有关；H表示观测墩高度，包括观测墩高出地表和地下延伸到基岩的部分；T(t)表示t时刻观测墩附近的地表温度；为观测墩附近的地表温度年平均值。T(t)和值从地表温度时间序列获得。

本实施例中，将T(t)近似为观测墩的内部温度。Δh(t)和ΔL(t)之和即测站热膨胀效应的总影响，也即热膨胀引起的垂直方向位移时间序列。

步骤3，从全球地球物理流体中心GGFC(Global Geophysical Fluid Center)获取大气负载、非潮汐海洋负载以及大陆水储量负载等环境负载引起的IGS基准站垂直方向位移时间序列，记为序列B。

下面将本步骤进行详细说明：

环境负载造成的三维地表位移通过将描述地表质量分布的格网模型与描述地球的单位脉冲响应的格林函数求卷积获得，表现为负载及响应位置的函数。目前，GGFC提供了ATML、NTOL、CWSL三种负载变化造成的地表位移全球格网，用户通过内插格网得到任意指定经纬度处的负载位移时间序列。

其中，ATML根据NCEP提供的地表气压数据获得，空间分辨率为2.5°×2.5°，时间分辨率为6小时，具体为将残差气压格网与Farrell格林函数求卷积计算得到。NTOL位移时间序列采用ECCO模型提供的经卡尔曼滤波得到的海底压力(OBP)产品计算，空间分辨率同样为2.5°×2.5°，时间分辨率为0.5天。CWSL造成的地表位移采用GLDAS模型，包括月均雪水当量及土壤湿度值计算得到，空间分辨率为2.5°×2.5°，时间分辨率为30天，经过插值得到日解。三种负载引起的位移相加即为环境负载引起的测站地表形变。

步骤4，从全球数据中心SOPAC(Scripps Orbit and Permanent Array Center)获取ITRF2005框架下组合解，即“干净”的、去除趋势项的IGS基准站原始坐标时间序列，记为序列C。

步骤5，从步骤4获得的坐标时间序列中扣除步骤3获得的垂直方向位移时间序列，得经环境负载改正后的IGS基准站的坐标时间序列，记为序列C-B；从序列C-B中扣除步骤2所得的垂直方向位移时间序列，得经热膨胀效应改正后的IGS基准站坐标时间序列，记为序列C-B-A。

步骤6，采用频谱分析法确定序列C-B和序列C-B-A的周期，采用最小二乘拟合法估计振幅和相位。

序列的周期采用频谱分析法确定，假设有两个周期，分别记为D₁和D₂。采用式(9)所示的测站坐标运动经验模型，利用最小二乘拟合法估计各显著周期的振幅和相位信息。

式(9)中，p表示斜率；q表示截距；A₁、f₁分别表示周期D₁对应的振幅、相位和角频率；A₂、f₂分别表示周期D₂对应的振幅、相位和角频率；y_t表示IGS基准站在时刻t的坐标；ε表示噪声。

步骤7，采用数值指标量化热膨胀效应对GPS坐标时间序列的影响。

具体为：采用序列C-B和序列C-B-A之间的相关系数、加权均方根误差的改变程度PWRMS以及周期振幅的变化程度(即贡献量PE)三个指标，来量化热膨胀效应对GPS坐标时间序列的影响和贡献，结果见图3～6。

热膨胀效应对GPS坐标时间序列季节性信号的贡献评价指标如下：

(a)相关系数

计算序列C-B和序列C-B-A的相关系数，用于评价热膨胀效应引起的垂直位移和GPS垂向坐标时间序列间的相关性。两个相同长度的序列相关系数ρ_UH的定义为：

其中，n为序列长度，U_i和H_i分别为序列U和H中第i个元素，和分别为序列U和H的平均值。本发明中，序列U和H分别表示序列C-B和序列C-B-A。

相关系数可以用于评价两个序列的相关程度，相关系数为正代表两个序列正相关，为负代表两个序列负相关，为0代表两个序列互不相关。相关系数的取值范围为[-1,1]，相关系数的绝对值越大即相关系数越接近于1或-1，两个序列的相关程度越强；反之，相关系数越接近于0，两个序列的相关程度越弱。通常情况下，通过以下取值范围判断变量的正相关强度：若相关系数的取值范围为(0.8,1.0]，序列正相关性极强；若相关系数的取值范围为(0.4,0.8]，序列正相关强；若相关系数的取值范围为(0,0.4]，序列正相关性较弱。

(b)改正效果

改正效果采用加权均方根误差的改变程度PWRMS来量化。根据扣除热膨胀效应影响前后，WRMS(加权均方根误差)是否减小，以及减小的程度(记为PWRMS值)，可用来表征热膨胀效应对GPS垂直方向坐标时间序列季节性信号的改正效果。PWRMS值的定义为：

本发明中，WRMS(U)和WRMS(U-H)分别表示序列U和U-H的加权均方根误差，序列U和H分别表示序列C-B和序列C-B-A。

WRMS值用于评价序列观测值相对于真值的离散程度，通常用平均值近似代替真值。下面以序列U为例说明WRMS值的定义：

其中，WRMS(U)表示序列U的WRMS值，δ_u(i)表示序列U第i个元素的不确定度。

PWRMS值为正，意味着经过热膨胀效应的改正，垂直方向坐标时间序列的WRMS值减小，且PWRMS值越大，改正效果越好；PWRMS值为负则表示施加热膨胀效应改正会增加垂直方向时间序列的WRMS值。

(c)贡献量

根据扣除热膨胀效应影响前后垂直方向坐标时间序列的周期振幅，计算贡献量PE值，用于量化热膨胀效应对垂直方向季节性信号的贡献。PE值的定义为：

其中，k为周期数；PE(k)表示第k个周期下的贡献量；A_a(k)和A_b(k)分别为扣除热膨胀影响前后相应第k个周期的振幅。

如果PE值为正，则代表加入热膨胀改正后坐标时间序列的季节信号振幅减小，且PE值越大，代表热膨胀效应对于季节性振幅的改正效果越好；如果PE值为负，则意味着加入热膨胀改正后季节性振幅增大。

本实施例中，精确确定了热膨胀效应引起的全球560个IGS基准站垂直方向位移，并量化了热膨胀效应对GPS坐标时间序列季节性信号的改正效果和贡献量。全球560个IGS基准站中，由热膨胀效应引起的垂直方向位移时间序列的周年、半周年振幅最大分别可达5.83mm和0.51mm。同一IGS基准站半周年振幅量级较小，平均只占周年振幅的8.4％。周年及半周年振幅较大的测站主要分布在北美洲内陆及靠近极地的高纬度地区、欧洲中东部内陆区域和亚洲中部及东亚高纬度地区。

由模型计算的热膨胀影响时间序列与实际观测的GPS坐标时间序列的相关系数最大可达0.85，平均为0.32。其中，有41％的测站两序列强正相关，有50％的测站两序列弱正相关，只有9％的测站两序列负相关。扣除热膨胀效应影响后，有88.8％的测站垂直方向坐标时间序列WRMS值减小，最大可以减小3.43mm，改正效果最好可达27.0％。热膨胀效应平均可以解释17.7％和1.5％的GPS垂直方向坐标时间序列周年、半周年季节性信号。

由于环境负载影响较热膨胀影响更为显著，因此在量化热膨胀效应影响前应先予以扣除。结果表明，在待研究的49个IGS基准站中，扣除环境负载平均可使周年振幅降低1.09mm，在继续扣除热膨胀效应影响后，周年振幅总共可降低2.07mm，这意味着热膨胀效应平均可以使周年振幅降低约1mm，见图5。

为了更为直观的描述本发明中提出的量化指标的效果，图6给出了部分IGS基准站扣除热膨胀效应影响前后的垂直方向坐标时间序列，由上到下PWRMS值分别为21.3％、19.1％、18.0％。

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