2020, Vol. 41
2. 东北大学 深部金属矿山安全开采教育部重点实验室，辽宁 沈阳 110819
2. Key Laboratory of Ministry of Education on Safe Mining of Deep Metal Mines, Northeastern University, Shenyang 110819
与传统重力异常数据相比，重力梯度数据包含有更多的频率信息，能够更好地用于解译地下矿产的分布情况.利用重力梯度数据可以进行地下目标体的成像与反演，相关成像及其反演是常用的方法之一.该成像方法最初被应用于自然电场数据[1]，是一种利用观测异常与几何函数间的相关系数描述地质体空间位置和形态的算法，能够应用于不同类型的地球物理场数据，包括重力异常[2]、重力梯度数据[3]、磁异常[4]以及电阻率异常[5]等；在成像基础上加入物性参数的变化量能够反演地下目标体的物性分布，国内学者在该方向开展了较深入的研究并取得了良好的效果[6-8].然而，相关成像反演在深度方向上的分辨能力较低；根据其他反演方法的研究，引入深度加权函数虽可能改善结果的纵向分辨率，但当地下存在多个目标体的复杂情况时，结果仍与真实分布有一定差距.
因此，本文将联合多个重力梯度数据，使用深度加权函数与分区加权的方法提高相关成像反演的空间分辨率.在数据试验中，可证实提出的方法大幅提高了纵向分辨能力，且具有一定的可行性.
1 重力梯度数据联合相关成像反演原理设研究区域位于笛卡尔坐标系中，o为坐标原点，xoy为一水平面，z轴的正方向垂直向下；地下空间被划分为一系列点元，第j个点元在第i个观测点上引起的重力梯度Vαβ, j为
|
(1) |
其中:γ表示万有引力常数；ρj，νj分别为第j个点元的密度、体积；Bαβ, j表示第j个质点在第i个观测点上引起的重力梯度数据的几何函数，则第j个点元的重力梯度相关系数为
|
(2) |
由于相关成像只能计算相关系数，不能反映地下的物性分布.因此，使用相关成像反演方法得到地下物性参数分布情况.首先，建立地下密度分布密度的初始模型，设定密度变化量Δρ，在初始模型的基础上，利用相关系数与密度变化量的乘积，进行迭代计算，公式为
|
(3) |
其中:iter表示迭代的次数; ρj, iter表示第iter次迭代过程中第j个点元的密度值; Cj的意义同公式(2)；Δρ表示密度变化量，根据Hou等的研究[7]，密度变化量可利用观测数据的最大值除以G×C的最大值来设定:
|
(4) |
其中，G表示正演矩阵.
通过联合不同的梯度分量，能够更加全面地反映地下异常体信息，降低反演问题的多解性，增加结果的可靠性.根据Hou等的研究[7]，多分量重力梯度数据联合反演中的相关系数为
|
(5) |
当进行联合反演时，Δρ存在多个不同的值，通常选取最小值参与反演.
如前文提到的，相关成像反演在深度方向上的分辨率有待提高，特别是对目标体底部位置与形态识别的精度较差.下面引入深度加权函数来提高反演的纵向分辨率，本文采用Commer等提出的空间梯度加权函数[9-10]，该函数能够引入先验深度信息改善反演结果，其形式为
|
(6) |
其中:z表示点元中心深度；dz表示研究区域内地下深度范围；α为经验常数，一般为较小的正数；r表示缩放因子；zc1, zc2分别表示模型顶、底部深度, 根据先验地质信息取值.深度加权后的迭代反演公式为
|
(7) |
上述加权方法对地下空间存在孤立目标体时效果较好，如果地下不同深度上存在多个目标体，不能直接通过公式(7)加权，否则将导致所有异常体处于同一深度范围内.因此，在处理该情况时，首先根据先验信息划分研究区域，为不同区域分别设定该区内目标体的顶、底埋深，再应用公式(7)进行反演.
2 模型试验为验证提出方法的有效性，设计立方体组合模型进行数据试验.模型位于笛卡尔坐标系下，研究区域设定在坐标轴x，y方向范围均为0~2 000 m、z方向范围为0~1 000 m；等分地下空间为20×20×20个立方体，小立方体大小为100 m×100 m×50 m.组合模型由3个立方体A，B和C组成(图 1)，A立方体大小为200 m×300 m×300 m，剩余密度为0.5 g/cm3；B立方体大小为200 m×200 m×200 m，剩余密度为0.7 g/cm3；C立方体大小为500 m×500 m×200 m，剩余密度为0.5 g/cm3；A和B立方体的顶面埋深均为200 m，C立方体的顶面埋深为700 m.为组合模型产生的重力梯度数据加入5%的高斯噪声验证算法的抗噪能力；在试验中使用x=1 500 m，y=1 000 m，z=300 m三个剖面进行结果分析.
|
图 1 模型空间位置示意图 Fig.1 Spatial location of the model |
根据秦朋波的研究[11]，对比分析不同重力梯度分量的特点，可得出Vxz|Vyz|Vzz组合反演效果最好.因此，本文联合这三个梯度分量开展试验，并设定最大迭代次数为100次.为对比深度加权改善的效果，首先使用未改进的相关成像反演(见图 2).经16次迭代计算后，结果的最大值达到0.7 g/cm3，与理论模型一致.在图 2c中，可见反演的横向分辨率较高，水平方向上的形态和模型基本一致.但是由图 2a和图 2b中识别出的结果可见，深度方向上的分辨率较低，模型B的下延深度过大，而模型C几乎无法辨识出.
|
图 2 使用未改进的相关成像反演结果(矩形框为模型位置) Fig.2 Results from unimproved correlation imaging Inversion (rectangular box indicates the position of model) (a)—x=1 500 m; (b)—y=1 000 m; (c)—z=300 m. |
下面使用基于空间分区深度加权的方法进行反演.首先对地下空间进行分区，使分区后每个立方体独自占据一个区域，进而使用深度加权函数为立方体A，B，C分别赋值.立方体A，B的zc1均取200 m，由于下底面深度不同，则模型A的zc2=500 m，模型B的zc2=400 m；对于立方体C，zc1=700 m，zc2=900 m.加权函数中α=0.001，dz=1 000 m，r=20.在16次迭代后，反演结果见图 3，剩余密度最大值达到0.7 g/cm3，与理论模型保持一致.分别对比图 3a和图 2a，图 3b和图 2b，可见本文提出的算法能够较准确地计算出下底深度，并能同时计算出不同深度的目标体的分布范围，反演的纵向分辨率和未改进的算法相比大幅提升；结果的数值范围也与理论模型基本一致.此外，结果还反映出本文方法具有良好的抗噪能力，可用于模拟实测数据的反演研究.
|
图 3 使用基于深度加权的相关成像反演结果(矩形框为模型位置) Fig.3 Results from correlation imaging inversion based on depth weighting (rectangular box indicates the position of model) (a)—x=1 500 m; (b)—y=1 000 m; (c)—z=300 m. |
将本文方法应用于美国文顿盐丘地区的航空实测重力梯度数据以验证其可行性.许多学者都曾对文顿盐丘地区开展了大量研究[12-16].该地区地下存在一剩余密度约为0.55 g/cm3的盖岩，采集的异常数据便是由它产生的.根据Ennen等[13]、Geng等[14]的研究，盖岩的深度范围约在160~650 m.令数据位于WGS84坐标系下，x和y方向范围分别为440.5~444.5 km, 3 332.78~3 336.78 km，研究区域内均匀分布20×20个测点.设最大反演深度1 km，等分地下空间为20层.在反演前已对数据进行地形校正和高通滤波处理.试验将选取x=442.56 km，y=3 334.44 km和z=375 m三个剖面进行分析.
根据模型试验结论，仍使用Vxz|Vyz|Vzz组合进行联合相关成像反演.在深度加权中，设α=0.001，dz=1 000 m，r=20；结合文顿盐丘地区的先验信息[12], 设盖岩的顶部深度zc1=200 m，底部深度zc2=500 m；根据盖岩已知范围，将研究区域分成“回”字形的两部分，盖岩模型占据中间区域.经过17次迭代得到反演结果(见图 4)，其中最大剩余密度达到0.5 g/cm3，接近盖岩的已知密度信息；从图 4a和图 4b可见，盖岩的深度范围在200~500 m，其形态也与其他学者研究的结果一致[12-16].因此，认为本文提出的方法能够用于实测数据的反演解释.
|
图 4 使用基于深度加权的实测数据反演结果 Fig.4 Results from the inversion based on depth weighting with real data (a)—x=442.56 km; (b)—y=3 334.44 km; (c)—z=375 m. |
本文基于深度加权技术，提出了多分量重力梯度数据联合相关成像反演方法.在数据试验中，联合Vxz, Vyz和Vzz三个梯度分量实现反演，结果较为准确地反映出目标体的位置与形态，密度数值范围和理论值基本一致；对比未改进的方法，深度分辨率大幅提高；同时，验证了算法具有抗噪性.在文顿盐丘地区的实测数据试验中，反演出的盖岩形态与密度分布均和已知资料一致.上述试验证明了方法的有效性和可行性.
致谢 感谢Bell Geospace公司提供的实测重力梯度数据.
| [1] |
Patella D. Introduction to ground surface self-potential tomography[J]. Geophysical Prospecting, 1997, 45(4): 653-681. DOI:10.1046/j.1365-2478.1997.430277.x |
| [2] |
郭良辉, 孟小红, 石磊, 等. 重力和重力梯度数据三维相关成像[J]. 地球物理学报, 2009, 52(4): 1098-1106. (Guo Liang-hui, Meng Xiao-hong, Shi Lei, et al. 3D correlation imaging for gravity and gravity gradiometry data[J]. Chinese Journal of Geophysics, 2009, 52(4): 1098-1106.) |
| [3] |
Guo L H, Meng X H, Shi L. 3D correlation imaging of the vertical gradient of gravity data[J]. Journal of Geophysics and Engineering, 2011, 8(1): 6-12. DOI:10.1088/1742-2132/8/1/002 |
| [4] |
郭良辉, 孟小红, 石磊. 磁异常ΔT三维相关成像[J]. 地球物理学报, 2010, 53(2): 435-441. (Guo Liang-hui, Meng Xiao-hong, Shi Lei. 3D correlation imaging for magnetic anomaly △T data[J]. Chinese Journal of Geophysics, 2010, 53(2): 435-441.) |
| [5] |
Mauriello P, Patella D. Resistivity anomaly imaging by probability tomography[J]. Geophysical Prospecting, 1999, 47(3): 411-429. DOI:10.1046/j.1365-2478.1999.00137.x |
| [6] |
Liu G F, Yan H F, Meng X H, et al. An extension of gravity probability tomography imaging[J]. Journal of Applied Geophysics, 2014, 102: 62-67. DOI:10.1016/j.jappgeo.2013.12.012 |
| [7] |
Hou Z, Huang D, Wei X. Fast inversion of probability tomography with gravity gradiometry data based on hybrid parallel programming[J]. Journal of Applied Geophysics, 2016, 124: 27-38. DOI:10.1016/j.jappgeo.2015.11.009 |
| [8] |
Hou Z, Huang D, Wang E, et al. 3D density inversion of gravity gradiometry data with a multilevel hybrid parallel algorithm[J]. Applied Geophysics, 2019, 16(2): 141-152. DOI:10.1007/s11770-019-0763-4 |
| [9] |
Commer M. Three-dimensional gravity modeling and focusing inversion using rectangular meshes[J]. Geophysical Prospecting, 2011, 59(5): 966-979. |
| [10] |
Commer M, Newman G A, Williams K H, et al. 3D induced-polarization data inversion for complex resistivity[J]. Geophysics, 2011, 76(3): 157-171. |
| [11] |
秦朋波.重力和重力梯度数据联合反演方法研究[D].长春: 吉林大学, 2016. (Qin Peng-bo.A study on integrated gravity and gravity gradient data in inversion[D].Changchun: Jilin University, 2016. ) |
| [12] |
Coker M O, Bhattacharya J P, Marfurt K J. Fracture patterns within mudstones on the flanks of a salt dome:syneresis or slumping?[J]. Gulf Coast Association of Geological Societies, 2007, 57: 125-137. |
| [13] |
Ennen C, Hall S.Structural mapping of the Vinton salt dome, Louisiana, using gravity gradiometry data[C]//SEG Annual Meeting.San Antonio, 2011: 830-835.
|
| [14] |
Geng M, Huang D, Yang Q, et al. 3D inversion of airborne gravity gradiometry data using cokriging[J]. Geophysics, 2014, 79(4): G37-G47. DOI:10.1190/geo2013-0393.1 |
| [15] |
Thompson S A, Eichelberger O H. Vinton salt dome, Calcasieu Parish, Louisiana[J]. AAPG Bulletin, 1928, 12(4): 385-394. |
| [16] |
Oliveira V C, Jr Barbosa V C F. 3-D radial gravity gradient inversion[J]. Geophysical Journal International, 2013, 195(2): 883-902. DOI:10.1093/gji/ggt307 |