华北平原地壳垂向位移季节性和时空变化分析方法和系统技术方案

本发明专利技术公开了一种华北平原地壳垂向位移季节性和时空变化分析方法和系统，该方法包括：获取地壳形变序列、GRACE卫星同期RL06球谐系数产品和降雨数据产品；解算得到地壳负荷形变；分解得到两个模态分量；确定两个模态分量各自的季节项分量和趋势项；根据两个模态分量各自的季节项分量和趋势项进行重构，确定季节项分和趋势项的重构结果；根据重构结果，确定华北平原地区地壳形变的时空变化关系，以实现对当地水资源的有效管理。本发明专利技术分析了华北平原地壳垂向位移与降雨量的时空变化关系，对当地水资源的有效管理以及居民生活具有重要意义。

【技术实现步骤摘要】
华北平原地壳垂向位移季节性和时空变化分析方法和系统
本专利技术属于卫星重力学、卫星导航学、水文学等交叉
，尤其涉及一种华北平原地壳垂向位移季节性和时空变化分析。
技术介绍
大气、海洋、冰冻圈和陆地水储量的重新分布及其相互作用会引起地表载荷的变化，从而使地壳在水平和垂直方向上产生形变。地壳形变主要分为构造形变与非构造形变，构造形变是由地球内部构造运动所引起的，主要表现为地壳在水平(N、E)方向的线性运动，非构造形变主要由大气负载、积雪负载、土壤水及海洋非潮汐等因素引起的，周年性的水圈运动会引起岩石圈季节性变化，主要表现为地壳垂向(U)的季节性波动。华北平原地处印度板块、菲律宾板块和太平洋板块之间，且该区域岩石结构在横、纵方向存在不均匀性，导致华北平原地壳运动活跃。此外，由于农业灌溉以及工业用水等人为活动对水资源的影响，该地区近几年出现了诸如地面沉降、海水入侵和水质恶化等一系列环境地质灾害。因此，对华北平原地壳垂向位移季节性和时空变化的研究极为重要。
技术实现思路
本专利技术的技术解决问题：克服现有技术的不足，提供一种华北平原地壳垂向位移季节性和时空变化分析方法和系统，分析了华北平原地壳垂向位移与降雨量的时空变化关系，对当地水资源的有效管理以及居民生活具有重要意义。为了解决上述技术问题，本专利技术公开了一种华北平原地壳垂向位移季节性和时空变化分析方法，包括：获取地壳形变序列yGPS、GRACE卫星同期RL06球谐系数产品和降雨数据产品；根据负荷弹性形变理论、滑动窗拟合多项式法、高斯滤波及GAC，对GRACE卫星同期RL06球谐系数产品进行改正，得到地壳负荷形变采用变分模态分解VMD方法分解地壳形变序列yGPS和地壳负荷形变得到yGPS和各自的模态分量IMFGPS和IMFGRACE；根据能量谱值和L-S谱分析法，确定IMFGPS的季节项分量SEAGPS和趋势项TRGPS、以及IMFGRACE的季节项分量SEAGRACE和趋势项TRGRACE；分别对IMFGPS的季节项分量SEAGPS和趋势项TRGPS、以及IMFGRACE的季节项分量SEAGRACE和趋势项TRGRACE进行重构，确定季节项分和趋势项的重构结果；根据重构结果，利用加权均方根法进行分析，确定SEAGPS与SEAGRACE的季节项关系，并结合TRGPS、TRGRACE和降雨数据产品确定华北平原地区地壳形变的时空变化关系，以实现对当地水资源的有效管理。在上述华北平原地壳垂向位移季节性和时空变化分析方法中，获取地壳形变序列yGPS，包括：获取中国大陆环境监测网络CMONOC提供的CORS站坐标时间序列，将CORS站坐标时间序列作为GPS卫星观测数据U(t)；剔除GPS卫星观测数据U(t)中大于三倍标准差的异常值，得到地壳形变序列yGPS。在上述华北平原地壳垂向位移季节性和时空变化分析方法中，根据负荷弹性形变理论、滑动窗拟合多项式法、高斯滤波及GAC，对GRACE卫星同期RL06球谐系数产品进行改正，得到地壳负荷形变包括：采用滑动窗拟合多项式法、高斯滤波及GAC，对GRACE卫星同期RL06球谐系数产品进行改正，得到地球重力场球谐系数C的变化量ΔClm和地球重力场球谐系数S的变化量ΔSlm；根据负荷弹性形变理论，确定负荷弹性形变公式：通过负荷弹性形变公式，解算得到地壳负荷形变其中，θ和分别表示余纬度和经度，A表示地球半径，hl表示l阶引力潮作用下的径向勒夫参数，kl表示l阶引力潮作用下的水平勒夫参数；Wl表示高斯平滑核函数，表示完全规格化缔合勒让德函数，m表示固定球谐系数的阶次，l表示勒夫参数的阶次。在上述华北平原地壳垂向位移季节性和时空变化分析方法中，其中，r表示高斯滤波平滑半径，a表示地球的平均半径，γ表示球面上任意两点与之间角度。在上述华北平原地壳垂向位移季节性和时空变化分析方法中，采用滑动窗拟合多项式法、高斯滤波及GAC，对GRACE卫星同期RL06球谐系数产品进行改正，得到地球重力场球谐系数C的变化量ΔClm和地球重力场球谐系数S的变化量ΔSlm，包括：确定GRACE卫星同期RL06球谐系数产品中的C20项，将GRACE卫星同期RL06球谐系数产品中的C20项采用CSR机构提供的激光数据文件C20_RL05.txt进行替换，并读取得到地球重力场球谐系数C的原始变化量ΔC和地球重力场球谐系数S的原始变化量ΔS；采用固定球谐系数的阶次m、大小为w的滑动窗口，对各阶球谐系数Cl-2i,m，…，Cl-2,m，Cl,m，Cl+2,m，…，Cl+2i,m；Sl-2i,m，…，Sl-2,m，Sl,m，Sl+2,m，…，Sl+2i,m进行二阶多项式拟合，得到对应阶次的第一拟合值C和第二拟合值S，将第一拟合值C和第二拟合值S作为改正值；其中，下标(p，m)表示阶数，p＝l-2i、…l-2、l、l+2、…、l+2i；获取GAC球谐系数C的变化量ΔCGAC和球谐系数S的变化量ΔSGAC，并依次通过如下公式解算得到地球重力场球谐系数C的变化量ΔClm和地球重力场球谐系数S的变化量ΔSlm：ΔC0＝ΔC-C，ΔS0＝ΔS-SΔC′＝ΔC0+ΔCGAC，ΔS′＝ΔS0+ΔSGACΔClm＝ΔC′+ΔC1，ΔSlm＝ΔS′+ΔS1其中，ΔC0表示去相关后的地球重力场球谐系数C的变化量，ΔS0表示去相关后的地球重力场球谐系数S的变化量，ΔC1表示地球重力场球谐系数C的一阶改正，ΔS1表示地球重力场球谐系数S的一阶改正。在上述华北平原地壳垂向位移季节性和时空变化分析方法中，采用变分模态分解VMD方法分解地壳形变序列yGPS和地壳负荷形变得到yGPS和各自的模态分量IMFGPS和IMFGRACE，包括：构造变分方程：其中，yGPS和均可作为待分解的序列，yk表示模态函数，wk表示IMF1-～IMFk两侧对称频率，k表示设定的模态分量个数，δ(t)表示均段脉冲函数，表示泛函数对时间t的一阶偏导，j为虚数单位；确定拉格朗日增广表达式：其中，α表示惩罚因子，λ(t)表示拉格朗日乘子；利用交替方向乘子法进行求解，得到变分解其中，和表示y(t)、yk(t)和λ(t)的傅里叶变换结果；将滤波后的信号进行傅里叶逆变换处理，当时，得到yGPS时域下的模态分量IMFGPS，当时，得到时域下的模态分量IMFGRACE；其中，IMFGPS＝IMF1+IMF2+…+IMFn或IMFGRACE＝IMF1+IMF2+…+IMFn。在上述华北平原地壳垂向位移季节性和时空变化分析方法中，根据能量谱值和L-S谱分析法，确定IMFGPS的季节项分量SEAGPS和趋势项TRGPS、以及IMFGRACE的季节项分量SEAGRACE和趋势项TRGRACE，包括：求解各IMF分量IMF1、IMF2、…、IM本文档来自技高网...

【技术保护点】
1.一种华北平原地壳垂向位移季节性和时空变化分析方法，其特征在于，包括：/n获取地壳形变序列y

【技术特征摘要】
1.一种华北平原地壳垂向位移季节性和时空变化分析方法，其特征在于，包括：
获取地壳形变序列yGPS、GRACE卫星同期RL06球谐系数产品和降雨数据产品；
根据负荷弹性形变理论、滑动窗拟合多项式法、高斯滤波及GAC，对GRACE卫星同期RL06球谐系数产品进行改正，得到地壳负荷形变
采用变分模态分解VMD方法分解地壳形变序列yGPS和地壳负荷形变得到yGPS和各自的模态分量IMFGPS和IMFGRACE；
根据能量谱值和L-S谱分析法，确定IMFGPS的季节项分量SEAGPS和趋势项TRGPS、以及IMFGRACE的季节项分量SEAGRACE和趋势项TRGRACE；
分别对IMFGPS的季节项分量SEAGPS和趋势项TRGPS、以及IMFGRACE的季节项分量SEAGRACE和趋势项TRGRACE进行重构，确定季节项分和趋势项的重构结果；
根据重构结果，利用加权均方根法进行分析，确定SEAGPS与SEAGRACE的季节项关系，并结合TRGPS、TRGRACE和降雨数据产品确定华北平原地区地壳形变的时空变化关系，以实现对当地水资源的有效管理。


2.根据权利要求1所述的华北平原地壳垂向位移季节性和时空变化分析方法，其特征在于，获取地壳形变序列yGPS，包括：
获取中国大陆环境监测网络CMONOC提供的CORS站坐标时间序列，将CORS站坐标时间序列作为GPS卫星观测数据U(t)；
剔除GPS卫星观测数据U(t)中大于三倍标准差的异常值，得到地壳形变序列yGPS。


3.根据权利要求1所述的华北平原地壳垂向位移季节性和时空变化分析方法，其特征在于，根据负荷弹性形变理论、滑动窗拟合多项式法、高斯滤波及GAC，对GRACE卫星同期RL06球谐系数产品进行改正，得到地壳负荷形变包括：
采用滑动窗拟合多项式法、高斯滤波及GAC，对GRACE卫星同期RL06球谐系数产品进行改正，得到地球重力场球谐系数C的变化量ΔClm和地球重力场球谐系数S的变化量ΔSlm；
根据负荷弹性形变理论，确定负荷弹性形变公式：



通过负荷弹性形变公式，解算得到地壳负荷形变
其中，θ和分别表示余纬度和经度，A表示地球半径，hl表示l阶引力潮作用下的径向勒夫参数，kl表示l阶引力潮作用下的水平勒夫参数；Wl表示高斯平滑核函数，表示完全规格化缔合勒让德函数，m表示固定球谐系数的阶次，l表示勒夫参数的阶次。


4.根据权利要求3所述的华北平原地壳垂向位移季节性和时空变化分析方法，其特征在于，






其中，r表示高斯滤波平滑半径，a表示地球的平均半径，γ表示球面上任意两点与之间角度。


5.根据权利要求3所述的华北平原地壳垂向位移季节性和时空变化分析方法，其特征在于，采用滑动窗拟合多项式法、高斯滤波及GAC，对GRACE卫星同期RL06球谐系数产品进行改正，得到地球重力场球谐系数C的变化量ΔClm和地球重力场球谐系数S的变化量ΔSlm，包括：
确定GRACE卫星同期RL06球谐系数产品中的C20项，将GRACE卫星同期RL06球谐系数产品中的C20项采用CSR机构提供的激光数据文件C20_RL05.txt进行替换，并读取得到地球重力场球谐系数C的原始变化量ΔC和地球重力场球谐系数S的原始变化量ΔS；
采用固定球谐系数的阶次m、大小为w的滑动窗口，对各阶球谐系数Cl-2i,m，…，Cl-2,m，Cl,m，Cl+2,m，…，Cl+2i,m；Sl-2i,m，…，Sl-2,m，Sl,m，Sl+2,m，…，Sl+2i,m进行二阶多项式拟合，得到对应阶次的第一拟合值C和第二拟合值S，将第一拟合值C和第二拟合值S作为改正值；其中，下标(p，m)表示阶数，p＝l-2i、…l-2、l、l+2、…、l+2i；
获取GAC球谐系数C的变化量ΔCGAC和球谐系数S的变化量...

【专利技术属性】
技术研发人员：郑伟，尹文杰，沈祎凡，
申请(专利权)人：中国空间技术研究院，
类型：发明
国别省市：北京;11