位置：
2024-12-31
中国地震科学实验场
研发团队：中国地震局地球物理研究所重力与地壳形变研究团队、中国地震局第一监测中心重力测量团队
一、模型研究的目的、必要性和需求分析
目的：建立川滇地区高空间分辨率的时变重力场模型，提供比GRACE (Gravity Recovery And Climate Experiment) 卫星模型更高空间分辨率的重力场时变信息，推进实验场范围内小尺度质量变化和动力学过程的研究进展。
必要性：2008年汶川地震之后，以中国地震局为主的国内外多家研究机构，在中国大陆开展了大量的多学科、多物理场的监测和探测工作。2018年5月，在汶川地震十周年国际研讨会暨第四届大陆地震会议上，国务委员王勇代表中国政府向世界宣布建设“中国地震科学实验场”。2021年国家“十四五”规划中，明确将建设实验场作为国家重大基础科研设施项目之一。建立川滇地区高空间分辨率时变重力场模型，不仅有助于深入理解地震科学实验场前沿科学问题，也可有力支撑国家重大基础设施建设。
需求分析：GRACE (Gravity Recovery And Climate Experiment) 和GRACE Follow-On为代表的重力测量卫星，可为全球大范围时变重力研究提供时变重力场模型。但卫星时变重力场模型的空间分辨率较为粗略（约300km），难以满足川滇地区科学研究对高分辨率时变重力信息的需求。构建川滇地区更高空间分辨率的时变重力场模型，可为该区域科学研究提供必要的基础重力变化信息，推动实验场范围内小尺度质量变化及动力学过程的研究进展。
二、模型预期指标和技术路线
预期指标：建立川滇地区高空间分辨率的时变重力场模型。
技术路线：（1）收集与整理2014至2017年川滇地区重力测网的流动重力观测数据，采用贝叶斯平差方法（Chen et al., 2019; 软件已开源发布https://gitee.com/cea2020/geoist）完成重力网平差。（2）发展基于Slepian局部谱分析理论的时变重力场建模方法，并利用这一方法和平差后重力观测数据构建川滇地区2014至2017年150阶时变重力场模型，对应空间分辨率约120km.（3）将模型低阶部分与独立的GRACE卫星重力数据进行比较，验证建模方法的有效性。
三、模型数据来源和质量分析
流动重力观测数据来源：川滇地区重力测网2014.09-2017.09共8期观测成果。每期观测都覆盖的测点称为重复测点，空间分布如图1所示。
数据质量分析：采用贝叶斯重力平差方法，完成对收集的2014至2017年各期流动重力数据的平差计算，估计了多台相对重力仪权系数、非线性漂移以及格值系数造成的观测数据的不确定性，给出了绝对重力控制下的各期的重力点值序列，各期点值平均精度在5.3-7.0微伽范围内，1 微伽= 1×10-8 m s-2。通过残差分析以及绝对重力交叉验证等方式评估了平差后各期数据的可靠性，其中与绝对重力观测偏差均值为12.9微伽。
图1 重复测点空间分布图。其中，蓝色十字为重复测点位置，黄色虚线为实验场边界，绿色线为Han等2021年构建的120阶时变重力场模型覆盖区域的边界，黑色线为川滇地区150阶时变重力场模型覆盖区域的边界，洋红色线为活动断裂，背景为该区域高程。
四、研制关键技术方法
Han等 (2021) 构建了川滇地区120阶时变重力场模型，空间分辨率约150km，模型覆盖区域如图1中绿色曲线所示，模型时间跨度为2015.09-2017.09。我们在Han等 (2021) 成果基础上，进一步发展局部重力场建模方法，构建川滇地区具备更长时间跨度（2014.09-2017.09）、更大覆盖区域（图1中黑色曲线）以及更高展开阶数/空间分辨率（150阶/空间分辨率约120km）的时变重力场模型。
我们采用Slepian谱分析方法来构建川滇地区时变重力场模型。Slepian谱分析的几个关键参数（最大展开阶数、最佳基函数数量以及抑制宽频带效应的平均半径）通过模拟实验得到。
最大展开阶数和最佳基函数数量确定方法
以时变水文模型WaterGAP Global Hydrology Model (WGHM)来模拟重力场变化，构建不同阶的重力场变化（作为已知场真值）。WGHM 模拟了全球范围的水流和水储量变化，除了给出全球范围的地下水 (groundwater) 储量变化外，以及其他9种水储量变化（canopy, soil moisture, snow, local lake, global lake, local wetland, global wetland, reservoir 和 river）。WGHM目前最新版本为2.2d，全球格网数据的间隔为0.5度（相当于360阶球谐展开）。
首先将WGHM给出的2014.04水储量变化转化为等效水高 (equivalent water height)，如图2a所示，单位为cm. 然后将等效水高变化转换为球谐系数，最后根据转换的球谐系数计算不同阶的重力变化，展开阶数分别为360, 180, 150, 120, 90 和60 （图2b-2g）。其次，根据不同阶的WGHM重力变化，模拟得到重复测点处不同阶N0的重力观测值 (N0=180, 150, 120, 90, 60)。
然后利用N0阶模拟观测值和Slepian方法（选取Shannon数作为最佳基函数数量）构建N0阶重力场，最后将建模结果与初始N0阶真值进行比较。令N0取360, 180, 150, 120, 90 和60，上述过程执行6次，获得一系列的建模结果；计算建模结果与初始真值的相关系数 (correlation coefficient) 以及二者差值均方根 (root mean square, RMS)，通过综合比较相关系数和RMS两组信息，确定研究区域的最大展开阶数为150.
图2 川滇地区用来模拟已知场的WGHM模型 (2014.04)：(a) 原始EWH结果，单位为cm；(b-g) 转换为360、180、150、120、90和60阶的重力变化，单位均为微伽。其中黑色线为150阶时变重力场模型覆盖区域的边界。
宽带泄漏的抑制方法
若观测数据中包含了超过最大展开阶数N的信号时，观测数据中高于N阶的信号会影响解算得到的低阶系数，造成宽频带泄漏 (broadband leakage)。采用如下步骤削弱宽频带泄漏的影响：
首先还是利用360阶WGHM结果（图2b）模拟重复测点处的重力变化，假设
为移除高阶信号后的新的重力变化
| (11) |
当迭代了250次后，半径变化范围为40-299km. 对每一次迭代的RMS进行分析，我们发现RMS最小值出现在r=124km处，即确定Gaussian空间平均的最佳半径为124km。
五、完成情况与成果
获得了川滇地区2014年9月至2017年9月150阶时变重力场模型，空间分辨率约120km，图3给出了各月的重力变化空间分布。
图3 川滇地区150阶时变重力场模型给出的重力空间变化（2014.04-2017.09），黑色线为时变重力场模型覆盖区域的边界。
六、验证（测试）与精度评价
我们利用GRACE卫星模型结果对150阶模型的低阶部分进行了比较。对比的阶数选为GRACE的最大阶数，即60阶。图4给出了两组模型重力变化速率的空间分布。在处理方法一致（均采用Slepian方法处理）和截断阶数一致（均截断至60阶）的前提下，图4a中重力变化速率与图4b中速率在研究区内都呈现出减小趋势且减小的中心位于康定-石棉-雅安一带，二者在空间分布上吻合较好，图4a中速率的变化幅度较图4b中稍大一些，可能原因是图4a在抑制宽频带泄露时，仍有残余的高频信号，影响到了所构建的60阶模型结果。
图4 川滇地区60阶重力变化速率（2014.04-2017.04）：(a) 由地面数据和Slepian建模方法所得结果；(b) 由GRACE数据和Slepian方法所得结果。黑色线为150阶时变重力场模型覆盖区域的边界，洋红色线为活动断裂。
七、模型使用说明
关键配套文件
文件名：研究区域边界.gmt
数据格式：两栏数据，第一栏为经度（°），第二栏为纬度（°）
备注：研究区域边界，研究区域外部没有观测数据，因此基于地面重力观测数据构建的150阶时变模型仅在研究区域内部有效。
川滇地区150阶时变重力场模型
文件名：YYYY-MM_dNNN_JOO.txt，其中YYYY为年份，MM为月份，NNN为展开阶数，不足三位则在前面补零，OO为截断数，不足两位则在前面补零，如2014年4月150阶截断数为24的模型，文件名为2014-04_d150_J24.txt
数据格式：三栏数据，第一栏为经度（°），第二栏为纬度（°），第三栏为时变重力场模型值（微伽）
参考文献
Chen, S., J. Zhuang, X. Li, et al. (2019), Bayesian approach for network adjustment for gravity survey campaign: methodology and model test, Journal of Geodesy, 93, 681–700.
Han, J., S. Chen, Z. Chen, et al. (2021), Determination of the degree 120 time-variable gravity field in the Sichuan-Yunnan region using Slepian functions and terrestrial measurements, Earthquake Science, 34(3), 211–221.
