陆地高精度重力观测可以获得近地壳内部场源的重力信号．通过对陆地定点重力场重复观测，获得的时变重力场信息，是监测地下介质物性变化的重要地球物理观测手段，也称为地震流动重力监测．近年来，地震流动重力观测不断地在年尺度地震危险性判定中发挥作用[1-4]．中国大陆地震重力观测网络系统是目前国际上最大规模的陆基重力观测系统，已建成101个绝对重力基准点，采用FG-5，A10等绝对重力仪每年观测一到两期，以及超过 3 000余个相对重力联测点，采用LCR-G，Burris，CG-5，CG-6等相对重力仪每年观测两期，该系统定期可以产出多种时空尺度的数据产品．此地震重力观测网络中测点之间距离一般十几km到几十km不等，两点之间观测时间间隔需要几十分钟至数小时，有限的观测数据为仪器参数的估计带来挑战．另外，其获得的重力观测数据信噪比较低，重力信号随时间变化量仅为十几个或几十μGal(1 μGal=1×10-8 m/s2)，是典型的μGal级重力信号[5-6]，因此须要对重力仪器参数的不确定性进行量化估计，并对数据处理算法的适用性开展相关研究．针对相对重力仪参数的估计问题，贝叶斯重力平差方法[7-8]将仪器的漂移率和格值系数作为超参数，采用二阶光滑作为先验假设，并以测网中已知的多个绝对重力点作为约束条件，引入赤池贝叶斯信息准则[9](Akaike's Bayesian information criterion)估计参数最优值．贝叶斯重力平差方法已应用到中国大陆地震重力观测网的多个区域，包括大华北测网的多网多期数据的处理和精度估计[10]，南北地震带南部多期观测数据的处理[11]，以及华南地区的累积重力场的获取与重力观测数据质量评估[12]．现有的贝叶斯平差方法体系中主要通过三种手段评估结果的可靠性：通过分析重力测量段差的残差时间序列是否存在规则性或周期性的信号成分，用于检验仪器的非线性漂移特性是否被剔除；通过测网内已知绝对重力点的交叉验证实验，可以检验平差后得到的重力点值的可靠性；通过测网内相邻重力测点的段差与仪器间互差之间是否具有相关性，判断使用相对重力仪器格值系数的准确性．但是，在川滇重力测网中，很多重力测点间距较大，地形复杂且交通不便导致相对重力仪观测每天观测的段差数量有限；因此，依据有限的仪器读数或采样，贝叶斯重力平差方法能否有效估计仪器参数和分离观测噪声，以及在某些段差观测不准或测点缺失的情况下，是否仍然能够稳健准确估计出仪器每日的漂移率特性是考察的关键．本研究围绕贝叶斯重力平差方法的稳定性开展研究，并给出了解决方案．针对模型参数反演的精度评定和统计估计值的偏差估计问题，Bootstrap重采样方法[13]和Jackknife重采样方法[14]应用较为广泛．其中，已有一些学者将Bootstrap方法用于地球科学反演精度评估问题，采用贝叶斯方法和Bootstrap重采样方法对震源参数的不确定性进行了误差分析[15]；通过Bootstrap方法对发震断层位置的不确定性进行评估[16]；应用Jackknife重采样方法解决加权整体最小二乘方法的精度评定问题[14]．利用Bootstrap方法实现震源参数的反演和精度评定[17]．引入Bootstrap和Jackknife两种统计学检验方法用于贝叶斯重力平差算法的参数估计稳定性研究中，结合四川和云南重力测网实例给出了稳定性评估方案，并开展模拟推估和实测数据应用，其结果可以为评价贝叶斯重力平差算法的稳定性和有效性提供新的解决方案．1 方法原理1.1　Bootstrap检验方法Bootstrap方法也称为自助法，是将原始样本作为总体，通过重采样代替求导计算，利用有放回的随机抽样法从总体分布函数中得到一系列独立样本，并通过计算每个独立样本来获取未知统计量抽样分布的经验估计[13-18]．该方法的优势性主要表现在它无须对未知模型及分布做任何假设，也无须推导估计量的精确表达式，仅须通过检验样本内统计量的变化便能够估计未知参数的均值和方差[19]．为检测贝叶斯平差方法在观测段差一定的情况下，有限采样的重力观测是否能有效恢复误差分布，基于实际重力测网的仪器读数来生成模拟观测数据进行验证．假设观测读数F的个数为n，根据流动重力中每个测点的观测噪声随机分布的特点[20]，生成一个样本容量为n的高斯噪声序列ε(εi为样本值)，采用Bootstrap方法重采样噪声序列，获取M个样本容量仍为n的噪声序列自助样本，每次将自助样本加到模拟观测读数F中，组成含有噪声的观测读数，并结合观测仪器类型，观测时间，模拟漂移率，观测点号与测点点位，生成模拟的流动重力观测文件．贝叶斯重力平差算法通过读取观测文件计算得到重力点值、段差残差、漂移率和仪器格值系数．贝叶斯重力平差算法的参数估计稳定性和精度评定的Bootstrap方法完整计算步骤如下所述．a．假设噪声序列N=(ε1,ε2,⋯,εn)为随机样本ε，将1/n的概率分别设置在N中的每个样本值上．b．产生服从均匀分布的n个随机数(i1,i2,⋯,in)，其中，1≤ij≤n，j=1,2,⋯，n；对N中的噪声值采样，获取自助样本Nr=(εi1r,εi2r,⋯,εinr)，其中1≤r≤M．c．将Nr带入观测读数，并组成观测文件，采用贝叶斯平差处理，获得4个未知参数的估计值，段差残差标准差σr，漂移率变化方差σb,r2，格值系数lr，点值偏差RGV,r．d．循环步骤b和c一共M次，得到M组4类参数的估计值(σ1,σ2,⋯,σM；σb,12,σb,22,⋯, σb,M2；l1,l2,⋯,lM；RGV,1,RGV,2⋯,RGV,M)；e．分别获得4个参数的Bootstrap方法估值均值，并与模拟真实值对比．1.2　Jackknife检验方法Jackknife重采样的思想来源于“leave-one-out”的交叉验证法．Jackknife方法是指从样本中每次删除一组数据从而得到Jackknife样本，由此生成的Jackknife样本集之间的差异很小，每两个Jackknife样本中只有两个原始样本不同．Jackknife重采样方法基于样本数据的大小，并且采样样本即为观测数据样本，与估计量分布无关，计算流程较简单[21]．将Jackknife方法应用于贝叶斯重力平差算法中，检验贝叶斯平差在部分段差样本丢失时，是否仍可以稳定估计模型参数，具体计算步骤如下．a．根据实际测网的观测路线和时间，模拟相对重力仪的观测文件，获取测网中n个测点的点号N=(u1,u2,⋯,un)．b．对点号文件进行Jackknife的采样，与原点号文件对比，获得被Jackknife法删除的点号ui；将观测文件中与点号ui同一行的相关的观测数据删除，剩余的观测数据组成Jackknife样本，将其作为新的观测数据样本．c．利用贝叶斯平差方法读取Jackknife法处理的观测样本，获得4个未知参数的估计值，段差残差标准差σr，漂移率变化方差σb,r2，格值系数lr，点值偏差RGV,r．d．循环步骤b和c一共n次，得到n组4类参数的估计值(σ1,σ2,⋯,σn；σb,12,σb,22,⋯，σb,n2；l1,l2,⋯,ln；RGV,1, RGV,2,⋯,RGV,M)，其中每一种参数估计值的序号为删除相应序号的点号的观测文件，在贝叶斯平差处理之后得到的结果．e．计算Jackknife法参数的估计值，即为n个Jackknife方法样本参数估值的均值；并与模拟真实值进行对比，测试准确性．2 模型测试在中国地震重力监测网中，四川测网和云南测网重力点间距分布不均，受交通条件限制测点间观测时间间隔差异性较大，但绝对重力基准点数分布较多，有利于用作贝叶斯重力平差结果的稳定性检验，因此选取这两个测网的实际观测时间和观测路线来设计稳定性评估所需的模拟数据集．如图1所示，四川测网采用LCR-G型相对重力仪观测165个相对重力联测点，同时含有3个绝对重力基准点(郫县，西昌，攀枝花)；云南测网包括192个相对重力联测点，采用两台CG-5型相对重力仪观测，有6个绝对重力观测点(昆明，丽江，下关，乡城，瑞丽，泸州)．设计两个测网的模拟观测文件，首先根据每个测网观测仪器实际漂移特性来设计模拟漂移率；并且考虑实际测量噪声的影响，从而在相对重力仪的模拟观测读数中加入标准差为3 μGal的高斯白噪声；每台仪器的模拟格值系数设定为1.000 000．Bootstrap重采样方法中的采样次数均设定为经验值200[19]．10.13245/j.hust.230305.F001图1四川和云南测网中测点和观测路线分布2.1　Bootstrap测试通过设计标准差均为3 μGal的观测随机噪声，采用Bootstrap方法对噪声序列分别进行重采样，并按照1.1节中的方法分别生成200个观测文件；分别采用贝叶斯平差方法重复计算，统计结果中观测残差标准差、格值系数、漂移率变化及点值结果．图2是基于Bootstrap重采样方法得到的四川和云南测网的格值系数和残差标准差的统计结果，图中：灰色箱型图表示Bootstrap重采样的模拟真实噪声；红色和蓝色箱型图分别表示贝叶斯平差的残差和格值系数；阴影区域表示数据的上下四分位数所包含的50%的数据；中间横线表示数据均值．四川测网通过Bootstrap重采样得到的格值系数和残差标准差的均值分别为0.999 998和2.99 μGal，与模拟真实值的均值偏差为2×10-6和0.01 μGal．类似的，Bootstrap重采样得到的云南测网格值系数和残差标准差均值分别为1.000 001和2.94 μGal，与模拟真实值的偏差为1×10-6和0.06 μGal．而且，Bootstrap重采样的模拟真实噪声标准差(灰色箱型图)，均在在贝叶斯平差后的残差标准差(红色箱型图)的范围内，观测噪声得到了较好地恢复．10.13245/j.hust.230305.F002图2Bootstrap方法重采样四川和云南测网模拟数据得到的残差标准差和格值系数的箱形图统计图3(a)和(b)分别给出了通过Bootstrap重采样得到的四川测网和云南测网的漂移率的采样估值及其范围，可以看出：两个测网的模拟真实漂移率均分布在采样估值的范围内，而且云南测网的漂移率采样估值均值与模拟漂移率基本一致．10.13245/j.hust.230305.F003图3Bootstrap方法重采样四川和云南测网模拟数据得到的漂移率和真实值对比最后，根据Bootstrap方法重采样200次得到的四川和云南测网所有测点的平差点值，计算其与模拟真实值的偏差值的绝对值的平均值，统计结果如图4所示．橙色实线区域内是四川测网西部受噪声影响较大的测点．10.13245/j.hust.230305.F004图4Bootstrap方法重采样四川和云南测网模拟数据得到的平差点值与模拟值的平均绝对偏差值(色标单位：μGal)可以看出：在四川测网的整体偏差在10 μGal以内，其中在测网西部(橙色实线范围内)，由于测环闭合时间较长，且无绝对重力基准点的控制，使得该测段受观测噪声变化影响较大，误差累积，导致点值偏差增大至9 μGal，而包含绝对点的东部和南部测网的点值偏差均在5 μGal以内；云南测网的整体偏差较小，在5 μGal以内，这与测网内6个绝对点的基准控制作用密切相关．2.2　Jackknife测试根据四川和云南测网的实测数据，模拟生成两个测网的模拟观测文件，采用Jackknife方法对两个测网的测点点号分别进行重采样，按照1.2中的方法分别生成与测点个数对应的165和192个观测文件，并采用贝叶斯重力平差方法分别处理得到平差结果．同样采用箱型图统计Jackknife方法得到的四川和云南测网的格值系数和残差标准差的分布情况，如图5所示．四川测网通过Jackknife重采样得到的格值系数和残差标准差均值分别为1.000 000和3.17 μGal，与模拟真实值的偏差为0和0.17 μGal．Jackknife方法重采样得到的云南测网格值系10.13245/j.hust.230305.F005图5Jackknife方法重采样四川和云南测网模拟数据得到的观测残差标准差和格值系数的箱形图统计数和残差标准差均值分别为0.999 989和2.97 μGal，与模拟真实值的偏差为1.1×10-5和0.03 μGal．相较于图2中Bootstrap方法得到的箱型图结果，图5中箱型图箱体更扁，数据更为集中，说明Jackknife测试得到的格值系数和残差标准差的估值样本离散度较小．Jackknife方法在采样过程中将噪声水平固定，而格值系数估计的主要干扰因素来自于噪声水平，使得方法得到的格值系数估值更稳定．Jackknife重采样方法获得漂移率估值如图6所示，相较于图3中Bootstrap得到的结果，Jackknife方法的漂移率估计范围明显较小，说明个别采样点的差异，对贝叶斯方法的估计结果影响不大，不影响稳定性．因为噪声水平同样会影响漂移率的估计，Jackknife重采样方法中噪声水平变化较小，使得漂移率的估计基本稳定．10.13245/j.hust.230305.F006图6Jackknife方法重采样四川和云南测网模拟数据得到的漂移率和真实值对比测点的平差点值与模拟真实值的偏差值的绝对平均值统计结果如图7所示．四川测网的整体偏差较小，在7 μGal以内，尤其是联测绝对点的测环偏差较小；而在云南测网当中，由于东西两侧测网仅以昆明绝对点相关联，因此当该附近出现测段缺失时，平差后重力点值会产生较大误差，偏离真实值可达9 μGal，同时在攀枝花附近的测段的段差较大，受测段缺失产生的误差累积较大(橙色实线范围内)．10.13245/j.hust.230305.F007图7Jackknife方法重采样四川和云南测网模拟数据得到的平差点值与模拟值的平均绝对偏差值(色标单位：μGal)3 实测数据检验在实际观测数据的检验过程中，采用Jackknife方法对观测数据进行重采样，来验证贝叶斯重力平差方法处理实测数据的可靠性和稳定性．采用Jackknife方法对四川和云南测网实际观测的测点点号进行重采样，从观测文件中删除对应的观测数据，生成重采样样本，并采用贝叶斯重力平差方法分别处理，得到平差结果．为了检验平差结果的稳定性，在没有应用Jackknife方法重采样的情况下，采用贝叶斯重力平差处理了四川和云南测网的原始观测文件，得到相应的平差结果，四川测网的格值系数和残差标准差分别为1.000 508和6.36 μGal，相应的云南测网的格值系数和残差标准差分别为0.999 959和4.58 μGal．图8给出了残差标准差和格值系数的Jackknife10.13245/j.hust.230305.F008图8Jackknife方法重采样四川和云南测网实测数据得到的观测残差标准差和格值系数的箱形图统计方法的估值分布情况．四川测网通过Jackknife重采样得格值系数和残差标准差均值分别为1.000 508和6.36 μGal，与未应用Jackknife重采样处理的原始观测文件得到的结果的偏差均为0．Jackknife方法重采样得到的云南测网格值系数和残差标准差均值分别为0.999 960和4.58 μGal，与未应用Jackknife重采样处理的原始观测文件得到的结果的偏差为1×10-6和0 μGal．与2.2节中的模拟测试的图5相类似，箱型图的形状较扁，数据较集中．进一步，对比了未应用和应用Jackknife重采样方法获得的四川和云南测网的漂移率，如图9(a)和(b)所示．可以看出：Jackknife方法的估值均值与原始观测得到的结果基本一致．同样地，与2.2中图6类似，噪声水平变化较小的情况下，Jackknife方法的漂移率估计范围明显缩小，不确定性减小，使得漂移率的估计基本稳定．其中：四川测网的漂移率较小，可能是与四川测网的大部分重力测段的段差较小和观测时间较短有关；云南测网的漂移率较大且逐渐增长，可能与该测网中测段的段差较大，以及交通不变导致观测时间较长有关．10.13245/j.hust.230305.F009图9四川和云南测网实测数据未应用和应用Jackknife方法重采样得到的漂移率对比最后，统计了Jackknife方法得到的两个测网的采样估值与原始观测的平差值的偏差，如图10所示，两个测网的偏差的绝对值的平均都在2 μGal以内，偏差较小．因此，对于四川和云南两个重力测网，当某些重力测段的观测结果存在粗差，贝叶斯平差仍然可以有效估计平差结果，稳定性较好．10.13245/j.hust.230305.F010图10四川和云南测网实测数据应用和未应用Jackknife方法重采样得到的平差点值之间的平均绝对偏差值(色标单位：μGal)通过Bootstrap重采样方法的模拟测试结果可以得出：当测点噪声标准差为3 μGal时，四川测网和云南测网在有限采样的流动重力观测中均能够有效的恢复误差分布，偏差在0.1 μGal以内，且格值系数估值均值偏差在1×10-5以内；其中四川测网西部因无绝对点分布，受噪声影响较大，建议在该地区建立绝对重力基准点，对流动重力观测起到基准控制的作用．Jackknife重采样方法的模拟测试表明，在部分段差缺失的情况下，噪声水平估计偏差在0.2 μGal以内，格值系数估值的均值偏差在1.1×10-5以内，漂移率估计在误差范围内．点值偏差分析发现，云南测网由于东西两侧测网之间只有昆明绝对点相连接，在部分段差缺失的情况下，在该测网连接地区存在整体的系统偏差，可以在该地区形成多个测环连接东西两侧测网以减小偏差；另外，攀枝花地区测段的段差较大，可以在该地区加密观测，以减小段差误差累积．对四川和云南测网实测数据的Jackknife重采样方法测试进一步验证了模拟测试的结论，贝叶斯平差的稳定性较高，重采样的点值估值的偏差均在2 μGal以内．综上所述，对于中国大陆典型的云南和四川重力测网中存在的重力测点稀疏、观测时间间隔非均匀问题，贝叶斯重力平差结果在Bootstrap和Jackknife两种统计检验过程中表现出了稳健特性．该方法可以结合具体流动重力测网分布和实际观测的时空采样特点，对观测误差与点值采样间距不同引起的结果不确定性进行量化评估，有助于评估重力平差结果的可靠性，为重力测网优化改造提供科学依据．