刘惠敏，刘繁明，荆 心，夏琳琳

（1. 哈尔滨工程大学 自动化学院，哈尔滨 150001；2. 青岛农业大学 机电工程学院，青岛 266109；3. 东北电力大学 自动化工程学院，吉林 132012）

浅地表地下质量异常重力与磁探测数值模拟方法

刘惠敏1,2，刘繁明1，荆 心1，夏琳琳3

（1. 哈尔滨工程大学 自动化学院，哈尔滨 150001；2. 青岛农业大学 机电工程学院，青岛 266109；3. 东北电力大学 自动化工程学院，吉林 132012）

提出一种重磁观测值数值模拟方法，用于地下空穴测量数据的解释和反演算法研究。采用更具结构灵活性的旋转椭球体作为地下空穴模型，在重力异常和重力梯度异常计算表达式的基础上，推导其磁异常计算表达式，实现空穴重磁异常信号的分析计算；借助扰动位的功率谱密度函数实现重力异常谱和重力梯度异常谱的计算，经傅立叶逆变换到空间域实现背景场统计学模拟；采用高斯白噪声模拟测量仪器或其它随机噪声信号。将三部分数据叠加并进行数值仿真，证明该方法合理有效。

重力异常；重力梯度异常；磁异常；地下空穴；重磁背景场；数值模拟

地下空穴是威胁城市、交通、矿区、大坝和核电站等重要设施安全、稳定运行的重要因素之一。作为地下质量异常存在的一种形式，地下空穴相对于周围环境存在负密度差和负磁化强度差，引起重力异常、重力梯度异常和磁异常。这些信号可由高灵敏度重磁仪器测量，是一种有效且非破坏性地表探测方法。

目前，高精度重磁测量技术、探地雷达技术等已经成功用于地下空穴的探测和圈定过程中[1-5]。其中，对于已测得的重磁数据，探索和发展适用于地下空穴的数据解释理论和反演算法已经成为判断空穴存在状态和物性参数的迫切需要，决定测量工作的成效。然而，由于测量设备、测区环境、测量人员、改正算法等因素的影响，实际工作中很难得到大量准确的、足够精度的测量结果用于数据解释和反演算法研究。

Moritz指出，在一个平面或是球体上，重力场可以看作是一个随机过程，在去除长波长的情况下，可以保证该随机过程的平稳性[6]。Jekeli指出，在统计学中心极限定理的基础上，若局部重力场可以看成多个随机场的叠加，则可以假设其在局部区域范围内满足高斯性[7]。尽管学术界对于采用模拟测量值的有效性存在不同观点，但在大地测量学领域，多个应用已经证明采用数值模拟技术是一种合理可行的近似方法[8-9]。如Jekeli采用零均值的高斯平稳随机过程作为重磁背景场的模拟模型[10]，并验证了该模型的正确性[11,14-15]。Abt等采用数值模拟模型讨论了匹配滤波器在地下质量异常探测中的有效性[11]。本文针对天然地下空穴在测量数据解释和反演算法研究中的需要，采用统计学方法构建重磁模拟背景场和随机仪器噪声信号，提出采用更具结构灵活性的旋转椭球体取代棱柱体或球体作为地下空穴的模型，满足重磁观测的数值模拟要求。

1 重磁观测模拟值的组成

假设测量区域位于零高程水平面上，如图1所示。测区下方存在三个以P(x1,y1,z1)、Q(x2,y2,z2)和R(x3,y3,z3)为中心点的空穴。某测点的观测值Z(x,y,z)由三个部分组成：

式中：S(x,y,z)表示地下空穴产生的重磁异常信号；B(x,y,z)表示重磁背景场信号；N(x,y,z)表示仪器或其他随机噪声信号。

图1 地下空穴在观测区域的位置分布Fig.1 Distribution of underground voids in observation area

2 地下空穴产生的重磁异常信号模拟

采用旋转椭球体作为地下空穴模型，通过任意改变长、短轴的长度可以得到不同的旋转椭球体以匹配真实存在的空穴，与常采用的棱柱体和球体相比更具有结构灵活性。文献[12]将旋转椭球体分为绕短轴旋转型和绕长轴旋转型两种，并分别给出详细的重力异常、重力梯度异常的计算表达式（中心位于(0,0, h)(h＜0)处），可直接用于地下空穴产生的重力异常和重力梯度异常信号S(x,y,z)的计算。

Poission指出磁场与重力场具有相似结构。当磁力位和引力位由同一地质体产生，地质体具有均匀的磁化强度和密度且二者之比为常数时，磁位与重力位之间满足泊松关系式[13]。Abt等在局部小范围区域采用相对重力仪进行了重力梯度初步测量，利用泊松关系得到模拟磁力图，与质子旋进磁力仪实测的磁力图进行比较显示二者具有明显相关性[14]。这是因为在局部小范围区域内容易满足泊松关系式所需的两个条件。地下空穴探测往往在浅地表局部小范围区域内进行，可假设满足泊松关系式成立的两个条件，因此，在文献[12]给出的旋转椭球体重力梯度异常计算表达式的基础上，在不考虑剩磁影响的条件下，通过泊松关系式可以推导出旋转椭球体在中心正上方的磁异常计算表达式。

绕短轴旋转型旋转椭球体磁异常计算表达式：

绕长轴旋转型旋转椭球体磁异常计算表达式：

式(2)～(9)中：Hax、Hay和Za分别表示磁异常在X、Y、Z各方向上的分量；ΔT表示总磁异常；μ0为真空中的磁导率，m为磁矩（磁化强度与椭球体体积之积）；i为磁化倾角；δ为X轴正方向与磁化方向的水平分量的夹角；I0为地磁场倾角；A′ 为X轴正方向与地磁场北方向的夹角；式中的其他参数满足的条件见文献[12]。

3 重磁背景场信号的模拟

在整个观测值中，重磁背景场信号B(x,y,z)占最主要的部分。对重磁背景场的模拟基于扰动位T的协方差模型：

式中：dx、dy为平面两点X和Y坐标的差值；z和z′表示两点的垂直高度，由于忽略各测量点的垂直坐标；r表示平面上两点距离；αj和 σj为经验参数。对式(10)进行傅立叶变换，可以得到频域中扰动位T对应的功率谱密度函数：

式中：μ1、μ2为周期频率，μ1=x/(ΔxM)，μ2=y/(ΔyN)，Δx、Δy和M、N分别表示X方向和Y方向的网格间距和测点个数，且x=0,…,M-1，y=0,…,N-1。扰动位的谱T0可由下式计算得到：

为产生随机重磁背景场，引入随机数组变量Λ和Θ，其组成元素为零均值、正态分布且标准差分别为σΛ、σΘ的随机数。显然，在空间域中T0、Λ和Θ的维数与观测网格维数相同，都是 N×M维。重力和重力梯度的谱分别等于扰动位的谱的一阶和二阶导数，故可得重力和重力垂直梯度的谱为：

对各个谱分量求傅立叶逆变换，可以得到空间域中重力gx、gy、gz和重力梯度Γxx、Γyy、Γzz、Γxy、Γyz、Γxz分量，再利用泊松关系式得到空域中磁场背景场，从而得到不同数据类型的背景场信息。

4 仪器或其他随机噪声信号的模拟

用零均值高斯白噪声信号来模拟仪器或其他随机噪声信号N(x,y,z)，不同观测量类型取不同的标准差，如表1所示。

表1 仪器噪声信号的标准差Tab.1 Standard deviations for instrument noise

5 数值仿真

假设测量区域为零高程水平面上 N×M的方形区域。取测点个数N=M=101，间距1 m。测区下方存在三个旋转椭球体结构的空穴，具体参数如表2。

空穴相对周围岩石的密度差Δρ= -2.67 kg/m3，磁化强度差为ΔM= -28 nT，真空中磁导率μ0=4π×10-7N/A2，万有引力常数G=6.673×10-11m3/(kg·s2)。假设磁异常源在地磁场中被磁化，即磁化方向与地磁场方向一致，磁倾角i=I0=68°，磁偏角δ=A′= -7º（具体到特定日期的角度值可查询美国国家海洋和大气管理局网站）。

表2 旋转椭球体结构参数Tab.2 Structural parameters of rotational ellipsoids

将上述参数带入旋转椭球体重力异常、重力梯度异常及磁异常计算表达式，得到图2所示异常信号分量。

图2 地下空穴产生的异常信号Fig.2 Anomaly signals generated by underground voids

取随机数组变量Λ和Θ的标准差σΛ=σΘ=12。经验参数αj和σj取值参照文献[15]附录A。相关参数代入式(10)～(19)求取重力和重力梯度各分量的谱，傅立叶逆变换后得到空间域中重力和重力梯度背景场值，代入泊松关系式得到磁力背景场值，如图3所示。

图3 模拟背景场Fig.3 Simulated background field

将地下空穴产生的重磁异常信号 S(x,y,z)与背景场信号B(x,y,z)及仪器噪声信号N(x,y,z)叠加，得到总重磁观测模拟值Z(x,y,z)，如图4所示。可以看出，各测量数据在合理区间范围，空穴产生的影响在图中有所体现。

取空穴产生的重磁异常信号、模拟背景场信号、随机仪器噪声信号以及合成观测值信号中 y=30 m处的一行数据汇总在同一曲线图中，如图5所示，可更加直观地看出各部分数据所占比重。显然，模拟背景场所占数据比重最高，这与事实相符。另外，由于扰动位的谱与随机数组变量Λ和Θ有关，因此不同的随机数组变量会得到不同的随机背景场，产生不同的观测模拟值，有利于形成多个重磁背景场数据来验证数据解释和反演算法的正确性。

图4 总模拟观测值Fig.4 Total simulated observations

图5 模拟观测值各分量剖面曲线Fig.5 Profiles of simulated observation components

6 结 论

在地下质量异常探测中，科学、合理的重磁观测数值模拟方法是验证数据解释和反演算法正确性的重要前提。本文提出采用更具有结构灵活性的旋转椭球体作为地下空穴模型，在重力异常和重力梯度异常计算表达式的基础上，进一步推导出磁异常计算表达式，实现空穴重磁异常信号的计算。借助扰动位的功率谱密度函数实现频域中重力异常谱和重力梯度异常谱的计算，实现背景场的统计学模拟。通过添加高斯白噪声信号模拟测量产生的仪器或其他随机噪声信号。

数值仿真给出了实现效果，显示出该方法的合理性和有效性，为研究天然地下空穴的存在状态和物性参数提供了保障。

（References）：

[1] Ardestani V E. Detecting, delineating and modeling the connected solution cavities in a dam site via microgravity data[J]. Acta Geod Geophys, 2013, 48: 123-138.

[2] Mochales T, Casas A M. Detection of underground cavities by combining gravity, magnetic and ground penetrating radar surveys: A case study from the Zaragoza area, NE Spain[J]. Environ Geol, 2008, 53: 1067-1077.

[3] Leucci G, Giorgi L D. Microgravimetric and ground penetrating radar geophysical methods to map the shallow karstic cavities network in a coastal area (Marina Di Capilungo, Lecce, Italy)[J]. Exploration Geophysics, 2010, 41: 178-188.

[4] Dumrongchai P. Small anomalous mass detection from airborne gradiometry[D]. Columbus, Ohio: The Ohio State University, 2007.

[5] Styles P, Toon S, Thomas E, et al. Microgravity as a tool for the detection, characterization and prediction of geohazard posed by abandoned mining cavities[J]. First Break, 2006, 24(5): 51-60.

[6] Moritz H. Advanced physical geodesy[M]. Karlsruhe: Abacus Press, 1980: 56-59.

[7] Jekeli C. The statistics of the Earth’s gravity field, Revisited[J]. Manuscripta Geodaetica, 1991, 16: 313-325.

[8] Abt T, Huang O, Jekeli C. Combination of local gravimetry and magnetic data to detect subsurface anomalies using a matched filter[C]//International Association of Geodesy Symposia 136. Berlin Heidelberg: Springer-Verlag, 2012: 247-254.

[9] Cardarelli E, Marrone C, Orlando L. Evaluation of tunnel stability using integrated geophysical methods[J]. Journal of Applied Geophysics, 2003, 52(2): 93-102.

[10] Jekeli C. Statistical analysis of Moving-Base gravimetry and gravity gradiometry. OSU Report 466[R]. Columbus, Ohio: The Ohio State University, 2003.

[11] Jekeli C, Abt T L. The statistical performance of the matched filter for anomaly detection using gravity gradients. OSU Report 494[R]. Columbus, Ohio: The Ohio State University, 2010.

[12] 刘繁明, 刘惠敏, 荆心. 基于旋转椭球体的等轴状场源体重力异常场模型构建方法[J]. 中国惯性技术学报, 2016, 24(1): 26-29. Liu Fan-ming, Liu Hui-min, Jing Xin. Modeling of isometric field-source in gravity anomaly field based on rotational ellipsoid[J]. Journal of Chinese Inertial Technology, 2016, 24(1): 26-29.

[13] Jekeli C, Abt T L. Locating anomalies by gravity gradiometry using the matched filter with a probabilistic assessment[J]. Studia Geophysica et Geodaetica, 2011: 55(1): 1-19.

[14] Dumrongchai P. Small anomalous mass detection from airborne gradiometry. OSU Report 482[R]. Columbus, Ohio: The Ohio State University, 2007.

[15] Abt T L. Detection of a local mass anomaly in the shallow subsurface by applying a matched filter. OSU Report 497[R]. Columbus, Ohio: The Ohio State University, 2011.

Numerical simulation method for gravity and magnetic observations of mass anomaly in shallow subsurface

LIU Hui-min1,2, LIU Fan-ming1, JING Xin1, XIA Lin-lin3
(1. College of Automation, Harbin Engineering University, Harbin 150001, China; 2. College of Mechanical and Electrical Engineering, Qingdao Agricultural University, Qingdao 266109, China; 3. School of Automation Engineering, Northeast Dianli University, Jilin 132012, China)

A numerical simulation method for gravity and magnetic observations of mass anomaly in shallow subsurface is proposed for studying the data interpretation and inversion algorithm of underground voids, assuming that the gravity and magnetic anomaly signals generated by the underground voids can be described by the signals of rotational ellipsoids. One important aspect of the study is to derive the magnetic anomaly expressions of rotational ellipsoids with the known formulas of gravity and gravity gradient. The power spectral density function of disturbing potential is used to calculate the gravity anomaly spectrum and gravity gradient anomaly spectrum. The inverse Fourier transformation of the above spectra provides the gravity and magnetic background filed signals in the space domain. Furthermore, the random instrument noise signals are simulated by Gaussian white noise. Superposing these three parts of data and then carrying out numerical simulation for them with specific parameters, the results show that the method is reasonable and effective.

gravity anomaly; gravity gradient anomaly; magnetic anomaly; underground voids; gravity and magnetic background field; numerical simulation

P312.1

：A

2016-04-14；

：2016-07-29

国家自然科学基金（60834005）；国家自然科学基金青年基金项目（61503073）

刘惠敏（1983—），女，博士研究生，从事微弱信号测量与处理技术方向研究。E-mail: lhmgct@126.com

联 系 人：刘繁明（1963—），男，教授，博士生导师。E-mail: Hrblfm407@126.com

1005-6734(2016)04-0490-06

10.13695/j.cnki.12-1222/o3.2016.04.013