CN111797491B - 华北平原地壳垂向位移季节性和时空变化分析方法和系统 - Google Patents
华北平原地壳垂向位移季节性和时空变化分析方法和系统 Download PDFInfo
- Publication number
- CN111797491B CN111797491B CN202010365182.8A CN202010365182A CN111797491B CN 111797491 B CN111797491 B CN 111797491B CN 202010365182 A CN202010365182 A CN 202010365182A CN 111797491 B CN111797491 B CN 111797491B
- Authority
- CN
- China
- Prior art keywords
- gps
- imf
- grace
- crust
- spherical harmonic
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 230000001932 seasonal effect Effects 0.000 title claims abstract description 87
- 238000000034 method Methods 0.000 title claims abstract description 79
- 238000006073 displacement reaction Methods 0.000 title claims abstract description 28
- 230000001360 synchronised effect Effects 0.000 claims abstract description 19
- 230000008859 change Effects 0.000 claims abstract description 15
- 238000000354 decomposition reaction Methods 0.000 claims description 32
- 238000005259 measurement Methods 0.000 claims description 26
- 238000004458 analytical method Methods 0.000 claims description 21
- 238000012937 correction Methods 0.000 claims description 18
- 230000005484 gravity Effects 0.000 claims description 18
- 230000005489 elastic deformation Effects 0.000 claims description 17
- 238000001228 spectrum Methods 0.000 claims description 14
- 238000001914 filtration Methods 0.000 claims description 12
- 238000010183 spectrum analysis Methods 0.000 claims description 11
- 238000009499 grossing Methods 0.000 claims description 7
- 238000012544 monitoring process Methods 0.000 claims description 7
- 230000009471 action Effects 0.000 claims description 6
- 238000004364 calculation method Methods 0.000 claims description 6
- 230000007613 environmental effect Effects 0.000 claims description 5
- 101000872083 Danio rerio Delta-like protein C Proteins 0.000 claims description 4
- 230000002159 abnormal effect Effects 0.000 claims description 3
- 230000003416 augmentation Effects 0.000 claims description 3
- 230000002146 bilateral effect Effects 0.000 claims description 3
- 230000007246 mechanism Effects 0.000 claims description 3
- 238000012545 processing Methods 0.000 claims description 3
- 238000012216 screening Methods 0.000 claims description 3
- 208000028257 Joubert syndrome with oculorenal defect Diseases 0.000 claims 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 abstract description 20
- 238000011160 research Methods 0.000 description 9
- 230000033001 locomotion Effects 0.000 description 6
- 230000009467 reduction Effects 0.000 description 6
- 230000000694 effects Effects 0.000 description 5
- 230000002262 irrigation Effects 0.000 description 3
- 238000003973 irrigation Methods 0.000 description 3
- 238000012986 modification Methods 0.000 description 3
- 230000004048 modification Effects 0.000 description 3
- 238000002474 experimental method Methods 0.000 description 2
- 239000003673 groundwater Substances 0.000 description 2
- 238000011084 recovery Methods 0.000 description 2
- 239000011435 rock Substances 0.000 description 2
- 239000002689 soil Substances 0.000 description 2
- 102100022536 Helicase POLQ-like Human genes 0.000 description 1
- 101000899334 Homo sapiens Helicase POLQ-like Proteins 0.000 description 1
- 238000012300 Sequence Analysis Methods 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000006866 deterioration Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000001704 evaporation Methods 0.000 description 1
- 230000008020 evaporation Effects 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 239000008235 industrial water Substances 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 230000009545 invasion Effects 0.000 description 1
- 239000003621 irrigation water Substances 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 239000013535 sea water Substances 0.000 description 1
- 239000002352 surface water Substances 0.000 description 1
- 230000003313 weakening effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q10/00—Administration; Management
- G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Engineering & Computer Science (AREA)
- Business, Economics & Management (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Strategic Management (AREA)
- Human Resources & Organizations (AREA)
- Economics (AREA)
- General Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Marketing (AREA)
- Geometry (AREA)
- Game Theory and Decision Science (AREA)
- Evolutionary Computation (AREA)
- Computer Hardware Design (AREA)
- Entrepreneurship & Innovation (AREA)
- Development Economics (AREA)
- Operations Research (AREA)
- Quality & Reliability (AREA)
- Tourism & Hospitality (AREA)
- General Business, Economics & Management (AREA)
- Position Fixing By Use Of Radio Waves (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开了一种华北平原地壳垂向位移季节性和时空变化分析方法和系统，该方法包括：获取地壳形变序列、GRACE卫星同期RL06球谐系数产品和降雨数据产品；解算得到地壳负荷形变；分解得到两个模态分量；确定两个模态分量各自的季节项分量和趋势项；根据两个模态分量各自的季节项分量和趋势项进行重构，确定季节项分和趋势项的重构结果；根据重构结果，确定华北平原地区地壳形变的时空变化关系，以实现对当地水资源的有效管理。本发明分析了华北平原地壳垂向位移与降雨量的时空变化关系，对当地水资源的有效管理以及居民生活具有重要意义。
Description
技术领域
本发明属于卫星重力学、卫星导航学、水文学等交叉技术领域，尤其涉及一种华北平原地壳垂向位移季节性和时空变化分析。
背景技术
大气、海洋、冰冻圈和陆地水储量的重新分布及其相互作用会引起地表载荷的变化，从而使地壳在水平和垂直方向上产生形变。地壳形变主要分为构造形变与非构造形变，构造形变是由地球内部构造运动所引起的，主要表现为地壳在水平(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、…、IMFn的能量谱指数值，筛选得到第一个能量谱指数的极小值点，将所述第一个能量谱指数的极小值点之前的IMF分量确定为残差序列；
根据确定的残差序列，通过L-S谱分析法，判断IMF分量中季节项部分并对季节项部分进行重构，将IMF分量中除季节项部分之外的其余部作为即趋势项部分；其中，当IMFGPS＝IMF1+IMF2+…+IMFn时，得到IMFGPS的季节项分量SEAGPS和趋势项TRGPS；当IMFGRACE＝IMF1+IMF2+…+IMFn时，得到IMFGRACE的季节项分量SEAGRACE和趋势项TRGRACE。
本发明还公开了一种华北平原地壳垂向位移季节性和时空变化分析系统，包括：
数据获取模块，用于获取地壳形变序列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和降雨数据产品确定华北平原地区地壳形变的时空变化关系，以实现对当地水资源的有效管理。
本发明具有以下优点：
本发明公开了一种华北平原地壳垂向位移季节性和时空变化分析方案，融合了GPS、GRACE和降雨数据，利用GRACE改正GPS时间序列，研究GPS与GRACE序列在季节性方面的一致性关系，并结合气象数据网提供的降雨产品，分析了研究区地壳垂向位移与降雨量的时空变化关系，对当地水资源的有效管理以及居民生活具有重要意义。
附图说明
图1是本发明实施例中一种华北平原地壳垂向位移季节性和时空变化分析方法的流程示意图；
图2为BJFS测站GPS垂向时间序列VMD分解结果；
图3为华北平原GPS和GRACE垂向序列VMD分解所得季节项和长趋势项的示意图；
图4为华北平原地壳非构造形变与降雨量的时间序列的示意图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明公开的实施方式作进一步详细描述。
本发明公开了一种华北平原地壳垂向位移季节性和时空变化分析方法和系统，选取华北平原为研究区，利用全球定位系统GPS(Global Positioning System)、重力场恢复与气候实验GRACE(Gravity Recovery and Climate Experiment)和降雨数据分析了该地区垂向形变季节性特征和时空变化关系。
GPS具有高效率、高精度、全天候等特性。因此，国内外很多学者利用国际GPS服务IGS(International GNSS Service)和连续运行参考站CORS(Continuously OperatingReference Stations)对地壳形变进行监测并取得了显著成果，如利用GPS坐标时间序列分析亚欧大陆和喜马拉雅山地壳形变等。然而，GPS测站空间分辨率较低，不能对全球各个地区进行全覆盖。如何在保证精度的前提下提高地壳监测的空间分辨率，成为近些年的研究热点。2002年GRACE重力卫星的发射，开创了高精度重力场观测的新纪元。GRACE卫星可以监测到地面及以下所有深度的水储量变化，其中包括积雪、地表水、土壤水和地下水等组分。水储量的季节性变化会引起地壳负载的变化，导致地壳在U方向产生形变以及周围重力场的改变。联合GPS和GRACE研究地壳形变可以充分发挥两种监测手段的优势，该研究近些年在全球很多地区得到了广泛的应用。
在研究地壳形变在U方向的位移序列时，最常用方法是利用最小二乘拟合来提取信号中的季节项和趋势项，但该方法会遗漏信号中的隐匿值。鉴于上述方法的不足，学者进一步利用经验模态分解EMD(Empirical Mode Decomposition)研究地壳形变，该方法可以自适应地提取出信号中的季节项和趋势项，但在分解时会出现模态混叠和端点效应的问题。随着对EMD算法的优化，一些学者提出了集合经验模态分解EEMD(Ensemble EmpiricalMode Decomposition)、完备总体经验模态分解CEEMD(Complete Ensemble EmpiricalMode Decomposition)、变分模态分解VMD(Variational Mode Decomposition)等方法。虽然EEMD和CEEMD可以有效地抑制模态混叠的现象，但大规模的计算同时带来了冗余信息。康佳星等(2016)利用VMD方法对地震信号进行分解，结果表明该方法不仅分解精度优于EMD，并且分解效率也大幅度提高。本发明首次将VMD方法应用于GPS和GRACE序列信号分解的研究中，提取出序列的季节项和趋势项，并与EMD提取结果进行了对比分析。
不同于前人的已有研究，本发明首次利用VMD分解方法提取GPS和GRACE时间序列的季节项和趋势项。利用GRACE改正GPS时间序列，研究GPS与GRACE序列在季节性方面的一致性关系。并结合气象数据网提供的降雨产品，分析了研究区地壳垂向位移与降雨量的时空变化关系，对当地水资源的有效管理以及居民生活具有重要意义。
如图1，在本实施例中，该华北平原地壳垂向位移季节性和时空变化分析方法，包括：
步骤101，获取地壳形变序列yGPS、GRACE卫星同期RL06球谐系数产品和降雨数据产品。
在本实施例中，可以利用利用中国环境监测网络CMONOC(Crustal MovementObservation Network of China)提供的CORS站坐标时间序列进行分析：获取中国环境监测网络CMONOC提供的CORS站坐标时间序列，将CORS站坐标时间序列作为GPS卫星观测数据U(t)；剔除GPS卫星观测数据U(t)中大于三倍标准差的异常值，得到地壳形变序列yGPS。其中，华北平原共有10个CORS站点。
步骤102，根据负荷弹性形变理论、滑动窗拟合多项式法、高斯滤波及GAC，对GRACE卫星同期RL06球谐系数产品进行改正，得到地壳负荷形变
在本实施例中，可以采用滑动窗拟合多项式法、高斯滤波及GAC，对GRACE卫星同期RL06球谐系数产品进行改正，得到地球重力场球谐系数C的变化量ΔClm和地球重力场球谐系数S的变化量ΔSlm；然后，通过负荷弹性形变公式，解算得到地壳负荷形变
优选的，可以根据负荷弹性形变理论，确定负荷弹性形变公式：
其中，θ和分别表示余纬度和经度，A表示地球半径，hl表示l阶引力潮作用下的径向勒夫参数，kl表示l阶引力潮作用下的水平勒夫参数；Wl表示高斯平滑核函数，/>表示完全规格化缔合勒让德函数，m表示固定球谐系数的阶次，l表示勒夫参数的阶次。
优选的，本实施例利用300km的高斯滤波扣除南北条带误差，高斯滤波属于各项同性滤波的一种，其平滑核函数只取决于球面的两点/>和/>之间角度γ的函数关系，其中：
其中，r表示高斯滤波平滑半径，a表示地球的平均半径。
优选的，地球重力场球谐系数C的变化量ΔClm和地球重力场球谐系数S的变化量ΔSlm的解算可以如下：
确定GRACE卫星同期RL06球谐系数产品中的C20项。
由于GRACE卫星对C20项不敏感，遂本实施例将GRACE球谐系数数据中的C20项用CSR机构提供的激光数据文件C20_RL05.txt进行替换，并读取得到地球重力场球谐系数C的原始变化量ΔC和地球重力场球谐系数S的原始变化量ΔS。
本发明采用Swenson提出的滑动窗拟合多项式法扣除相关误差，以达到减弱阶次的相关性来削弱条带噪声的影响。具体的：采用固定球谐系数的阶次m(60次)、大小为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
由于考虑大气、海洋非潮汐因素的影响，在计算时应对重力场的球谐系数加入GAC改正，即在去相关后的地球重力场球谐系数C的变化量ΔC0和去相关后的地球重力场球谐系数S的变化量ΔS0基础上，加上GAC球谐系数C的变化量ΔCGAC和球谐系数S的变化量ΔSGAC，得到ΔC′和ΔS′：ΔC′＝ΔC0+ΔCGAC，ΔS′＝ΔS0+ΔSGAC。
由于卫星重力场测量的坐标系是地球质心坐标系，在上述得到的ΔC′和ΔS′的基础上加上地球重力场球谐系数C的一阶改正ΔC1和地球重力场球谐系数S的一阶改正ΔS1后，即可得到ΔClm和ΔSlm：ΔClm＝ΔC′+ΔC1，ΔSlm＝ΔS′+ΔS1。
步骤103，采用变分模态分解VMD方法分解地壳形变序列yGPS和地壳负荷形变得到yGPS和/>各自的模态分量IMFGPS和IMFGRACE。
Dragomiretskiy等提出了变分模态分解概念，它是一种递归自适应分解方法，可将复杂信号分解为具有中心频率和带宽的IMF分量。在迭代计算过程中，中心频率与有限带宽值持续变化。基本思路为先构造变分方程，再进行解分。
在本实施例中，构造如下变分方程：
其中，或/>yGPS和/>均可作为待分解的序列；yk表示模态函数，wk表示IMF1-～IMFk两侧对称频率，k表示设定的模态分量个数，δ(t)表示均段脉冲函数，/>表示泛函数对时间t的一阶偏导，j为虚数单位。
在变分方程基础上，利用拉格朗日因数和二次惩罚因子，将约束问题转化为非约束问题。信号中存在高频噪声，因此利用二次惩罚因数保证重构的准确性，使用拉格朗日因数确保约束的严格性，拉格朗日增广表达式如下：
其中，α表示惩罚因子，保证信号的重构精度；λ(t)表示拉格朗日乘子，使约束条件更为严格。
利用交替方向乘子法ADMM(Alternate Direction Method of Multipliers)进行求解，得到变分解yi n+1(t)：
其中，和/>表示y(t)、yk(t)和λ(t)的傅里叶变换结果。
将滤波后的信号进行傅里叶逆变换处理。其中，当时，得到yGPS时域下的模态分量IMFGPS，当/>时，得到/>时域下的模态分量IMFGRACE；其中，IMFGPS＝IMF1+IMF2+…+IMFn或IMFGRACE＝IMF1+IMF2+…+IMFn。
步骤104，根据能量谱值和L-S谱分析法，确定IMFGPS的季节项分量SEAGPS和趋势项TRGPS、以及IMFGRACE的季节项分量SEAGRACE和趋势项TRGRACE。
在本实施例中，求解各IMF分量IMF1、IMF2、…、IMFn的能量谱指数值，筛选得到第一个能量谱指数的极小值点，将所述第一个能量谱指数的极小值点之前的IMF分量确定为残差序列；根据确定的残差序列，通过L-S谱分析法，判断IMF分量中季节项部分并对季节项部分进行重构，将IMF分量中除季节项部分之外的其余部作为即趋势项部分。
其中，当IMFGPS＝IMF1+IMF2+…+IMFn时，得到IMFGPS的季节项分量SEAGPS和趋势项TRGPS；当IMFGRACE＝IMF1+IMF2+…+IMFn时，得到IMFGRACE的季节项分量SEAGRACE和趋势项TRGRACE。
以BJFS测站为例，利用VMD方法对GPS和GRACE信号分解处理，分解结果如图2所示。图2中，IMF1分量为提取出的序列长趋势项和季节项的叠加，其余部分均为高频残差时间序列。通过与EMD分解结果对比可知，利用VMD分解方法可有效去除序列的高频残差部分，获得序列的季节项(SEAGPS、SEAGRACE)与趋势项(TRGPS、TRGRACE)，并且可以削弱噪声对序列分析带来的影响。
步骤105，分别对IMFGPS的季节项分量SEAGPS和趋势项TRGPS、以及IMFGRACE的季节项分量SEAGRACE和趋势项TRGRACE进行重构，确定季节项分和趋势项的重构结果。
步骤106，根据重构结果，利用加权均方根法进行分析，确定SEAGPS与SEAGRACE的季节项关系，并结合TRGPS、TRGRACE和降雨数据产品确定华北平原地区地壳形变的时空变化关系，以实现对当地水资源的有效管理。
在上述实施例的基础上，下面对一种华北平原地壳垂向位移季节性和时空变化分析系统进行说明。
在本实施例中，该华北平原地壳垂向位移季节性和时空变化分析系统，包括：数据获取模块，用于获取地壳形变序列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和降雨数据产品确定华北平原地区地壳形变的时空变化关系，以实现对当地水资源的有效管理。
在上述实施例的基础上，下面对分析结果进行是说明。
GPS与GRACE监测结果比较分析
利用VMD分解方法，对10个站点的GPS和GRACE形变序列的趋势项和季节项进行提取，结果如图3所示。据图3可知，华北平原地区的地壳形变(GPS)与非构造形变(GRACE)在相位及振幅间存在一定差异。地壳形变序列的振幅在2.3～12.5mm之间，而非构造形变的变化范围为2.1～3.6mm。整体而言，GPS的形变量大于GRACE的结果，但HELQ和HAHB测站表现出GRACE形变量略大于GPS，究其原因可能是这两个测站位于华北平原西南部的边缘地区，此区域地壳构造运动不明显，地壳相对稳定，主要以地壳非构造形变为主。此外，地壳非构造形变呈现出显著的季节变化，主要是由于在春末和夏末时，地下水主要用于农业灌溉，而秋冬两季的雨水除了蒸发和径流外，其余部分会储存至第二年春季，故春季非构造形变量往往呈现负值的状态。
为了进一步验证地壳形变(GPS)与地壳非构造形变(GRACE)在季节项上的一致性，利用序列均值根误差WRMS(Wighted Root Mean Square)对其进行定量分析。
其中，n为每日解个数；ci为去趋势项后的时间序列；σ为标准误差；为GPS去除GRACE的时间序列；WRMSGPS和WRMSGPS-GRACE为GPS和去除GRACE后的加权均方根误差。
根据式(2)求出加权均方根误差后，同样也对照式(2)计算出GPS坐标时间序列扣除自身拟合后的WRMSGPS-GPSfit，利用式(3)定量评估GPS测站的WRMS小百分比WRMSreduction。
其中，WRMSreduction值可以表明序列之间在周期、振幅及相位等方面的关系。当WRMSreduction为1时，说明两者拟合得到的周年参数完全相同，利用式(9)计算得到各个测站的WRMSreduction值。对WRMSreduction值进一步统计可知，每个GPS测站的WRMSreduction值均在0.5以上，10个GPS测站改正平均值达到0.69；其中BJFS、TJWQ及HELY测站的WRMS改正量均达70％以上，说明GPS与GRACE周年结果一致性较好，且利用GRACE改正GPS行之有效。
华北平原地壳形变的时空变化特征
地壳垂向的非构造形变主要与水资源储量有关，当水量增多时，地壳产生向下的垂向形变；反之，地壳产生向上的回弹形变。为了分析华北平原地区地壳负载与非构造形变间的关系，本发明利用中国气象数据网CMA(China Meteorological Administration)提供的0.5度网格产品，提取出整个研究区的平均降雨量。利用VMD分解方法提取出GRACE反演地壳形变在各个时段的趋势信号，结果如图4所示。
据图4可知，在这13年间，华北平原的地壳非构造形变以0.20±0.07mm/yr的速率抬升。并且，地壳形变的季节性波动与降雨存在密切的联系。当夏季降雨较多时，地壳形变处于波谷；随着降雨的减少，形变量逐渐反弹，并在相位上存在2个月的滞后。2013-2015年地壳非构造形变呈明显的抬升趋势，斜率为1.66±0.62mm/yr，这期间的降雨量也明显低于多年的平均水平，表明降雨的变化特征对华北平原的地壳形变量具有重要影响，与之前研究结论一致。
整体而言，2003-2015年间华北平原北部地区降雨呈明显的增加趋势，而南部呈相反的下降趋势。南部地区水资源补给量减少，导致地壳在垂向上呈明显抬升趋势，这表明降雨量对地壳非构造形变具有重要影响；然而，降雨量下降趋势的最大漏斗位于华北平原南部(山东省中部)，而地壳抬升速率最快的区域位于华北平原的西南部(河南省境内)，两者空间分布差异表明农业灌溉等人为因素对西南部地区的垂向位移起主导作用。
本发明虽然已以较佳实施例公开如上，但其并不是用来限定本发明，任何本领域技术人员在不脱离本发明的精神和范围内，都可以利用上述揭示的方法和技术内容对本发明技术方案做出可能的变动和修改，因此，凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化及修饰，均属于本发明技术方案的保护范围。
本发明说明书中未作详细描述的内容属于本领域专业技术人员的公知技术。
Claims (2)
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和降雨数据产品确定华北平原地区地壳形变的时空变化关系；
获取地壳形变序列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机构提供的激光数据文件进行替换，并读取得到地球重力场球谐系数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、…、IMFn的能量谱指数值，筛选得到第一个能量谱指数的极小值点，将所述第一个能量谱指数的极小值点之前的IMF分量确定为残差序列；
根据确定的残差序列，通过L-S谱分析法，判断IMF分量中季节项部分并对季节项部分进行重构，将IMF分量中除季节项部分之外的其余部作为即趋势项部分；其中，当IMFGPS＝IMF1+IMF2+…+IMFn时，得到IMFGPS的季节项分量SEAGPS和趋势项TRGPS；当IMFGRACE＝IMF1+IMF2+…+IMFn时，得到IMFGRACE的季节项分量SEAGRACE和趋势项TRGRACE。
2.一种用于实现如权利要求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和降雨数据产品确定华北平原地区地壳形变的时空变化关系。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010365182.8A CN111797491B (zh) | 2020-04-30 | 2020-04-30 | 华北平原地壳垂向位移季节性和时空变化分析方法和系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010365182.8A CN111797491B (zh) | 2020-04-30 | 2020-04-30 | 华北平原地壳垂向位移季节性和时空变化分析方法和系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111797491A CN111797491A (zh) | 2020-10-20 |
| CN111797491B true CN111797491B (zh) | 2024-04-12 |
Family
ID=72806094
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010365182.8A Active CN111797491B (zh) | 2020-04-30 | 2020-04-30 | 华北平原地壳垂向位移季节性和时空变化分析方法和系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111797491B (zh) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113496083B (zh) * | 2021-04-13 | 2023-05-05 | 中国地震局地震预测研究所 | 一种gps流动站垂向速度场优化方法及装置 |
| CN116307853B (zh) * | 2023-02-21 | 2024-03-22 | 广东海洋大学 | 一种东南印度洋南赤道流的年际变化规律分析方法 |
| CN117556218A (zh) * | 2023-12-25 | 2024-02-13 | 北京建筑大学 | 一种线性工程季节性形变监测方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2015042754A1 (zh) * | 2013-09-29 | 2015-04-02 | 清华大学 | 一种低低星星跟踪卫星重力场测量性能解析计算方法 |
| CN107240228A (zh) * | 2017-05-23 | 2017-10-10 | 华东交通大学 | 一种基于倾斜仪和gps的火山监测与预警系统及方法 |
| CN110175214A (zh) * | 2019-02-01 | 2019-08-27 | 中国空间技术研究院 | 一种利用重力卫星数据监测极端气候变化的方法和系统 |
-
2020
- 2020-04-30 CN CN202010365182.8A patent/CN111797491B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2015042754A1 (zh) * | 2013-09-29 | 2015-04-02 | 清华大学 | 一种低低星星跟踪卫星重力场测量性能解析计算方法 |
| CN107240228A (zh) * | 2017-05-23 | 2017-10-10 | 华东交通大学 | 一种基于倾斜仪和gps的火山监测与预警系统及方法 |
| CN110175214A (zh) * | 2019-02-01 | 2019-08-27 | 中国空间技术研究院 | 一种利用重力卫星数据监测极端气候变化的方法和系统 |
Non-Patent Citations (3)
| Title |
|---|
| Research on Early fault diagnosis of rolling bearing based on VMD;Tao Zan等;2018 6th International Conference on Mechanical, Automotive and Materials Engineering (CMAME);地41-45页 * |
| 利用ＧＰＳ和ＧＲＡＣＥ分析四川地表垂向位移变化;丁一航;地球物理学报;第4777-4788页 * |
| 顾及GRACE季节影响的华北平原水储量变化反演;李圳;章传银;柯宝贵;刘阳;李婉秋;尹财;;测绘学报(07);第940-949页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111797491A (zh) | 2020-10-20 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111797491B (zh) | 华北平原地壳垂向位移季节性和时空变化分析方法和系统 | |
| Liu et al. | Spatiotemporal distributions of surface ozone levels in China from 2005 to 2017: A machine learning approach | |
| CN113297528B (zh) | 一种基于多源大数据的no2高分辨率时空分布计算方法 | |
| Birol et al. | Coastal applications from nadir altimetry: Example of the X-TRACK regional products | |
| Zhao et al. | Rapid glacier mass loss in the Southeastern Tibetan Plateau since the year 2000 from satellite observations | |
| Wei et al. | Evaluation of seventeen satellite-, reanalysis-, and gauge-based precipitation products for drought monitoring across mainland China | |
| Jiang et al. | Monitoring time-varying terrestrial water storage changes using daily GNSS measurements in Yunnan, southwest China | |
| CN108764688B (zh) | 基于星地多源降水数据融合的冬小麦湿渍害遥感监测方法 | |
| Shen et al. | Groundwater depletion in the Hai River Basin, China, from in situ and GRACE observations | |
| Yin et al. | Improved water storage estimates within the North China Plain by assimilating GRACE data into the CABLE model | |
| Tang et al. | Evaluating suitability of multiple precipitation products for the Lancang River Basin | |
| CN114201922A (zh) | 基于InSAR技术的动态滑坡敏感性预测方法及系统 | |
| Shen et al. | Feature extraction algorithm using a correlation coefficient combined with the VMD and its application to the GPS and GRACE | |
| Jiang et al. | Insights into hydrological drought characteristics using GNSS-inferred large-scale terrestrial water storage deficits | |
| Liu et al. | Analysis of groundwater changes (2003–2020) in the North China Plain using geodetic measurements | |
| Zhu et al. | Evaluation of a new satellite‐based precipitation data set for climate studies in the Xiang River basin, southern China | |
| Feng et al. | Groundwater storage change and driving factor analysis in north china using independent component decomposition | |
| Liao et al. | Numerical simulation and forecasting of water level for Qinghai Lake using multi-altimeter data between 2002 and 2012 | |
| Li et al. | Long-term and inter-annual mass changes of Patagonia Ice Field from GRACE | |
| Wang et al. | Tracking the source direction of surface mass loads using vertical and horizontal displacements from satellite geodesy: A case study of the inter-annual fluctuations in the water level in the Great Lakes | |
| Liu et al. | Unrevealing past and future vegetation restoration on the Loess Plateau and its impact on terrestrial water storage | |
| Li | Research on the Characteristics and Effects of Climate Extremes on Multi-spatial-temporal Scales in the Mongolian Plateau | |
| CN112285808B (zh) | 一种aphrodite降水数据的降尺度方法 | |
| Yang et al. | Dynamic evolution of recent droughts in Central Asia based on microwave remote sensing satellite products | |
| CN116068511B (zh) | 一种基于深度学习的InSAR大尺度系统误差改正方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |