近地表流体迁移、地震、壳内物质运动和圈层相互耦合等不同动力学过程，都可能引起不同时空规模的重力场变化，已成为研究地球动力学过程的最基本和最直接的物理量[1]．陆地高精度微重力测量观测获取的时变微重力变化信号可用于研究地壳变形[2-3]、地下水储量变化[4-5]、冰川质量变化[6]、油气及藏矿勘探[7]、火山灾害[8]、地震及前兆[9-12]等科学问题．地表获取的时变微重力数据通常须要经过平差计算以获得测点的重力值，平差方法的选择及仪器参数的不确定性因素都是影响数据平差结果的重要因素，主要包括仪器漂移及格值系数．目前Lacoste-G型及CG-5型弹簧相对重力仪[13-14]在我国高精度陆地流动重力测量中发挥着重要作用，而CG-5型仪器由于其石英弹簧的迟滞性，日漂移量高达200 µGal[15]，尤其对于观测周期较长、空间跨度大，测网分布不均匀的云南测区，仪器漂移率变化并非线性．而仪器格值系数也会随着时间发生变化，有较大系统误差[16]，目前中国地震局重力学科组每4年组织一次长基线标定，但长基线标定不仅耗时耗财耗力，且标定周期难以满足每年两期的观测任务需求．经典平差方法通常采用最小二乘法进行平差计算[17]，但该方法并未考虑到仪器的非线性漂移特征，大量学者也提出不同的方法来消除仪器的非线性漂移，并估算仪器格值系数[18-20]，但都没有提出一套较为完整的鲁棒性较好的平差方案．而这里采用最新的贝叶斯重力平差方法[21-22]，考虑到仪器的非线性漂移，基于赤池贝叶斯信息准则(ABIC)[23]假设仪器漂移率光滑作为先验，从而对格值系数进行优化估计，同时获得仪器的漂移率变化特征并对不同仪器进行自动定权．而受到断裂构造、地形及交通道路等客观条件限制，地表流动重力测点分布并不均匀．此外，地表获取的时变重力信号很难进行精细分离，容易受到邻区局部场源及环境变化等因素干扰，这些不确定性因素都可能会影响重力观测系统的场源时空分辨率．因此，如何从非均匀分布的重力测网中构建合适的重力场源模型，并估算陆地流动重力测网的场源分辨能力，提取有效的场源变化信号，仍然是一个具有挑战性的难题．本研究以地震频发且构造活动强烈的云南地区为例，首先通过对比贝叶斯平方法[21-22](BAY)与经典平差方法(CLS)中仪器的段差残差、段差互差，并采用绝对重力交叉验证量化了相对重力仪漂移及格值系数等不确定性因素带来的影响；其次，采用球面六面体单元(Tesseroid模型)构建不同大小的重力异常模型检测板，通过等效源反演方法对异常源进行恢复，获取了测网的场源分辨能力．1 研究区概况云南位于青藏高原东南缘，受印度板块与欧亚板块的强烈碰撞密切，是中国大陆地震活动最强烈最频繁的地区之一．区域内构造活动背景复杂，强震频发，孕育了小江断裂、红河断裂、则木河断裂、安宁河断裂、澜沧江断裂、楚雄-建水断裂等主要活动断裂．测网经多次优化改造已形成迄今为止较为完整的、覆盖测区主要断裂的、绝对重力控制下的相对重力联测网络(图1)，共245个相对重力测点，10个绝对重力测点，其中最大段差约为305 mGal(10-5m/s2)，最小段差不到100 µGal(10-8m/s2)，平均点间距20~50 km．为去除季节性降雨变化影响，每年固定上下半年进行两期观测，采用两台CG-5型相对重力仪CG-5#1169及CG-5#1170(为简便，后文写为CG-169，CG-170)进行测量，观测周期约为90 d．10.13245/j.hust.220908.F001图1云南地区构造背景及重力测网分布图相对重力测量是测定地表两点之间的相对重力差值的过程，一般为沿一条测线的同点往返观测，而相邻两测点间的重力段差是测量的基本单元和重力网平差最直接的元素，段差残差可以直观地反映平差后段差与原始段差的差异，因此这里主要基于重力段差及段差残差进行精细化分析．2 方法及原理2.1　贝叶斯重力平差方法及原理由于陆地相对重力观测一般采用两台及以上的相对重力仪器同时进行，含有冗余观测，经典平差方法通常采取最小二乘法建立平差模型，而实际测量中，求解方程多是欠定问题，导致解的非唯一性，因此引入ABIC准则作为约束进行优化求解，通过求解ABIC最小化，在模型复杂度与模型对数据的描述能力之间寻求最佳平衡，从而避免过拟合．该方法将仪器噪声方差、漂移方差及仪器格值系数一并作为超参数加入平差方程中进行解算．首先假设观测误差、绝对重力测点误差及漂移误差均服从正态分布，则有Ax+Dv-y~N(0,σ2)，(1)Gx-g~N(0,σg2)，(2)Bv~N(0,σb2)，(3)式中：矢量x，v，y，g分别为待估重力点值、待求漂移率、测段观测段差及绝对点值；σ2，σg2，σb2分别为观测噪声方差、绝对重力观测方差及仪器漂移率方差；N(∙)为正态分布；矩阵A，D，G，B分别为观测顺序矩阵，观测时间矩阵，绝对点矩阵及二阶光滑矩阵．若令：S=ADG0，X=xv，Y=yg，W˜=W00Wg=σ-200σg-2，则有SX=Y，X由未知重力点值及漂移率构成，Y包括重力段差及已知重力值a1．该平差问题的最小二乘解则由下式求解U(X)=(SX-Y)TW˜(SX-Y)，(4)式中：Y=Δyl-ΔTα-ΔPβ；Δy为测段重力读数差；l为待求仪器的格值系数；ΔT为测段固体潮差值；α为潮汐因子；ΔP为大气压力的重力负荷差值；β为气压到导纳系数．通过求上述方程的最大似然估计，可以得到x，v的最优估计，多台仪器下相对重力仪和绝对重力观测的联合概率密度分布为L=[det(2πW˜-1)1/2]⋅exp[-(SX-Y)TW˜(SX-Y)/2]. (5)基于贝叶斯公式，后验概率函数则可表示为P=∫L(yv,x)ε(v)dv/∫ε(v)dv，(6)式中ε(v)为漂移率的先验分布，具体为ε(v)=∏i=1Pexp-12σbi2Bvi2．(7)采用赤池贝叶斯信息准则(ABIC)原理推导后得的最优值为μABIC=-2maxlog(L)+2H，(8)式中：L为似然函数；H为超参数个数，当对格值系数进行标定时，H=3P(P为仪器数)，且须代入至少两个以上独立观测的绝对重力点进行优化求解．2.2　球面六面体单元场源模型构建由于流动重力测网覆盖范围较广，达到几百km，而地球的曲率变化不容忽视，因此这里采用球六面体单元(Tesseroid)模型[24]来模拟场源，该模型基于球坐标系，而无须进行坐标转换，从而减小地球曲率影响及坐标投影变换带来的误差，适用于区域及全球空间尺度的三维模型构建．云南地区地震多发生于10 km左右，因此假设地下埋深10 km处存在正负密度相间的密度异常体，其等效密度异常为±1×10-3 g/cm3(地壳平均密度约为2.7 g/cm3)，场源模型的等效厚度为1 km(图2)．利用检测板模型实验，构建不同分辨率的重力异常场网格模型，再利用等效源反演方法[25]评估非均匀分布的重力测网对场源的空间分辨能力．10.13245/j.hust.220908.F002图2Tesseroid单元检测板模型3 测网数据质量分析这里重新整理了2018~2020年云南测区连续6期区域时变重力场数据，并以测点较多、测网较为完整的四期数据2018年第一期(2018-C1)、2019年第二期(2019-C2)、2020年第一期(2020-C1)及2020年第二期(2020-C2)数据进行精细化分析．3.1　仪器漂移率估计图3为经典平差方法(CLS)和贝叶斯平差方法(BAY)评估的两台仪器的漂移率变化特征．明显能看出在同一测量周期内经典平差方法的漂移率被视为常数，而贝叶斯平差获取的漂移率随着测量时间呈现出明显的非线性变化特征，且两台仪器漂移率变化存在显著差异，因此对于测量周期较长的测区须考虑仪器的非线性漂移特征．同时可看出CG-169仪器可能随着观测时间、地形、气压及路线等变化较为敏感，漂移率呈明显的非线性变化特征，而CG-170仪器的漂移率变化则相对稳定，基本呈现出稳态上升趋势，两台仪器漂移率目前基本稳定在30~40 µGal/h．10.13245/j.hust.220908.F003图3云南测网CG-5相对重力仪漂移率变化特征通过段差残差对比(图4)发现大部分残差控制在±15 µGal以内，且贝叶斯平差的段差残差波动小于经典平差．可看出经典平差的段差残差中带有明显的非随机变化信号，整体存在波动性，且这种波动性与图3中仪器的漂移率变化趋势相同，而贝叶斯平差后段差残差基本呈现为随机白噪声信号，与漂移无明显的相关性，也即是贝叶斯平差能去除由于仪器的非线性漂移带来的误差，更符合实际观测结果．同时可以看出：CG-169的段差残差总体在10.13245/j.hust.220908.F004图4贝叶斯平差与经典平差的段差残差时序图±40 µGal范围内波动，而CG-170仪器的段差残差相对较小，基本在±20 µGal范围内波动，即是CG-170仪器性能优于CG-169仪器．从段差残差标准差结果(表1)也可看出：贝叶斯平差后的标准差均小于传统平差的结果，且随着测量期次的增加，仪器逐渐趋于稳定，两种平差方法的段差残差标准差均得到提高，且贝叶斯平差的结果达到6.5 µGal．10.13245/j.hust.220908.T001表1段差残差的标准差对比6.5平差方法2018-C12019-C22020-C12020-C2经典平差(CLS)9.38.68.17.0贝叶斯平差(BAY)8.28.07.56.5µGal3.2　仪器格值系数估计应用贝叶斯重力平差优化模型，以测期准同步绝对重力为约束以估计格值系数．其中2018-C1初始格值系数为采用2018年初长基线标定的结果，2019-C2和2020-C2测期初始格值系数为云南测区内短基线标定的结果，其余测期初始格值系数均为采用该期实测数据进行标定．由于一年仅有一期绝对重力观测值，因此2018-C1和2018-C2贝叶斯格值系数估计均以2018年上半年的绝对重力观测为基准，2019-C1，2019-C2和2020-C1均以2019年上半年的绝对重力为基准，2020-C2测期以2020年下半年的绝对重力为基准．首先通过分析两台仪器间的段差互差相关性来确定格值系数标定是否准确，通常格值系数若准确，其相关性较弱．图5分别拟合了四期段差互差间的线性相关性，可见初始格值段差与互差间存在较为显著的线性相关性，而贝叶斯优化格值结果则显示互差相关性较弱，拟合线基本趋于水平．可见贝叶斯平差方法标定的格值系数较为准确，略优于长基线及短基线标定的结果．10.13245/j.hust.220908.F005图5段差互差图对测区2018~2020年6期格值系数进行贝叶斯优化估计并与初始格值系数对比得到表2．可见两台仪器贝叶斯格值优化前后的格值系数基本呈现同步变化，且不同测期优化后两台仪器的格值系数呈现不同程度的偏差，偏差最大值达32.8×10-5，最小值仅为2.9×10-5，平均偏差约为16.5×10-5．而整个测网观测的最大段差达305 mGal，即最大格值偏差将引起约90 µGal的误差，因此格值系数的偏差将可能影响整个测网的平差结果．同时发现：漂移率的变化对格值系数的标定也会带来影响，图3(b)中仪器漂移率变化较为平缓，表2中也显示2019-C1测期格值优化后偏差较小，而其余三期仪器漂移率变化明显不同步，且波动性较大，因此格值优化后偏差较大，最大达0.000 3．10.13245/j.hust.220908.T002表2初始格值系数与贝叶斯优化格值系数的偏差观测时间CG-169CG-170BAY格值系数格值偏差/10-5BAY格值系数格值偏差/10-52018-C12018-C22019-C12019-C22020-C12020-C20.999 9560.999 9911.000 1261.000 0971.000 0941.000 14729.730.0+3.75.513.918.80.999 7520.999 7380.999 8990.999 8350.999 8290.999 81331.432.8-2.92.912.117.43.3　绝对值交叉验证进一步结合测网独立观测获取的高精度绝对重力点值来验证格值系数的准确性．通过将一部分绝对重力测点代入平差计算，其余绝对点作为检验点，最后将解算得到的绝对重力点值与实测绝对重力值进行对比，以评估平差结果的优劣．以2020-C2为例，当选取昆明点为控制点时(图6(a))，对未参与平差的其余9个绝对点进行交叉检验，发现优化格值系数得到的绝对重力偏差较小．10.13245/j.hust.220908.F006图6经典平差及贝叶斯平差的绝对点点值误差进一步以丽江、下关、富宁为起算点(图6(b))，发现点值偏差值相比于仅以昆明为起算点的差值均有所减小，且贝叶斯方法减小更为明显，尤其对于在测网边界的控制点．最后采用9个绝对点以检测下关点时(图6(c))，发现随着控制点的增多，偏差越小，且贝叶斯平差结果的点值差均优于经典平差方法的结果，优化后的绝对重力点值差可减小约20 µGal的重力变化．3.4　点值结果对比最后对两种平差方法的点值结果进行比较得到图7，发现以昆明为起点，两种平差结果的点值差为-60~40 µGal；以昆明、下关和富宁为控制点时，点值差为-50~40 µGal；所有绝对点参与平差时，点值差为±30 µGal，两种平差方法获取的点值差随着控制点的增多而减小．且从图7(a)中不难看出：最大点值差主要集中在蒙自以南、攀枝花、乡城至丽江等观测时间较长或段差较大的区域，也说明了贝叶斯方法能在一定程度上减小漂移及格值系数不准带来的误差．10.13245/j.hust.220908.F007图7贝叶斯平差与经典平差点值结果对比图4 测网场源分辨能力评估设地下埋深10 km处存在正负密度相间的密度异常体，等效密度异常为±1×10-3 g/cm3，通过设计不同大小的网格单元模型，正演不同网格模型正演得到的地表实际测点处的理论重力异常场，由表3可知：当模型单元为1°×1°时，地表可观测到的重力异常约为±32 µGal，随着场源体的水平尺度的减小，可观测的重力信号越弱，即在设定场源参数下，由于重力异常的积分特性，不同尺度检测板场源体模型与对应地表可观测的重力异常量级呈正比．利用等效源反演方法，获取了无噪声下的场源密度结果(图8)，可见除测网边缘由于无测点分布而导致无法恢复场源参数外，测网基本具备0.75°×0.75°以上的场源监测能力；而0.5°×0.5°的场源反演结果显示，测网在点位分布密集的昆明周围小江断裂带附近及大理区域具有较好的场源分辨能力，而在测点相对较稀疏的滇东富宁至水城、滇西南的思茅附近、红河断裂带中段以及滇东北区域反演结果稍差，但整体基本具备0.5°×0.5°的场源分辨能力；0.25°×0.25°的反演结果显示，仅在测点极为密集的滇西大理测网及昆明小江断裂带附近附近有较好的分辨率，而其余测区结果较差，尤其在测网并不成环且测点稀疏的思茅-蒙自一带．10.13245/j.hust.220908.T003表3不同Tesseroid网格模型可观测到的重力异常Tesseroid网格模型模型正演重力异常/µGal最大值最小值平均值1°×1°32.02-31.92-1.370.75°×0.75°28.83-28.581.320.5°×0.5°23.24-23.000.140.25°×0.25°10.90-10.98-0.0510.13245/j.hust.220908.F008图8无噪声下的模型反演结果加入10%正态分布的随机噪声反演(图9)，可见在1°×1°和0.75°×0.75°的场源分辨率模型下，测网基本具备较好的测网分辨能力；在0.5°×0.5°的场源分辨率模型下，除蒙自-思茅一带由于测点相对较为稀疏，测网基本具备0.5°×0.5°的场源分辨能力；0.25°×0.25°的网格模型显示，仅在昆明小江断裂带附近及大理周边的部分区域由于点位分布密集，基本具备0.25°×0.25°的场源分辨能力．10.13245/j.hust.220908.F009图9加入10%高斯噪声下的反演结果通过研究得出以下主要结论：a．两台CG-5型相对重力仪呈明显的差异性漂移特征，CG-169漂移率随时间呈现明显的非线性变化，而CG-170漂移率变化则较为平缓；b．经典平差的段差残差中带有明显的非随机变化信号，而贝叶斯平差结果为随机白噪声信号，与漂移无明显的相关性，即贝叶斯平差能去除仪器非线性漂移带来的误差；c．相较于初始格值结果，贝叶斯优化格值系数获取的段差互差相关性明显减弱，即贝叶斯优化平差方法能获取较为准确的格值系数；d．贝叶斯优化后两台仪器格值系数基本呈同步变化，且格值系数变化与仪器漂移具有相关性；e．绝对重力交叉验证表明贝叶斯优化格值系数后的平差结果优于初始格值系数平差结果，优化后的绝对重力点值差可减小约20 µGal的重力变化，且随着控制点的增多，点值差逐步减小；f．在设定场源参数下，由于重力异常的积分特性，不同尺度检测板场源体模型与对应地表可观测的重力异常量级呈正比关系；g．加入10%的高斯噪声进行检测板模型测试，得出在测网主要块体边界带和断裂附近的场源分辨能力以0.50°×0.50°为主，而在昆明小江断裂带及大理红河断裂带北段等区域，由于测点分布相对密集，达到了0.25°×0.25°的场源分辨能力．贝叶斯重力平差方法也可推广到陆地区域重力测网的多网联合解算，提高联测解算精度，对评价陆地时变重力数据质量，研究地球内部介质迁移信号提供可靠的数据基础．通过对测网场源分辨能力的评估，有助于对重力测网监测能力进行评价进而优化测网布局，同时为后期更好地开展重力场地震监测、地震重力前兆信号分析、地壳深部质量变化和物质运移等研究提供了较好的数据及模型基础．