用于重力制图的量子传感

 

Ben Stray, Andrew Lamb, Aisha Kaushik, Jamie Vovrosh, Anthony Rodgers, Jonathan Winch, Farzad Hayati, Daniel Boddice, Artur Stabrawa, Alexander Niggebaum, Mehdi Langlois, Yu-Hung Lien, Samuel Lellouch, Sanaz Roshanmanesh, Kevin Ridley, Geoffrey de Villiers, Gareth Brown, Trevor Cross, George Tuckwell, Asaad Faramarzi, Nicole Metje, Kai Bongs & Michael Holynski (https://www.nature.com/articles/s41586-021-04315-3/metrics)

摘要

重力感应已成为地球物理学应用中的一种工具，例如工程和气候研究^1、2、3^，包括监测含水层^4^和大地测量^5^的时间变化。然而，由于去除振动噪声^6^需要很长的测量时间，使用重力制图来解决米级地下特征是不切实际的。本文通过设计量子重力梯度传感器克服了这一限制。该设计抑制了微震和激光噪声、热和磁场变化以及仪器倾斜的影响。该仪器的统计不确定度为 20 E (1 E = 10 ^-9^ s ^-2^ ) ，并用于在 8.5 米长的线路上执行 0.5 米的空间分辨率测量，检测信噪比为 8 的 2 米长隧道。使用贝叶斯推理方法，测得中心水平为±0.19米，中心深度为（1.89 -0.59/+2.3）米。消除振动噪声可以提高仪器性能，从而直接减少绘图时的测量时间。传感器参数适用于测绘含水层和评估对地下水位的影响7考古学^8,9,10,11^，土壤性质^12^和含水量^13^的测定，提供一个观察地下情况的途径，并降低在关键能源、交通和公用事业基础设施建设中不可预见的地面条件的风险^14^。

正文

量子重力梯度传感器使用原子干涉测量法^15^，它已在实验室实验中用于提供重力的灵敏测量^16^，研究等效原理^17^、精细结构常数^18^和牛顿引力常数^19^，促使人们渴望将这些传感器转变为在现实环境中使用的实用设备^20^。例如，已经创建了重力传感器，可用于火山和山地环境^21、22^，以及用于空中^23^、海上^24^和火箭^25^的测量。这些设备中，一种典型方法是使用光脉冲驱使原子中的双光子激发拉曼跃迁，并使之产生不同动量和能量状态物质波的叠加。由此产生的原子波包沿着两个空间分离的轨迹移动，然后被重新组合和干涉。这就产生了马赫-曾德干涉仪的物质波模拟。所得干涉图案中的相位差与局部重力场成比例。然而，与其他重力仪一样，由于需要微地震振动的平均值^26^，因此，此类设备的测量时间上会受到限制。这成为实现高空间分辨率的重力图的主要障碍。

为了实现重力制图及其在相关应用条件下的操作，我们实施了一个“沙漏”配置冷原子重力梯度仪^27^。这可以对由垂直基线^28^分开的两个原子云进行稳健的耦合差分测量。两个反向单光束磁光阱（MOT）允许普通拉曼光束通过以进行干涉测量（图1a）。对齐测量轴，以测量重力制图中最大和最相关的重力梯度张量（3×3）的垂直分量Gzz。差分操作抑制主要噪声源（振动和微地震）、系统偏移（例如倾斜）和用于驱动拉曼跃迁的光束之间的光路长度变化^29^。单光束磁光阱 6的冷却光束的激光强度噪声的共性使原子云温度波动稳定在几百纳开尔文内（图1b，顶部面板），限制交流斯塔克位移的影响，与传统的六光束方法相比，原子云质心运动降低了一个数量级（参见 方法）。由此产生的基线变化低于 75 ppm（图1b，底部面板），系统误差小于 0.1 E 。

图 1：沙漏梯度仪。

a，沙漏梯度仪，使用两个反向定向的单光束磁光阱，使用镜子组件（蓝色）实现。初始原子云（绿色）在局部重力加速度g下下降，然后受到时间间隔T的光脉冲以实现原子干涉仪（紫色）。光束传输用箭头表示（参见方法详情）。冷却光束（红色）被真空镜（蓝色）偏转以提供全方位的冷却，每个输入光束的中心部分穿过镜子之间的孔径，为对面的磁光阱提供最终的冷却光束。原子干涉光束（黄色箭头）具有较小的束腰，因此它们通过反射镜孔径而没有明显的削波。每个干涉仪同时操作，垂直基线间隔为 1 m。b，使用飞行时间和1m云分离基线全部时间的相对变化（下图）去测量每个捕获区域（上图）的原子云温度的时间变化^41^，（实线：在每次测量 4 s 时，50 次测量的 平均数据；阴影区域：平均数据的 σ 范围），根据抵达时间决定。C, 在靠近传感器（开放点）或远离传感器（填充点）两个位置之间，测量由测试质量在两个位置之间移动引起的重力梯度变化。每个测量数字代表测试质量的特定位置，奇数测量数字代表测试质量靠近传感器的数据。每个数据点由 8 次重力梯度测量的平均值形成，每个数据点包含来自原子干涉仪的 25 次发射，每次发射耗时 1.5 秒。每个数据点的误差条形图是八个重力梯度读数的标准误差。测试质量大约每 20 分钟移动一次，变化为 ±3.5 分钟，其位置可重复到大约 1 厘米。重力梯度信号变化的模型投影，ΔG zz, 显示为红色。

沙漏配置有几个实际的好处（见方法）。无需离轴光束，创建了一个坚固而紧凑的光学传输装置，在野外可以进行数月的操作，而无需校正对准。配置具有紧凑磁屏蔽，可以衰减25 dB从而抑制外部磁场的影响，阻止磁场影响原子云的生成。光束配置与坚固的全光纤激光系统相结合，能够独立控制反向传播的拉曼光束，促进光脉冲方向的反转^30^。在每个方向上的交叉测量会抑制若干系统效应，包括残差磁场效应引起的测量精度的降低。此外，寄生拉曼跃迁的相移和对比度损失阻止了通过每个方向独立传送拉曼光束^31^，因而无需锁相。

为了测量重力梯度（参见方法），在使用亚多普勒将Rb^87^冷却将云温度降低到微开尔文状态之前，每个磁光阱先加载 1 到 1.5 秒。接着，将原子云放下的同时进行原子干涉测量序列。使用荧光测量每个干涉仪的输出，以检测两个相关原子基态的种群比率，每个原子干涉仪中大约有 10^5^个原子，典型测量速率为 0.7 Hz。然后使用上部与下部原子干涉仪输出的利萨如图来提取不同的相位，从中确定重力梯度（图2)，插图）^32^。 在实验室条件下，通过调节传感器附近已知测试质量的位置以改变G zz来验证传感器（图1c）。与 202 E 的建模信号相比，会导致 (205 ± 13.1) E 的测量变化。

图 2：梯度计性能。

在室外运行约 8.6 小时期间，仪器输出的艾伦偏差（具有重叠平均值），以百分比误差显示。插图：用于艾伦偏差的椭圆数据（300 个点）的典型子集，显示了上下干涉仪输出信号的李萨如图，用于提取梯度相位。椭圆与圆的偏差会导致围绕椭圆极值点的点聚集。

同样，仪器的灵敏度和稳定性也在室外环境中进行了评估。相位数据的艾伦偏差^33^ (图2 ) 显示出 (466 ± 8) E/ $$\sqrt{H~z~}$$ 的平均短期灵敏度和在测量的 10 分钟内的统计不确定性为20E 。

为了展示重力制图的潜力，在现存的多用途隧道上方沿着 8.5 米测量线上进行了 0.5 米空间分辨率测量。该隧道的内部横截面为 2 米乘 2 米，其钢筋混凝土墙的墙壁大约 0.2 米厚。它位于两座多层建筑之间的路面下方。勘察现场附近的建筑物和地形提供了更多用来掩盖兴趣目标的信号^34^。为了估计隧道的预期信号，考虑了当地建筑物和地形，应用空气/土壤对比无限长方体空隙构建了现场模型，同时。该模型的参数是通过建筑平面图（CAD 文件），用探地雷达和现场测量的辅助数据（如地形扫描）的交叉检查而获取到的。隧道的估计峰值信号150 E ，而原子干涉仪的相位变化为17.5mrad。图3a显示了现场模型和原子干涉仪数据之间的比较，表明测量数据与预期的重力梯度异常与隧道的预期位置和尺寸一致。图3b展示了位置和隧道的比例表示，包括本地建筑物和位置拓扑。

图 3：实地调查。

a，具有标准误差（黑色散点）和总推断不确定性的重力梯度数据，包括模型不确定性（蓝色散点）、1 σ（深蓝色阴影）和 2 σ（浅蓝色阴影）可信区间带，以及场地模型（虚线） ）（见方法）。b , 现场比例示意图，显示隧道的尺寸和位置以及对重力信号的主要贡献。坐标原点（红点）在垂直方向上由测量线上的最低点定义，在水平方向上由隧道中心的预期位置定义。c，从重力梯度仪数据推断的开挖概率POE(probability of excavation)（蓝色轮廓）和预期的隧道位置（虚线）。d，通过转移焦点到使用已知的隧道几何知识和梯度仪数据推断的结果来获取土壤密度估计，显示有68%的最高密度间隔 (HDI)。

为了在实际应用中使用，以一种合理的方式解释数据和输出信息是非常重要的，从而能被用户所采纳和使用。为此，我们开发了一种贝叶斯推理方法，并将其应用于梯度计数据，并采用掩埋长方体^35^的数据生成模型为先验假设模型。使用梯度计数据以及对场地和地球物理参数的估计（详见扩展数据表2）来对深度和空间范围的异常进行定量预测。在此，假设在调查地点土壤类型密度在预期范围内，可以使用高斯分布，其平均值为 -1.80gcm^-3^，表示周围土壤中的空隙，标准偏差为 0.10gcm^-3^。推理过程使用开挖概率 (POE) 来度量^36^生成的隧道位置、深度和横截面积分布（图3c）。由于测量不确定性以及指定深度、面积和密度的模型参数之间存在的模糊性（特别是从可能实地数据推测出来的数据），可以预测观察到的 POE 是发散的^37^。用信噪比为8的信号对数据估计，得出隧道水平中心是沿测量线的(0.19 ± 0.19) m处，垂直中心在中心的 (1.89 -0.59/ +2.3) m处（见方法）。

此外，通过假设隧道几何形状的先验知识并包括调查地点的地形信息，推理的重点转移到推断土壤密度（图3d）。这导致密度参数接近高斯后验分布，平均值为 -1.80gcm^-3^，标准偏差为 0.15gcm^-3^。

原型仪器在静态操作期间展示出统计的不确定性（对于 20 E 梯度计的不确定性，相当于其两个重力仪中的每一个的 1.4 ng 不确定性）超过商业重力仪用于测量应用的报告性能 1.5 –4倍（参考文献38）。在使用量子重力传感器进行亚米分辨率映射的首次演示中，调查期间原型的可重复性与商业重力仪相似，并且受到系统效应的限制（见方法)，例如由于科里奥利效应，这可以通过进一步的工程来解决。此外，传感器可以在 75 秒内从一个空间位置移动到另一个空间位置，甚至可以垂直对齐到 1 毫度以内。如果解决这些方面的问题，例如通过在铁路或车辆上的操作，目前的仪器性能原则上将可以10 点线扫描和以信噪比为3实现15 分钟的总测量时间检测隧道，或者出现类似的异常。

隧道的检测用来评估仪器在一系列潜在应用中的性能。图4a显示了各种应用领域的典型信号大小范围与传感器的统计不确定性相比，当前仪器可检测到高于此范围内的特征。在土木工程中，这种特性可以减少地面条件的不确定性，并可用于检查城市中可开发地段、寻找隧道和大型或近地表公用设施，甚至检测落水洞的浸蚀特征。这种也可用于考古应用（例如，能够检测坟墓或隐藏的空洞，并探寻先前文明如何使用地下基础设施）。此外，该传感器在含水层测绘中特别有用，了解水其对环境的影响，进而能更好地优化水的使用。它还可用于测量地表的密度分布。根据推断的标准偏差基准，土壤密度提取方法取当前均值的 10% 范围变化，原则上这可用于区分干燥或饱和的土壤，或用于调查局部土壤压实（例如，在精准农业中）。这些应用的典型预期信号，以及我们的传感器在 10 分钟测量时间内实现的 20 E 统计不确定性，如图 1 所示。土壤密度提取方法目前对平均值中 10% 的水平变化敏感，这意味着原则上这可以区分干燥或饱和的土壤，或用于调查局部土壤压实（例如，在精准农业中）。这些应用的典型预期信号，以及我们的传感器在 10 分钟测量时间内达到的 20 E 统计不确定性，如图 1 所示。土壤密度提取方法目前对平均值中 10% 的水平变化敏感，这意味着原则上这可以区分干燥或饱和的土壤，或用于调查局部土壤结块状况（例如，在精准农业中）。这些应用的典型预期信号，以及我们的传感器在 10 分钟内测量的统计不确定性为 20 E ，如图4b所示。

图 4：应用相关性。

a，各种应用中典型重力梯度信号大小的等值线图，eötvös 中的等高线图的密度对比为1.8gcm^-3^。显示了检测到特征（红点）的参数和传感器统计的不确定性（实线）。实线上方区域中的特征用当前不确定性是可检测的。b，在一个区域未来的重力透视图用在一个区域内以 0.5 m 的空间分辨率使用，不确定性级别为 20 E。显示了一系列应用的预期信号大小。

与重力仪的情况相比，振动噪声的消除意味着未来仪器灵敏度的提高可以直接转化为测量时间的减少或不确定性的提高。对传感器进行进一步的科学改进，包括例如使用大动量分束器^39、40^，有可能使仪器的灵敏度进一步提高10到100 倍，从而可以更快地绘制或检测较小的和更多的特征。在未来 5-10 年内，预计这种性能实用仪器将会实现。

方法

实验概述

光脉冲原子干涉仪在概念上类似于光学干涉仪，只是光与原子的角色交换而已。作为物质波的原子受到一系列光脉冲的影响，这些光脉冲赋予它们动量，方式类似镜面反射和分束器。施加适当时长的光脉冲将导致原子的基态和激发态之间的跃迁，伴随着光子的吸收和受激发射。这种脉冲，通常称为 π 脉冲，由于传递动量而相当于原子光学镜。同样地，调整光脉冲通过向原子概率分布的一半进行能量刺激，使其只有 50% 的跃迁概率（通常称为 π/2 脉冲），因此也能用来充当分束器。通过在物质波脉冲上施加演化时间为T的 π/2- π - π/2的脉冲序列，可以创建光学Mach–Zehnder 干涉仪的等效作用，在脉冲之间。由此产生的由原子轨迹包围的时空区域（扩展数据图1b）与局部重力加速度成正比，然后可以从最终脉冲后两个原子态的相对布居数来测量。

重力梯度仪利用两个这样的干涉仪垂直偏移并用相同的脉冲序列同时探测。这抑制了共模效应，例如振动噪声或由于相对于地球重力的倾斜引起的相位变化，根据爱因斯坦的等效原理，这些效应与感兴趣的重力异常是无法区分的。我们的设备由两个子单元（扩展数据图1a）、一个传感器头和一个控制系统组成，光和电信号通过 5 米线传输。在扩展数据图2中总体显示了梯度计的尺寸、重量和功率特性。

传感器头具有真空系统，双 MOT 如同沙漏形式的配置，所有光通过轴上反向望远镜传递到原子。光在每个方向上传递，使用真空镜将部分光束重定向到原子俘获区域，在每个 MOT 中形成径向冷却光束。中心部分穿过，使得每个输入在给定方向为两个 MOT 提供垂直激光冷却光束。这使得径向冷却光束具有同样的波动强度（防止横向偏移），并且通过使用高斯光束形状，为垂直光束提供更高的强度，使该方向上的辐射压力更加饱满。提高了激光冷却过程的稳定性和稳健性，减少了温度或原子云位置的波动（图1b），不需要额外的激光源从而导致无法在野外操作使用的情况出现。在一个可比较的测试系统中，与六光束 MOT 在类似条件中在一小时内的平均云质心运动从（1.19 ± 0.86）毫米减少到（0.14 ± 0.09）毫米。两个 MOT 区域都有两个线圈，每个线圈由 92 匝 1 毫米 Kapton 涂层铜线缠绕在铝制模型周围（使用环氧树脂固定），并带有一个狭缝以防止涡流。线圈的半径为 43 mm，间距为 56 mm，在 2.5 A 的驱动电流下产生 12.5 G cm^-1^的线性场梯度。它们位于真空系统周围，因此其四极场的强磁场轴与冷却束轴的行进方向相同。此外，位于 MOT 区域周围的两组矩形线圈对，每组有 20 匝。它们间距100 毫米，垂直长度为 320 毫米，水平长度为 90 毫米，可用于补偿残留的磁场或应用的偏差。实际上， 粘团阶段没有使用场补偿。在下室中，使用一对线圈施加 0.63 G 场，在 MOT 阶段将原子云水平位置调整约 0.5 mm，提高干涉仪的对比度。偏置线圈^42^位于系统周围，以定义量化轴并消除磁子级之间的简并，与其他线圈在磁光捕获阶段后关闭。具有可变的间距形状以解决边缘效应并改善原子干涉测量区域上场的均匀性。该系统被封闭在一个磁屏蔽中，可提供 25 dB 的外部磁场衰减。原位磁场轮廓（通过拉曼跃迁光谱）是可测量的，在整个原子干涉测量区域内均低于 5%，并受到真空泵内部磁场源的限制。

激光系统由倍频至 780 nm 的电信激光器组成，靠近铷-87的D2线（参考文献^43、44^）。用于激光冷却的光是由激光器通过电光调制器 (EOM)输出，并从载波输出频率约为 1.2 GHz 的边带信号。提供一个锁定的信号用于使铷-85({|F}=3 > \to |{F}^{{\prime} }=4 > )转换来提供锁定信号，将载波频率设置为可在与 rubidium-87 的({|F}=2 > \to |{F}^{{\prime} }=3 > )跃迁的共振附近进行调谐，以提供冷却光。一个单独的 EOM 用于提供与({|F}=1 > \to |{F}^{{\prime} }=2 > )转换。原子干涉测量是通过双光子受激拉曼跃迁实现的。用于驱动这些的拉曼激光器的线宽为 73 kHz，并以 1.9 GHz 的偏移量锁定到({|F}=2 > \to |{F}^{{\prime} }=3 > \ )过渡。第二个拉曼频率是使用一对工作在 6.835 GHz 的 EOM 产生的。执行差分测量可抑制由于两个拉曼光束之间的光程长度变化而可能产生的相位噪声（例如由于光纤的振动和热引起的折射率变化而引起的相位噪声）。这允许两个光束独立传输，而无需在它们之间进行锁相，从而有助于实现仅将调制光谱应用于其中一个输入光束的实现。这避免了在使用传统的基于调制的方案（例如包括回射光束^31^的方案）时引起系统偏移和移相的寄生拉曼跃迁。实现时空反转的实际实现^30^，也称为波矢量反转，系统在拉曼光束的每个输入方向都有一个 EOM，并且在每次测量中将调制信号施加到一个臂上。这允许在测量之间改变赋予原子的动量反冲的方向，通过使用射频开关改变调制信号应用于哪个臂（参见扩展数据图1）。由于重力加速度对干涉仪相位的贡献对光所赋予的反冲方向敏感，而由许多其他影响（例如磁场引起的影响）产生的影响则不敏感。这允许在两个反冲方向上执行交错测量时消除这些影响。

使用保偏光纤将光传送到传感器头，冷却和拉曼光束使用单独的光纤。这些光纤将光传送到光学望远镜，将光准直到所需的光束尺寸。冷却梁的腰围为 24 mm，典型的最大功率为 130 mW。这些撞击在真空镜上，这是 15 毫米直角棱镜（Thorlabs，MRA15-E03），以提供水平冷却光束。这些反射镜以交叉配置安装在钛结构（使用 Epo-Tek H21D 粘合剂连接）上，使得它们的中心有一个 15 毫米的孔径。冷却光束的中心部分穿过这些孔，以提供对面 MOT 所需的第六光束。拉曼光束使用偏振分束立方体与冷却光束重叠，这样它们就会沿着与冷却灯相同的光束轴传递。每个拉曼光束的典型最大功率为 300 mW，其腰部设置为 6.2 mm，以限制真空镜中心孔径上的孔径和衍射，允许拉曼光束通过系统而不会被棱镜。尽管当前仪器中反射镜的孔径受到限制，但在更紧凑的系统或旨在进一步减少由激光束不均匀性引起的相移的系统中，可能希望使用比孔径更大的拉曼光束。如果追求这一点，则需要适当考虑来自孔径的衍射，因为在这种情况下，由于一个干涉仪从镜面反射中看到额外的光场，可能会导致进一步的光偏移。通过使用用于传递光的上下望远镜中的电压控制可变延迟板，将光的偏振设置为适当的配置，以冷却或驱动拉曼跃迁。使用来自光电二极管的反馈来主动稳定拉曼光束的强度，以控制也用于产生激光脉冲的声光调制器。

实验序列首先在 1-1.5 秒内从背景蒸气中收集每个 MOT 中大约 10 8铷 87 原子。然后使用糖蜜冷却将上云和下云温度分别降低至（2.86±0.09）μK和（3.70±0.20）μK（见图1b）。温度差异源于局部剩余磁场的差异，主要源于磁屏蔽几何形状，以及光学对准的微小差异。执行光学状态和速度选择以仅选择({|F}=1,{m}_{{\rm{F}}}=0 > )中的原子磁子水平和所需的速度等级。这是通过应用 π 脉冲和一系列吹走脉冲来去除处于不希望状态和速度等级的原子来实现的。然后以T = 85 ms的脉冲间隔 和 4 µs 的 π 脉冲长度执行原子干涉测量。使用双态荧光检测读取干涉仪以确定 | 的原子态总体比率。F = 2> 和 | F = 1> 基态，其中 (2.7 ± 0.1) × 10 5和 (1.7 ± 0.1) × 10 5原子分别参与上下干涉仪，典型测量速率为 0.7 Hz。通过将上部干涉仪输出与下部干涉仪输出作图来提取从中导出重力梯度的微分相位，以形成如图2插图所示的利萨如图。除了振动产生的随机噪声外，我们还故意在干涉仪的最终脉冲中添加了一个介于 0 和 2π 之间的随机相位值。在不对应于圆的椭圆相位处，即使对于均匀噪声，椭圆极值点周围的点簇也是可见的。

基于参与原子数的系统量子投影噪声约为 44 E/√Hz。总噪声预算包括来自其他术语的贡献，并显示在扩展数据表1中，以及调查期间观察到的相关系统。通过噪声过程的计算机模拟，与实验数据和椭圆拟合进行比较，研究了噪声预算。

测量数据的调查实践和处理

对于调查中的每次测量，通常在一个位置使用传感器头进行 600 次原子干涉仪运行（使用全站仪 Leica TS15 测量水平位置，距离路面大约 0.5 m 的垂直位置下传感器为 1.5 m，上传感器为 1.5 m），在每个干涉仪方向上给出 12 个 25 点椭圆，因此有 12 个单独的重力梯度估计值。对每个测量位置进行重复测量，通常每个位置有三个点。在每个测量点之间的基站进行测量，将一个位置的最终基站测量用作下一个位置的第一个。使用误差度量来确定每个椭圆的拟合质量，(\varepsilon )，定义为

$$\varepsilon =\frac{\left(\frac{1}{a}+\frac{1}{c}\right)}{2}{\left(\frac{1}{N}\mathop{ \sum }\limits{i=1}^{N}{L}{i}^{2}\right)}^{\frac{1}{2}},$$

其中N是数据点的数量，L是每个数据点与椭圆上一个点之间的最小距离，a和c对应于由方程参数化定义的椭圆(x=a{\rm{\sin } }\theta +b)和(y=c{\rm{\sin }}\left(\theta +\varphi \right)+d)分别。椭圆拟合中的误差对椭圆张角的变化很敏感^47^。 根据数值模拟，我们估计这种影响不到千分之几；因此，100 E 的变化将受到小于 0.5 E 的误差的影响。因此，与其他误差相比，此类误差很小。这种 100 E 的梯度变化将对应于椭圆形状的 11.6 mrad 变化。该相移可以与 2π 测量范围进行比较，这意味着该配置中仪器的测量范围与大多数感兴趣的实际特征相关（这些通常低于 400 E）。

发现具有(\varepsilon > 0.05)的椭圆拟合被自动丢弃。这导致所有数据中有 98.4% 可用于正常操作，与类似的传统地球物理设备相比，这代表了有利的数据正常运行时间。

为了处理数据，将一条直线拟合到基站点，然后从所有数据点中减去这条线。这是消除地球物理调查中漂移的标准做法。漂移的主要来源被认为是由于交流斯塔克位移，这也与两个云的温度差异有关。然后将重力梯度值作为测量点的平均值，从而估计测量位置和基站之间的梯度差异。此外，数据点的变化用于估计差值的误差。当组合来自同一位置的多个测量值时，使用加权平均值，从而为具有较大误差的测量值赋予较小的权重^48^。如图3a所示，数据未针对地形或潮汐等影响进行校正。潮汐效应没有得到校正，通过重力梯度的差分测量可以忽略不计。

在调查的测量位置发现的重力梯度误差的平均值为 17.9 E。将其与 150 E 的近似信号大小进行比较，得出的近似信噪比为 8。

从重力梯度仪数据推断

贝叶斯推理是一个框架，在该框架内，可以使用数据中包含的信息更新先验信念。对于模型参数向量 ( (\theta ) ) 和数据向量 ( (d) )

$$p\left(\theta |d\right)=\frac{p\left(d|\theta \right)p\left(\theta \right)}{p\left(d\right)},

其中(p\left(d|\theta \right))是可能性，(p\left(\theta \right))是先验，(p\left(d\right))是一个归一化常数，(p\left(\theta |d\right))是后验分布。

似然函数提供了测量数据(d)和从模型参数向量(\theta )计算的建模数据值之间的失配。这里使用的模型是一个三维长方体^35^；自由模型参数显示在扩展数据图3中，以及各个先验分布的函数形式。扩展数据表2详细说明了所选先验分布背后的基本原理。每个测量点的总不确定性是使用标准误差和模型不确定性随机变量的毕达哥拉斯总和乘以所有测量位置的标准误差的平均值来计算的。

概率 Python 包 pymc3（参考文献 49）用于实现长方体模型、定义模型参数先验分布和对后验分布进行采样，使用不掉头的采样器^50^。扩展数据 图4显示了选择模型参数的贝叶斯后验分布。

给定测量数据，参数后验分布表示关于模型参数的更新信念。为了帮助解释后验分布，在给定模型和先验分布的情况下，计算了 POE^36^，它表示地下异常的空间概率（如图3c所示）。隧道中心水平位置沿测线确定为（0.19±0.19）m，分布近似高斯分布。从使用测量线上最低点在垂直方向定义的原点到中心的深度为 (1.7 -0.59/+2.3) m。在隧道水平位置，从原点到地表的距离约为0.19 m，即地表到隧道中心的总距离为（1.89 -0.59/+2.3）m。从隧道几何形状来看，这将隧道顶部置于距地表约 0.89 m 深度处。

由当地特征产生的信号用于创建独特的场地模型。这用于提供对场地上重力梯度信号的预期形状的估计，以便与推理输出进行比较。这些特征包括感兴趣的隧道、附近建筑物的地下室、墙壁和排水管。它们显示在图3b的比例图中。