中国地震局2008年以来陆续引入美国Microg-LaCoast公司的gPhone重力仪，截止目前我国86个重力站中有62个安装了该仪器．这些仪器为台风信号的检验、潮汐模型的检验、同震重力变化的观测、大震震源机制约束等科学研究[1-6]提供了丰富的观测数据．多年1 Hz采样的连续重力观测模式带来了系统元器件老化[5]．这会使观测和真值存在偏差．目前用一个无量纲的乘常数值(格值系数，scale factor，SF)来修正观测值和真实值之间的差异．并用相对误差(relative error，RE)来评估精度．长期潮汐观测重力仪的标定可用全球潮汐模型获得[7-8]，也可用带有基准的绝对或相对重力仪到被标仪器临近的观测室同址比测获得[9]．前者因标定精度不高只被应用在早期的弹簧重力仪中．后者具有原理简单、被标仪器无须移动的优点，在21世纪前后被广泛应用于超导重力仪的标定[10-11]．这些标定实例多为单站结果，比测策略和精度评估方法都因重力站研究需求不同而存在差异[12]．这都为确定我国gPhone重力仪(continuous relative gravimeter，CRG)观测能力，制定标定策略和观测数据高精度科学应用等带来了困难．2014年前，我国连续重力站仅关注单站潮汐相对变化的规律，只在重力站入网观测前在比测站进行标定[13]．此后由于重力仪搬运不便，且比测站标定格值是否适合其他重力站存疑，因此一直使用出厂格值分析数据．2014年后在重力站绝对重力测定时才兼顾解决连续重力站的长期观测带来格值系数变化问题．目前部分连续重力站开始进行了每年1次，每次1~3 d历时的比测工作．研究表明gPhone重力仪长期观测已可获得和早期超导重力仪相当的周日和半日波的潮汐观测能力[2]．然而中国地震局重力站多为野外站，在观测设备完备性以及比测时长的充分性上还无法达到超导重力站水平，限制了进行诸如海潮负荷等高精度潮汐观测信号的相关研究[12-14]．针对上述中国大陆gPhone重力仪占比最大，缺乏标定的问题，这里收集了2014~2021年我国12个重力站共计18点次的FG5和gPhone的同址观测资料，以及重力仪标定年份的重力固体潮和气压数据．在基于比测数据线性相关的基础上，提出理论潮汐模型约束的比测标定格值方法．与gPhone重力仪出厂格值、理论固体潮标定格值方法进行了重力残余振幅、重力残差振幅及gPhone仪器标称精度比较，分析了比测标定方法结果的可靠性和合理性．同时还讨论了比测策略对格值系数相对误差影响规律，推估gPhone重力仪获得优于0.001的格值系数相对误差所须具备的必要条件．研究结果可为利用弹簧重力仪进行潮汐模型的检验以及海潮负荷矢量的计算，以及未来我国重力台站的标定策略提供依据和参考．1 观测数据与标定方法1.1　比测标定数据及其预处理重力基准站会进行绝对重力测定．部分基准站还有gPhone重力仪进行高精度潮汐观测和模型建立．这些站绝对重力测定结果和gPhone重力仪潮汐观测结果共同组成了同址比测标定数据(简称比测数据)．选择其中绝对重力测定精度最高的FG5重力仪，在张家口等12个台站，共计18点次的比测数据进行gPhone重力仪的格值系数标定．表1为重力站比测的绝对重力测量下落次数、比测历时以及标定年份．10.13245/j.hust.220907.T001表1同址比测重力站基本情况重力站-年份下落次数比测历时/h重力站-年份下落次数比测历时/h襄阳-172 67826.28库尔勒-192 57825.28襄阳-162 49424.28太原-202 50024.28鹤岗-192 48224.28松潘-194 78624.77鹤岗-202 49624.28张家口-142 60025.27格尔木-202 49524.28张家口-182 49324.28阿勒泰-192 96229.28乌什-162 54326.28昆明-1912 32062.77乌什-172 45425.28蓟县-202 99729.28郑州-194 99924.78蓟县-213 49934.28郑州-202 59925.28受到激光问题、时钟漂移、环境和微震噪声影响，比测数据中会有极大偏离背景值、比测数据差值过大的时刻[15]．部分粗差可在线性回归迭代剔除，剩余部分还会对标定结果产生影响．为消除这些粗差影响，首先对gPhone重力1 Hz重力固体潮数据采用平均滤波器消除微震、地震和环境噪声的影响，然后利用gPhone的和FG5的比测数据的差进行粗差探测(图1)．10.13245/j.hust.220907.F001图1张家口站2014年05月17日比测数据预处理过程图1为张家口的2014年05月17日时段的同址比测标定数据，时间系统为世界时，图中：g为重力值；t为时间．gPhone重力仪和FG5绝对重力仪(absolute gravimeter，AG)的观测数据中都包含约为±15×10-8~±25×10-8 m/s2的噪声信号(图1(b))．利用平均滤波方法不仅降低了gPhone重力的观测噪声，而且改正了gPhone重力观测数据中的噪声产生的粗差(图1(a)和(b))．滤波值(low filter，LF)的噪声水平已经低于FG5的．利用AG和LF的差极易识别出绝对重力测量数据中的粗差(图1(b)和(d))．在预处理中，通过缺记粗差下落时刻的绝对重力数据(raw，简称R)或通过设定AG和LF的差的阈值范围预处理比测标定数据(preprocess，简称P)(见图1(c)和(d))．1.2　单站比测标定法和精度评定表1统计的比测时长显示比测数据时长最短为1 d，最长为2.5 d．这个时间跨度的gPhone重力仪零漂为线性，且和绝对重力仪相对变化一致．选择最小二乘法(least-squares，LS)估计格值系数及格值系数的(相对)误差．由于比测时段多数在1 d以上，如果标定年份重力固体潮观测数据中观测系统零漂大时就进行零漂改正，否则不考虑零漂．1.3　潮汐观测模型及其残差规律的检测单站的潮汐观测模型中各潮波的参数应该和全球潮汐模型值相当．选择DDW-NHi模型作为检验标准[2]．此外，标定后气压改正和固体潮改正后的重力残差矢量，以及标定后气压改正和海潮负荷改正后的重力残余矢量[9,11]，还应该符合重力残差振幅大于重力残余振幅，以及残余振幅应达到gPhone重力仪标称精度2个规律．对年度观测数据进行格值系数、气压和海潮负荷改正(海潮负荷改正模型选择FES.2004)，利用重力残差振幅与重力残余振幅规律检验1.1节比测数据预处理策略是否正确．不正确的情况下重复1.1节，1.2节过程使得重力站的潮汐观测结果符合规律，实现全球潮汐模型对标定结果的约束．1.4　理论潮汐模型的标定方法为和比测标定法进行比较，采用DDW-NHi模型作为参考，利用模型中观测振幅最大的M2波和O1波标定格值系数．以M2波标定方法为例，如fM2=δDDW-NHi/δM2，(1)式中：fM2为格值系数；δM2为观测的M2波潮汐因子；δDDW-NHi为DDW-NHi模型的M2波潮汐因子．2 标定结果和检验分析当分析标定结果的有效性时，还收集了Molodensky，DDW-He和Mathews全球理论潮汐模型[2]，同时利用观测潮汐模型的M2波和O1波潮汐因子作为标定参数进行标定．2.1　LS方法的标定结果比测标定计算的相对误差≤0.002 0时就能获得重力残余振幅≤0.7×10-8 m/s2的结果，若要进行全球潮汐模型或海潮模型检验则要求更高(≤0.001)[12]．图2为张家口站预处理前后标定结果，图中gCRG和gAG分别表示两个重力仪比测数据的变化值．预处理前比测格值系数和相对误差为1.011 6±0.002 1．和出厂格值系数相比有约0.011 6的差异．预处理后比测格值系数和相对误差为1.013 0±0.001 9．和预处理前比测格值系数有0.001 4的差异．经预处理相对误差从±0.002 1降至±0.001 9，这表明当前数据量难以估计出相对误差优于0.001 0的格值系数，可以达到0.002 0的水平．10.13245/j.hust.220907.F002图2张家口站2014年比测数据预处理前后的标定结果有研究认为高于M3波的最高频率(短周期部分)为仪器的高频噪声[10]．针对重力残余(residual gravity，RES)信号，取合成潮M3波的最高频率3.081 周期/d(频率单位为1 周期/d=1/86.4 kHz)为截止频率高通滤波，滤波值的均方根(root mean square，RMS)作为gPhone重力仪的观测噪声(见图3)．同时为分析重力残余的周日和半日波信号，还选择重力固体潮(earth tide，ET)周日波的最高频率1.000周期/d作为截止频率进行高通滤波，滤波值RMS为gPhone重力仪的观测精度(见图3)．10.13245/j.hust.220907.F003图3张家口站2014年标定结果的潮汐改正和检验和出厂格值相比，比测格值计算的仪器观测噪声均为0.24×10-8 m/s2，未发生明显变化(图3(c)和(d)．这说明在时域内格值系数改正对gPhone重力仪的观测噪声影响有限．与文献[10]iGrav-007 (0.038×10-8 m/s2)相比要低一个数量级．gPhone重力仪的观测噪声比超导重力仪大．这和文献[17]的结论也是一致的．从观测精度看，比测标定格值极大改善了残余信号周日和半日波信号的影响(图3(a)和(b))．出厂格值系数计算的残余信号中还存在1×10-8~2×10-8 m/s2的周日和半日波信号(图3(c))，比测标定结果这部分信号几乎全部改正(图3(d))．可见不同的标定方法对仪器的周日半日波信号的影响可达到1×10-8~2×10-8 m/s2．该影响已经达到了中国中部地区重力站(九峰站)的海潮负荷影响[18]，可见精密测定格值系数对海潮负荷的研究有巨大影响．2.2　基于理论模型标定结果检验除了与出厂格值进行比较之外，本研究还计算了基于M2和O1波的标定格值(见式(1))，并对表1中各重力站标定年份观测数据标定，气压和海潮负荷改正后计算潮汐观测模型．图4为观测模型、4个理论潮汐模型M2波，O1波和K1波潮汐因子按照重力站纬度顺序排列的比较．10.13245/j.hust.220907.F004图4观测潮汐模型和理论模型潮汐因子的比较图4(b)表明：O1波标定计算的潮汐因子普遍小于理论值．与出厂格值和比测格值相比，M2波标定格值计算的潮汐因子虽分居模型值两侧，然多数低于理论值(图4(a))．与比测格值相比，出厂格值计算的K1和O1波分居模型两侧(图4(b)和(c))，而M2波普遍偏小，且部分重力站偏差过大(图4(a))．只有比测格值的结果不仅分居理论模型值两侧，而且差异(≤1%)较其他标定结果(≤1.5%))更小．比测标定法计算的重力残差振幅和残余振幅见表2和表3．和出厂格值相比，比测标定法和M2波标定法重力残余振幅都减小了(表2)．和M2波的6次，O1波的12次大于仪器性能指标(1×10-8 m/s2)相比，比测标定法的重力残余振幅几乎都优于仪器性能指标．10.13245/j.hust.220907.T002表2不同标定方法计算重力残余振幅重力站-年份出厂格值法M2波标定法O1波标定法比测标定法乌什-170.620.410.040.30乌什-160.680.280.150.39库尔勒-191.130.370.270.23阿勒泰-190.560.420.100.12格尔木-201.690.740.470.28昆明-191.050.441.080.86松潘-190.100.191.940.08襄阳-172.272.182.621.66襄阳-162.222.132.611.53太原-203.081.121.950.69郑州-190.591.072.410.60郑州-200.731.142.490.47张家口-142.240.941.790.30张家口-182.870.811.730.94蓟县-200.210.441.450.23蓟县-210.631.031.510.66鹤岗-193.090.563.140.89鹤岗-202.550.593.140.9010-8 m/s210.13245/j.hust.220907.T003表3重力残差振幅和重力站经度关系重力站-年份经度/(°)重力残差振幅/(10-8 m∙s-2)M2波标定法比测标定法乌什-1779.210.380.35乌什-1679.210.330.35库尔勒-1986.170.470.45阿勒泰-1988.130.590.41格尔木-2094.870.440.45昆明-19102.750.870.92松潘-19103.600.970.79襄阳-17112.412.132.38襄阳-16112.412.052.39太原-20112.600.841.63郑州-19113.570.661.11郑州-20113.570.571.35张家口-14114.920.591.31张家口-18114.920.690.56蓟县-20117.531.401.64蓟县-21117.530.961.29鹤岗-19130.142.302.03鹤岗-20130.142.582.29将重力站按照经度至西向东排列，比较M2波标定和比测标定法的重力残差、重力残余振幅．比测标定法结果不仅满足重力残余振幅小于重力残差的规律(表2和3)，还满足重力残差振幅随着经度的增加而增大的变化特征(表3)．前者因为FES.2004全球海潮模型的负荷改正了残差矢量中的海潮负荷影响，后者则和我国东部重力站受到西太平洋海潮负荷的影响较突出有关．M2波标定法结果虽也具备一定的重力残差至西向东逐渐增加的特征，在110°≤经度≤120°的重力站残余信号比残差信号还大，这不符合海潮负荷改正规律，因此可认为利用比测法进行重力站观测数据的标定有助于中国大陆重力站海潮负荷影响的相关研究．2.3　比测标定结果综合比测标定法和理论模型法标定结果检验，各重力站标定年份格值系数和相对误差如表4所示．10.13245/j.hust.220907.T004表4比测标定法测定的格值系数和相对误差重力站-年份格值系数相对误差重力站-年份格值系数相对误差襄阳-171.006 70.002 7库尔勒-191.005 80.001 6襄阳-161.007 70.001 3太原-201.019 50.002 8鹤岗-191.014 90.003 8松潘-191.000 40.001 9鹤岗-201.011 60.003 4张家口-141.013 00.001 9格尔木-200.990 80.002 5张家口-181.012 40.002 7阿勒泰-191.003 00.001 4乌什-160.998 20.001 6昆明-191.001 30.001 5乌什-170.998 00.002 8蓟县-200.999 80.002 1郑州-191.000 00.002 3蓟县-210.999 80.001 3郑州-201.002 10.002 9和出厂格值系数相比，长期的连续重力观测表明，gPhone比测标定结果中格值系数发生了0.02%~1.90%的变化．该方法测定格值系数相对误差在±0.001 3~±0.003 8之间，虽无法达到高精度重力仪标定误差指标0.001[11]，但依然不妨碍在全球潮汐模型约束下检测到和海潮负荷有关的重力信号，而其他标定方法无法做到．3 讨论在理论固体潮模型约束下，利用比测标定法格值可以计算出符合中国大陆海潮负荷规律的重力残差结果．只是该方法还无法达到进行海潮负荷研究所需的≤0.001的相对精度[11]．为获得适合海潮负荷研究所需要的格值系数，从比测策略、最优比测时长两个方面分析影响格值系数相对误差的原因．3.1　比测下落数和标定精度的关系讨论为获得相对误差优于0.001的格值系数，绝对和超导重力仪比测时长应≥5 d[15]．现代iGrav型仪器可缩至3 d[10]．比测采用的测量策略是以小时为测组，每测组进行等采用间隔指定下落数测量．待该测组测完后到下一个测组开始，会有测组间断．这里利用2.1节方法分析比测数据中绝地重力测量下落次数、比测历时、测组间断时长和格值系数相对误差的关系，见图5．10.13245/j.hust.220907.F005图5比测策略和格值系数相对误差关系结果表明：当下次落数和比测历时增大时，格值系数的相对误差越小．当进行2 500~5 000次下落数，比测历时50 000~100 000 s就可达到0.002的相对误差．相对误差在0.002~0.001范围则以指数形式变化(图5(a)和(b))．测组间断时长的变化(1 800 s或3 600 s)不是影响相对误差的主要原因(图5(c))．3.2　最优标定时长的推估绝对重力基准资源有限，长期在野外进行连续重力站的比测会带来高成本挑战．进行时变重力研究时格值系数相对误差至少要优于0.002，而进行海潮负荷和全球潮汐模型的研究时则至少要优于0.001[12]．估计最优比测时长不仅可以降低观测成本，而且可以根据研究目的设计测量策略．以格值系数相对误差优于0.001为标准，分析张家口2014年和襄阳2016年的格值相对误差和比测下落数的关系．利用线性函数推估达到优于0.001相对误差所需要的比测下落数(见图6)．10.13245/j.hust.220907.F006图6格值系数相对误差和下落数的线性关系推估襄阳站(图6(b))下落次数在2 500~2 700次就已接近0.001．线性推估认为要达到比测结果相对误差优于0.001，需要下落次数在3 500~6 500次之间(约2~3 d)．张家口站(图6(a))相对误差较大，线性推估相对误差优于0.001就需要下落次数为 4 500~13 500次(约1.7~5.6 d)．这还仅是线性推估的结果．如果呈非线性的变化，那么比测时间还会更长．分析结果认为gPhone重力仪获得优于0.001的格值系数相对误差最佳条件下需要3~6 d，该比测时长大于iGrav达到相同相对误差所需时长[10]，和早期超导重力仪所需时长相当[15-16]．12个重力站18次的26~63 h的比测数据标定结果可见：比测数据仅能提供相对误差≤0.002的格值系数，不能提供可直接进行海潮负荷研究的相对误差≤0.001的格值系数．当利用全球理论固体潮时，重力残余振幅应小于重力残差振幅，且重力残余振幅应优于仪器标称精度规律的约束下，可获得符合海潮负荷至西向东逐渐增大规律的格值系数．