超导重力仪能够精确监测地球表面重力场变化，它具有极高的灵敏度(可检测到振幅为0.2×10-9 m/s2的微弱信号)、极低的漂移率(一年只有几微伽)等优良特性，并且线性振幅响应频带范围为1 000 s到几年[1-3]，记录到的信号从数秒到数年这一很宽频带范围内的信噪比都非常高[4]．然而在将超导重力记录的数据用于科学研究之前，需要以0.1%的振幅精度和0.01 s的相位精度进行校准[4]．0.1%的振幅精度的确定是通过使用绝对重力仪进行同址观测获得的[5-10]；0.01 s的相位精度的确定是通过将正弦波或者阶跃信号输入传感器的反馈线圈中实现的[2,4]．例如，文献[4]通过对比利时Membach台站的GWR-C021注入正弦波和阶跃信号，获得了精度优于0.01 s的相位延迟．超导重力仪的主要用途是精确观测重力潮汐信号[1,11]．只有从超导重力仪的记录中准确去除潮汐信号后才能很好地获得其他信号，如大气、地下水、极移、地震及降雨等引起的重力变化[12-14]．在潮汐频段(周期从几小时到几年不等)，超导重力仪的振幅响应等于1.0，因此0.1%的格值因子精度足以精确预报局部重力潮汐．对于相位偏移(时间延迟)，通常在调和分析中不予考虑．因为若在调和分析中考虑相位偏移，则在计算重力合成潮(使用调和分析得到的潮波参数重新计算的局部重力潮汐，被用于从重力观测值中移除重力潮汐信号)时也必须考虑．但是相位偏移在用分析获得的潮汐因子进行地球动力学应用(例如研究液核自由章动)时，是必须要考虑的[15]．除了精确获得局部重力潮汐外，超导重力仪的观测记录还可进一步用于研究其他非潮汐信号．近年来，超导重力仪记录的背景噪声被用于研究高频信号，例如地球“嗡鸣”(Earth's hum)[16-19]．周期小于1 000 s时，超导重力仪的振幅响应与频率有关[2,4]，所以对于高频信号的研究，必须考虑仪器的振幅响应．振幅响应是超导重力仪传递函数的一部分信息，另一部分信息是相位响应．相位响应通常用以秒为单位的时间延迟来表示．由地球滞弹性引起的潮波时间延迟的准确估算至关重要，潮波的时间延迟可以用来解释地球内部的耗散机制[20]．尽管超导重力仪的传递函数非常重要，但仪器在出厂后未给买家提供仪器传递函数，因此超导重力仪在新台站进行安装后须要进行传递函数的测试实验．本研究通过输入阶跃信号来测定武汉台OSG-065的传递函数．实验步骤按照GWR公司(超导重力仪的生产厂商)提供的仪器操作说明进行．传递函数是通过TSoft软件的“Transfer function”命令获得的，使用快速傅里叶变换(fast Fourier transform，FFT)谱分析功能[21]．由于在数据处理过程中须要使用格值因子获取重力潮汐时间序列，因此在详述测定实验之前，首先简单介绍了OSG-065格值因子的标定方法；另外，由于所有超导重力仪重力潮汐观测数据都用GGP1滤波器进行低通滤波，也对该滤波器进行了简要介绍．1 确定OSG-065格值因子在进行传递函数的测定实验前，应精确测定超导重力仪的格值因子．OSG-065格值因子的测定是通过与FG5-112在2013年6月20日—24日期间的同址观测完成的．OSG-065和FG5-112记录的重力潮汐数据如图1所示，估算的格值因子为-92.392 7±0.037 3×10-8 m·s-2/V，相对精度为0.04%，优于0.1%，满足超导重力仪格值测定要求．10.13245/j.hust.220906.F001图1OSG-065和FG5-112观测记录2 GGP1滤波器超导重力仪记录的是重力采集卡的反馈电压．在获得数字信号之前，通常需要一个模拟低通滤波器起到防止频谱混叠的作用，以获得数字重力记录，GGP1滤波器就是这样一种抗混叠滤波器．国际地球潮汐中心(International Center for Earth Tides，ICET)认为该滤波器能满足全球地球动力学计划(Global Geodynamics Project，GGP)的要求[21]．从1997年7月1日开始，GGP组织了世界上大部分的超导重力仪台站，通过数据共享来研究全球地球动力学问题，此后GGP1滤波器被用于GGP所有超导重力仪上．该滤波器是由GWR 公司设计的8阶贝塞尔低通滤波器，其将频率大于0.5 Hz的信号衰减100 dB，拓扑结构为Sallen-Key，采样频率为1 Hz，拐角频率(-3 dB)为61.5 mHz，衰减(极限)为-160 dB/a，时间延迟为0.034 °(周期/d)-1，即在周日频段的相位延迟为0.034°，单位增益的1%范围为低于0.01 Hz的频段衰减量不超过±0.086 dB，单位增益的4%范围为低于0.02 Hz的频段衰减量不超过± 0.086 dB．武汉台超导重力仪OSG-065也安装了GGP1滤波器，所有记录到的重力潮汐数据都已采用该滤波器进行数据滤波．3 传递函数实验实验的所有操作均按照GWR公司提供的步骤执行．超导重力仪的输出信号分为 “G1-Sig”(未经过GGP1滤波)和“Grav-1”(经过GGP1滤波)两种类型．仪器引起的时间延迟可以通过对未经过GGP1滤波的信号来分析确定，GGP1滤波器引起的时间延迟可以通过计算经过GGP1滤波的信号与未经过GGP1滤波的信号的时间延迟之差得到．OSG-065传递函数测试实验于2015年5月6日实施，实验持续时间为UTC时间05:25:00到15:44:59，经过GGP1滤波和未经过GGP1滤波的阶跃响应数据如图2所示．通过使用格值因子标定后的数据减去合成潮汐计算得到重力残差序列，其中使用的格值因子是由上述绝对重力仪FG5-112经过同址观测实验获得．计算重力潮汐参数来合成潮汐，重力潮汐参数是通过对GWR-C032的长期(约14 a)观测资料进行调和分析得到．最后选择了31个间隔为10 min扣除潮汐后的阶跃响应信号(见图3)，单位由V转换为10-8m/s2．10.13245/j.hust.220906.F002图2经过和未经过GGP1滤波的OSG-065阶跃观测数据10.13245/j.hust.220906.F003图3OSG-065记录到的扣除潮汐后的阶跃观测数据4 传递函数计算方法按照GGP推荐的Michel Van Camp提供的TSoft中“Transfer function”命令计算传递函数[22]．为了使得时间延迟结果达到0.01 s或更高精度，至少需要对10个阶跃响应信号求取平均值，计算平均时所需要的阶跃响应信号的个数取决于台站的噪声水平．本研究对10个阶跃响应信号取平均来减少噪声影响，然后用Michel Van Camp的最小二乘滤波器(http://www.eas.slu.edu/GGP/phasecal.html)滤除高频噪声．在滤波过程中对未经过GGP1滤波的重力信号的截止频率设置为0.008 33 Hz，窗宽为61；经过GGP1滤波的重力信号的截止频率设置为0.2 Hz，窗宽为3．由于仪器对阶跃信号的响应主要集中在最开始的几分钟内，为了选择正确的时间段获得最佳估计，对于每个阶跃响应信号，将在阶跃响应开始后的2.0，2.5，3.0，3.5，4.0，4.5和5.0 min进行测试．运行“Transfer Function”命令后，出现“Transfer function (normalized amplitude)”和“Transfer function (group delay)”两个窗口．第一个窗口给出了用最低频率处的振幅来归一化振幅响应的结果；第二个窗口给出了以曲线光标所在处为阶跃响应起始时刻的群延迟(以s为单位)，其中大于拐角频率处的时间延迟是没有意义的．5 OSG-065的振幅响应基于未经过GGP1滤波的数据确定的OSG-065的归一化振幅响应如图4(a)和(b)所示，其中图4(b)是使用了最小二乘滤波后的结果．可以明显看出：对于周期小于1 000 s的信号，振幅响应是与频率有关的．比较图4(a)和(b)可以发现：最小二乘滤波后的结果比用原始数据计算的结果平滑得多，这表明最小二乘滤波可以改善振幅响应的计算结果．但是对于周期小于1 000 s的信号，经过最小二乘滤波后振幅响应曲线有明显的变化，因此不建议用最小二乘滤波器滤波后的数据来计算振幅响应．对比具有不同时间长度的阶跃响应信号的计算结果，表明阶跃响应信号的长度对振幅响应的计算几乎没有影响．10.13245/j.hust.220906.F004图4不同阶跃响应信号长度下的归一化振幅基于经过GGP1滤波的数据计算得到的归一化振幅响应如图4(c)和(d)所示，其中图4(d)是使用了最小二乘滤波的结果．可以清楚地看到周期小于1 000 s时振幅响应对频率的依赖性．对比图4(c)和(d)可见：最小二乘滤波的结果与未使用最小二乘滤波的结果类似，说明GGP1滤波器滤掉了高频噪声，最小二乘滤波对结果的影响可以忽略．6 OSG-065的时间延迟基于未经过GGP1滤波的数据计算得到的OSG-065的时间延迟如图5(a)所示．可以看到：对于周期小于4 000 s的信号，时间延迟是与频率相关的．虽然使用不同长度的观测数据，时间延迟达到恒定时对应的信号周期变化较大，但是在4 000 s之后时间延迟几乎是恒定的．所以最终选择了 8 000~16 000 s之间的时间延迟来确定OSG-065的时间延迟，并在图5(b)中进行了放大显示．图5(b)的结果表明：改变阶跃响应信号的长度会获得不同的时间延迟大小，阶跃响应信号长度取2.0 min时显得太短，但当使用更多数据时，计算中则可能会混入噪声，并且噪声会使得时间延迟的计算变得不稳定．基于图5(b)结果计算的OSG-065的时间延迟约为3.60±0.20 s，取的是去掉最大和最小值后的加权平均结果，结果不是特别准确，后续研究须采用更好的数学分析方法以获得更为准确的结果．10.13245/j.hust.220906.F005图5不同阶跃响应信号长度下未经GGP1滤波的时间延迟为了提高对时间延迟的估计，采用最小二乘滤波降低高频噪声对计算结果的影响．图5(c)显示了对应的8 000~16 000 s周期之间的时间延迟．比较 图5(b)和(c)可见：最小二乘滤波后时间延迟变得更加稳定，但其值从3.6±0.2 s变为18.81±0.01 s，说明此时最小二乘滤波不能用于时间延迟的确定，因为最小二乘滤波会使时间延迟发生改变．使用未经过GGP1滤波的数据确定的是从传感器到重力卡的时间延迟，大小为3.60±0.20 s．而使用经过GGP1滤波后的数据计算得到的OSG-065的时间延迟为10.16±0.01 s(图6(a))．两种结果之间的差异表明：GGP1滤波器会引起额外的约6.56 s的时间延迟．GWR公司给出的GGP1线性时间延迟为0.034 °(周期/d)-1，为8.16 s，与上面测定的结果6.56 s有差别，分析可能的原因有如下两个方面：未经过GGP1滤波的时间延迟3.60 s可能因为高频噪声的影响测定不准确；给定的GGP1滤波器的线性时间延迟不准确．10.13245/j.hust.220906.F006图6不同阶跃响应信号长度下经过GGP1滤波的时间延迟当数据经过GGP1滤波后，再使用最小二乘滤波对时间延迟的计算结果影响很小，如图6所示．这些数值结果表明：GGP1滤波器对数据具有时间延迟作用，且此时无须再对数据进行最小二乘滤波消除高频噪声的影响．也可以说明在数据降采样中采用最小二乘滤波对数据的时间延迟的影响非常小，这是在调和分析获得1 h采样数据后采用最小二乘滤波器的根本原因．7 结语本研究在国内首次通过在仪器中输入阶跃信号来测定武汉台OSG-065的传递函数(包括振幅响应和时间延迟)．结果表明：对于OSG-065记录的周期小于1 000 s的信号而言，其振幅响应与周期相关，这与前人研究结果一致[2,4]．为此，若要使用超导重力仪记录的周期小于1 000 s的信号，则须首先测定仪器的振幅响应．对于未经过GGP1滤波的OSG-065数据来说，其时间延迟为3.6±0.2 s，对于经过GGP1滤波后的数据来说，其值为10.16±0.01 s．对未经过GGP1滤波的数据使用最小二乘滤波会改变时间延迟，这表明使用未经过GGP1滤波的数据时不能再使用最小二乘滤波来计算时间延迟．但是当使用经过GGP1滤波后的数据时，最小二乘滤波对时间延迟的计算影响不大，两种情况下OSG-065的时间延迟均为10.16±0.01 s．因为OSG-065的观测数据都进行了GGP1滤波器滤波，因此数据分析时都须要考虑进行10.16 s的时间延迟改正．由于超导重力仪的优良特性，我国目前正布设多个超导重力仪台站，届时将会在国内组成一个超导重力仪观测网络．为研究背景噪声中的高频信号，其传递函数的测定是必不可少的工作，本研究为后续国内超导重力仪观测数据的分析应用奠定了基础．