-
相对于传统的重力测量, 重力梯度数据只有很窄的带宽, 浅层和中层地质特征所产生的高频信号在重力梯度的带宽中有着更显著的体现[1], 在探测浅层异常体中具有更高的分辨能力和更精确的灵敏度,在区域地质调查和矿藏探测等方面都发挥着重要的作用。重力梯度丰富的张量数据可用于探测异常体的边缘、拐点和中心点,是提高异常体解释精度的强有力工具[2]。作为一种前沿性的重力测量方法,许多学者对重力梯度数据的后期处理也进行了大量的研究。Zhang等[3]利用张量欧拉反褶积法解释了重力梯度全张量数据,在知道场源体的前提下,利用重力梯度的6个分量,得到更加明确紧凑的地质结构的解;Zhdanov等[4]提出了一种将最小梯度支撑函数作为前提条件,使用重加权正则化共轭梯度的反演算法,其反演结果具有较好的边界特征;Martinez等[5]给出了基于平滑约束的重力梯度张量数据的反演;蒋甫玉等[6]通过三度体球冠模型对重力梯度张量进行了研究分析,为寻找储油构造提供了可靠的理论基础;孙中科等[7]使用拟BP神经网络算法以及粒子群算法对重力梯度张量进行了反演。随着卫星及航空重力梯度测量系统投入使用,越来越多的实测重力梯度数据需要处理、分析和解释。但已有的研究成果多集中于使用算法对设计的正演模型进行数值模拟,当算法应用于实测数据时,由于实际情况的复杂性,会导致反演结果的不确定性。
本文基于以上问题,对预条件共轭梯度反演算法加以改进。当算法仅应用于正演模型计算时,拉格朗日经验参数值变化较小,但对于实际复杂情况,参数值不易确定,故使用L曲线的拐点值代替原有反演算法中的拉格朗日经验参数作为正则化参数;当研究的区域较大且观测数据量丰富时,反演中存在病态性问题且核函数快速衰减,将地下模型改进为不等间隔模型; 设计更接近实际复杂情况的真实模型进行数值模拟,充分利用重力梯度的5个分量,使用预条件共轭梯度算法进行联合反演;将改进后的算法应用于澳大利亚西澳洲Kauring地区实测重力梯度场,给出了测区内密度异常体的分布信息。
-
首先确定观测数据与模型参数之间的函数关系,构建目标函数;为使改进后的算法有效应用于实测数据的反演,用L曲线的拐点值代替拉格朗日经验参数;为对模型空间结构加以约束,引入粗糙度矩阵,从而改善反演中存在的病态性问题;引入深度加权函数,使核函数可以反映出各个深度真实的权值;最终利用改进后的预条件共轭梯度法作为反演算法,高效地完成反演计算。
-
反演中,目标函数φ由数据拟合差函数φd和模型目标函数φm构成,解决反演问题即在数据拟合差函数满足条件的前提下,寻找模型向量m使求得的目标函数φ值最小[1]:
$$ \left\{ {\begin{array}{*{20}{c}} {\min :\varphi = {\varphi _d} + \mu {\varphi _m}\;\;\;\;\;\;\;\;\;\;\;\;\;\left( 1 \right)}\\ {{\varphi _d} = \sum\limits_{i = 1}^N {{{\left( {\frac{{\Delta {d_i} - {G_i}\Delta \mathit{\boldsymbol{m}}}}{{{\sigma _i}}}} \right)}^2} = \;\;\;\;\;\;\;\;\;} }\\ {{{\left\| {{\mathit{\boldsymbol{W}}_d}\left( {\Delta \mathit{\boldsymbol{d}} - \mathit{\boldsymbol{G}}\Delta \mathit{\boldsymbol{m}}} \right)} \right\|}^2}\;\;\;\;\;\;\;\;\;\;\;\left( 2 \right)} \end{array}} \right. $$ 式中,Δm=m-m0代表模型参数向量m与初始模型m0之间的改正量;G代表核函数,即表示从模型单元到观测点的线性投影算子;Δd代表相对应的观测值的改正量;Wd为一个对角矩阵; σi为第i个数据的标准偏差。
-
预条件共轭梯度反演方法中,使用拉格朗日乘子作为正则化参数。当求解大规模矩阵时,对拉格朗日乘子的求解需要大量运算,采用经验值作为正则化参数的相对最佳取值。由于正则化参数起到平衡数据拟合泛函与模型拟合泛函的作用,取值过大,会导致模型响应与观测值相差较大; 取值过小,会导致模型与实际情况不符, 使反演结果不可靠。
本文使用L曲线的曲率拐点值作为正则化参数。L曲线是基于实际数据拟合与模型目标函数之间进行权衡的一种准则,适合大规模问题的求解。当把模型目标函数与数据拟合函数画在对数-对数图上,所呈现的曲线有一个明显的拐点; 当正则化参数值大于拐点时, 模型目标函数变化不大, 数据拟合差函数变化很大;当正则化参数值小于拐点时,模型目标函数快速增长,而数据拟合差函数呈缓慢递减趋势[8]。所以将拐点作为最佳的正则化参数,L曲线的曲率k定义为:
$$ k = \frac{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \rho } '\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \eta } '' - \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \rho } ''\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \eta } '}}{{{{\left[ {{{\left( {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \rho } '} \right)}^2} + {{\left( {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \eta } '} \right)}^2}} \right]}^{3/2}}}} $$ (3) 式中,$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \rho } = \ln {\varphi _d};\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \eta } = \ln {\varphi _m}; $符号“′”和“″”分别表示函数的一阶和二阶导数。因此反演中会根据实际数据所构建的模型寻求L曲线的最大曲率。
-
为了使反演图像能够在3个方向上连续变化,光滑过渡,将粗糙度矩阵引入到重力梯度张量的三维成像反演中[9]。定义三维模型的粗糙度矩阵R为模型向量m在x、y、z 3个方向上一阶偏微分的平方和:
$$ \begin{array}{l} \mathit{\boldsymbol{R}} = {\left\| {{\partial _x}\Delta \mathit{\boldsymbol{m}}} \right\|^2} + {\left\| {{\partial _y}\Delta \mathit{\boldsymbol{m}}} \right\|^2} + {\left\| {{\partial _z}\Delta \mathit{\boldsymbol{m}}} \right\|^2} = \\ \int {{{\left( {\frac{{\partial \Delta \mathit{\boldsymbol{m}}}}{{\partial x}}} \right)}^2}{\rm{d}}v} + \int {{{\left( {\frac{{\partial \Delta \mathit{\boldsymbol{m}}}}{{\partial y}}} \right)}^2}{\rm{d}}v} + \int {{{\left( {\frac{{\partial \Delta \mathit{\boldsymbol{m}}}}{{\partial z}}} \right)}^2}{\rm{d}}v} \end{array} $$ (4) $$ \mathit{\boldsymbol{R}} = \Delta \mathit{\boldsymbol{m}}\left( {\mathit{\boldsymbol{R}}_x^{\rm{T}}{\mathit{\boldsymbol{R}}_x} + \mathit{\boldsymbol{R}}_y^{\rm{T}}{\mathit{\boldsymbol{R}}_y} + \mathit{\boldsymbol{R}}_z^{\rm{T}}{\mathit{\boldsymbol{R}}_Z}} \right)\Delta \mathit{\boldsymbol{m}} $$ (5) 式中,v为模型向量所在的体积空间。
-
由于重力梯度值与场源到观测点距离的立方成反比,重力梯度数据没有固有的深度分辨率,所构建的核函数矩阵为线性的;随着深度的增加,核函数快速衰减,使反演结果集中于地表附近,难以反映出深部异常体的真实位置。通过引入深度加权因子,可以抵消核函数的快速衰减,采用文献[10]方法引入深度加权函数如下:
$$ W\left( z \right) = \frac{1}{{{{\left( {Z + {Z_0}} \right)}^{\beta /2}}}} $$ (6) 式中,Z0取决于观测面的高度以及单元体的尺寸大小;Z为每层网格单元中心点的埋深;β为常数,对于重力梯度数据β=3,函数图像与重力梯度核函数图像衰减一致。
在网格单元划分时,浅层采用高分辨率格网,深层为较低分辨率格网,可以有效增加深层块体在地表处的重力梯度效应,减弱核函数的快速衰减,同时也可以减少反演模型的参数个数,改善反演中存在的病态性问题。
-
将粗糙度矩阵以及深度加权函数代入到模型目标函数中[11]:
$$ \begin{array}{l} {\varphi _m}\left( \mathit{\boldsymbol{m}} \right) = {\alpha _s}\int {_V{{\left( {\partial \mathit{\boldsymbol{W}}\left( z \right)\Delta \mathit{\boldsymbol{m}}} \right)}^2}{\rm{d}}v} + {\alpha _x}\int {_V{{\left( {\frac{{\partial \mathit{\boldsymbol{W}}\left( z \right)\Delta \mathit{\boldsymbol{m}}}}{{\partial x}}} \right)}^2}{\rm{d}}v} + \\ \;\;\;\;{\alpha _y}\int {_V{{\left( {\frac{{\partial \mathit{\boldsymbol{W}}\left( z \right)\Delta \mathit{\boldsymbol{m}}}}{{\partial y}}} \right)}^2}{\rm{d}}v} + {\alpha _z}\int {_V{{\left( {\frac{{\partial \mathit{\boldsymbol{W}}\left( z \right)\Delta \mathit{\boldsymbol{m}}}}{{\partial z}}} \right)}^2}{\rm{d}}v} \end{array} $$ (7) 式中,$ {\alpha _s}、{\alpha _x}、{\alpha _y}、{\alpha _z}$为模型目标函数中各项的加权因子。
通过利用有限差分法代替微分,将式(7)改写成矩阵形式:
$$ \begin{array}{l} {\varphi _m}\left( \mathit{\boldsymbol{m}} \right) = \Delta {\mathit{\boldsymbol{m}}^{\rm{T}}}\left( {\mathit{\boldsymbol{W}}_s^{\rm{T}}{\mathit{\boldsymbol{W}}_s} + \mathit{\boldsymbol{W}}_x^{\rm{T}}{\mathit{\boldsymbol{W}}_x}{\rm{ + }}\mathit{\boldsymbol{W}}_y^{\rm{T}}{\mathit{\boldsymbol{W}}_y}{\rm{ + }}} \right.\\ \;\;\;\;\;\;\left. {\mathit{\boldsymbol{W}}_z^{\rm{T}}{\mathit{\boldsymbol{W}}_z}} \right)\Delta \mathit{\boldsymbol{m}} = \Delta {\mathit{\boldsymbol{m}}^{\rm{T}}}\mathit{\boldsymbol{W}}_i^{\rm{T}}{\mathit{\boldsymbol{W}}_i}\Delta \mathit{\boldsymbol{m}} \end{array} $$ (8) $$ {\mathit{\boldsymbol{W}}_i} = {\alpha _i}{\mathit{\boldsymbol{R}}_i}\mathit{\boldsymbol{D}}, i = s, x, y, z $$ (9) 式中,αi为目标函数中各项的权重系数;Ri为每一个分量的差分算子;D代表离散化的深度加权函数矩阵;s、x、y、z分别表示式(7)中加权因子参与计算的4个分量。将构建的模型函数代入到目标函数中:
$$ \begin{array}{l} \varphi = {\left( {\mathit{\boldsymbol{d}} - \mathit{\boldsymbol{G}}\Delta \mathit{\boldsymbol{m}}} \right)^{\rm{T}}}\left( {\mathit{\boldsymbol{d}} - \mathit{\boldsymbol{G}}\Delta \mathit{\boldsymbol{m}}} \right) + \\ \;\;\;\;\;\;\mu \Delta {\mathit{\boldsymbol{m}}^{\rm{T}}}\mathit{\boldsymbol{W}}_i^{\rm{T}}{\mathit{\boldsymbol{W}}_i}\Delta \mathit{\boldsymbol{m}} \end{array} $$ (10) 改写成方程组形式:
$$ \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{G}}\\ {\sqrt \mu {\mathit{\boldsymbol{W}}_i}} \end{array}} \right]\Delta \mathit{\boldsymbol{m}} = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{d}}\\ {\bf{0}} \end{array}} \right] $$ (11) 用$\mathit{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{G}}\\ {\sqrt \mu {\mathit{\boldsymbol{W}}_i}} \end{array}} \right]$代表雅克比矩阵,$\mathit{\boldsymbol{b}} = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{d}}\\ {\bf{0}} \end{array}} \right]$,则式(11)写成:
$$ \mathit{\boldsymbol{A}}\Delta \mathit{\boldsymbol{m = b}} $$ (12) 位场反演中存在病态性与非唯一性,且经常需要对大规模问题求解,而共轭梯度法(conjugate gradient algorithm, CGA)只需利用函数的一阶导数信息,存储量小,算法简便,收敛速度比最速下降法快,适合求解大规模问题。CGA主要用于求解二次方程$\varphi = \frac{1}{2}{\mathit{\boldsymbol{m}}^{\rm{T}}}\mathit{\boldsymbol{Am}} - {\mathit{\boldsymbol{m}}^{\rm{T}}}\mathit{\boldsymbol{b}}$,解为$\mathit{\boldsymbol{m = }}{\mathit{\boldsymbol{A}}^{ - 1}}\mathit{\boldsymbol{b}}$,与方程(12)等价。该算法利用初始点的梯度作为最初的共轭方向,沿着这组方向而非负梯度方向去搜索目标函数的极小点,改善解的稳定性[12]。反演过程如下:
1) 已知观测值dobs和初始模型m0k, 其中k=1, 2 …Nmax, Nmax为外循环最大迭代次数。
2) 计算观测数据和由初始模型计算得到的理论数据之差Δdk=dobs-Gm0k。
3) 设定初始变量: $\mathit{\boldsymbol{A}} = {\left[ {\mathit{\boldsymbol{G}}\sqrt {{\mu ^{ - 1}}} {\mathit{\boldsymbol{W}}_i}} \right]^{\rm{T}}}, {\mathit{\boldsymbol{b}}^k} = \left[ {\Delta {\mathit{\boldsymbol{d}}^k}\;\;0} \right], \Delta \mathit{\boldsymbol{m}}_0^k = {\bf{0}}, {\mathit{\boldsymbol{r}}_0} = {\mathit{\boldsymbol{A}}^{\rm{T}}}\left( {\mathit{\boldsymbol{b}} - \mathit{\boldsymbol{A}}\Delta \mathit{\boldsymbol{m}}_0^k} \right), {\mathit{\boldsymbol{z}}_0} = \mathit{\boldsymbol{S}}{\mathit{\boldsymbol{r}}_0}, {\mathit{\boldsymbol{p}}_0} = {\mathit{\boldsymbol{z}}_0}, \mathit{\boldsymbol{S}}$代表预条件矩阵。
4) 设定i=1, 2 …NCGmax, 其中NCGmax为内循环最大迭代次数,进行预条件共轭梯度算法计算,其中qi=Api, αi=(ri-1Tzi-1)/qiTqi, Δmik=Δmi-1k+αipi, ri=ri-1-αiATqi, zi=Sri。zk为改善的搜索方向; βi=(riTzi)/(ri-1Tzi-1), pi=zi-1+βipi-1。
5) 共轭梯度反演内循环计算结束, 设定m1k=m0k+Δmk。
6) 设定m0k=m1k, 重建初始模型并返回步骤2),继续迭代计算,直到满足收敛条件。
预条件矩阵S=(ATA)-1可以改善条件数,有效提高迭代速度,但求解大规模矩阵时对矩阵S的求解较为困难,在实际中用深度加权函数W(z)作为预条件因子,同样可以达到改善条件数的作用[11]。
-
为验证改进后的预条件共轭梯度算法应用于密度结构反演的有效性,利用加入随机噪声的重力梯度响应作为观测值,并将地下模型改进为不等间隔模型,经过单分量反演后,将结果和真实模型进行比较分析。为充分利用重力梯度张量数据,改善反演中解的非唯一性,将5个独立分量进行联合反演,与单分量的反演结果对比,突出联合反演的优势。为表明将反演算法应用于航空重力梯度数据的可靠性,设计的真实模型中包含多个形状大小、密度值和质心埋深各异的异常体,对其进行联合反演并对结果分析。
-
为了测试模型目标函数中依次加入正则化参数、粗糙度矩阵以及深度加权函数对反演结果的影响,根据构建的模型进行反演试算。图 1(a)代表所构建的地下模型是1 000×1 000 m2区域范围内埋藏的单个异常体,记为异常体A。将地下模型划分为10层,使浅层模型单元的中心点与观测点所在的位置对应,深层模型单元的宽度逐渐变大;1~3层模型单元的行数及列数均为25,4~5层为20,6~8层为15,9~10层为10;每层模型单元大小从40 m递增到100 m,共由3 550个模型单元构成。模型尺寸统计信息见表 1。地表观测点的间隔为50 m,总共有400个观测点,观测面高度为80 m,初始模型矩阵m的值为0,将加入5%随机噪声的Tzz作为观测值,对反演结果进行比较分析。
表 1 模型统计信息
Table 1. Statistics of Model Information
密度异常体 尺寸(L×W×H)/m3 密度差/(kg·m-3) 质心平面/m 质心位置/m 红色异常体A 200×200×200 1 000 (500, 500) 200 红色异常体B 200×200×200 800 (-125, 0) 200 蓝色异常体C 200×200×200 -800 (375, 0) 200 图 1(b)仅加入正则化参数,结果显示异常体集中于地表附近。图 1(c)再加入粗糙度矩阵,可以明显改善模型单元体的平滑度,避免相邻单元密度值的突跳现象。图 1(d)在此基础上加入深度加权函数,克服重力场固有的“趋肤”效应,使异常体归置于相应深度。将三者结合后,反演效果最佳,异常体集中,更接近于真实模型。
-
重力勘探中,重力梯度数据在浅层比重力数据具有更高的分辨能力。虽然从理论上同一方法不同分量无法改变由位场性质决定的多解性,但在同一高度不同分量的观测数据可以使数据更加丰富完整,改善反演中的多解性[13]。模型统计信息见表 1红色异常体B和蓝色异常体C。初始模型矩阵m的值为0,在观测值中加入5%的随机噪声。图 2(a)代表所构建的地下模型是1 000×1 000 m2区域范围内埋藏的两个并排的正负异常体,图 2(b)~2(g)列出了重力梯度张量单分量的反演结果。
图 2 单分量反演与多分量联合反演结果及其分析
Figure 2. Inversion of Independent Component and Joint-Inversion of Five Independent Components
由图 2(b)至图 2(g)可知,与真实模型相比,所得异常体的剩余密度值、形状、位置、边界均存在不同程度的偏差,Txx、Tyy、Txz、Tyz、Tzz为5个单独分量。利用Txx和Tyy得到的反演结果在水平方向上互补,但反演得到的剩余密度值偏差较大,振幅可达±0.5 g/cm3;利用Txz和Tyz得到的反演结果集中于异常周围, 可以反映地下密度异常体的边界信息,得到的剩余密度值偏低; Tzz反演结果与真实模型对应较好, 但异常体中心位置偏下,得到的剩余密度值偏低。为充分利用各单分量反演结果的优势并改善反演中解的不唯一性,利用重力梯度的5个分量联合反演的结果如图 2(h)所示,得到的反演结果更精确,异常体在深度方向上更集中于中心,剩余密度值与真实模型密度值更接近,边界明显更容易确定出异常体的形状。为进一步分析反演结果精度,图 2(i)代表包含噪声的观测值Tzz,等值线代表由联合反演得到的剩余密度值计算出的Tzz′,图形吻合较好,差值较小; 计算Tzz和Tzz′的残差值,图 2(j)代表两者之间的残差直方图,残差值的标准差为±0.977 7 E,表明反演结果可靠、有效。
-
为了增强航测重力梯度数据反演结果的可靠性,增加真实模型的复杂性,设计的模型由多个形状大小各异的异常体组成。图 3(a)代表所构建的地下模型是1 000×1 000 m2区域范围内埋藏的多个异常体。初始矩阵m的值为0,在观测值中加入5%随机噪声。利用重力梯度的5个分量进行联合反演,图 3(b)~3(c)表示真实模型在不同深度的横截面;图 3(d)~3(f)表示联合反演结果及其在不同深度的横截面;图 3(g)中的等值线为联合反演结果计算得到的Tzz′;图 3(h)为联合反演结果计算得到的Tzz′与观测值Tzz的残差直方图。
图 3 原始模型与反演结果对比及其分析
Figure 3. Comparison Between Original Model and Inversion Result with Analysis of Inversion Results
由图 3(e)可以看出, 反演结果在浅部反映出异常体的形状和位置更接近真实模型。随着深度的增加,由于梯形异常体由几个异常块体叠加而成,形状较为复杂,由图 3(f)可以看出反演结果的形状、位置和真实模型稍有偏差,小的异常体反演结果的形状、位置和真实模型较为相近。整体反演结果可以准确反映出异常体的数量、形状、大小及其所处的中心位置。对反演结果精度的分析表明,图 3(g)中观测值Tzz和剩余密度值m计算出的Tzz′之间的差值较小,图 3(h)中残差值的标准差为±1.003 1 E,在有效范围内,表明算法可用于复杂模型的反演计算。
-
实测数据来源于澳大利亚西澳洲地区2009年建立的航空重力/重力梯度测量Kauring试验场,西澳洲矿产资源丰富,以富产铁矿、金矿、氧化铝等而全球闻名,且具有储藏量大、分布广、埋藏浅和易开发等特点。该试验场位于珀斯简达科特机场东北方向115 km处(图 4),该地区地势平坦,有利于试验的开展。
图 4 Kauring重力梯度试验场位置及该地区地形图
Figure 4. Location of Kauring Gravity Gradient Test Site and the Topographic Map of Area
Fugro Airborne Surveys于2011年7月至2012年2月间,利用FALCON航空重力梯度仪(airborne gravity gradient, AGG)在该地区进行了航空重力梯度测量,在Kauring试验场飞行的总测线长度为1 265.3 km,飞行垂直高度为80 m,使用的坐标系统为参心大地坐标系。实测重力梯度数据已经过地形改正,为突出浅部异常,经过参数试验并结合之前的研究成果[14],图 5(a)~5(f)使用相匹配的滤波法消除由深部异常产生的影响。北向范围为6 469~6 473 km,东向范围为504~508 km,所选测区范围大小为4 000 m×4 000 m。
图 5 巴特沃斯四阶低通滤波后的重力梯度全张量
Figure 5. Gravity Gradient Tensor by Using Butterworth Four Order Low Pass Filter
实测航空重力梯度数据的分辨率为5.1 m,精度为5 E;在采集数据中共选取80×80个点作为观测数据,观测点间隔为50 m,利用5个独立分量进行联合反演,故由观测值组成的矩阵大小为32 000×1;将模型划分为10层,每层厚度为80 m。参照数值模拟中模型设计的方法,地下模型单元大小从80 m递增到150 m,地下模型共由11 546个模型单元组成,核函数G为32 000×11 546的矩阵,初始模型矩阵m的值设为0。图 6(a)给出反演结果的三维模型图,图 6(b)至图 6(d)给出异常体中心位置的横纵截面。为突出剩余密度异常值较大的异常体,在三维模型图上仅给出密度值大于1.0×103 kg/m3的异常块体。相比之前的结果,在反演得到的三维模型图 6(a)中发现中心异常块体周围还存在着多个大小不等的密度异常体单元,与之前已发现的中心块体沿线呈一定角度分布,这些异常体的密度值在1.25×103 kg/m3左右。在y=1 700 m处给出纵截面图,图 6(b)表明反演结果的中心异常块体位于地下深度240~400 m范围,反演出异常体的中心位置与刘金钊等[15]之前对该地区密度异常深度探测的结果相吻合。反演中心块体位置位于坐标系正东方向505.11~505.45 km范围内,正北方向6 470.45~6 470.89 km范围内,反演出中心块体的密度值在2.3×103 kg/m3左右,与之前用航空重力反演该地区异常体的密度值相近,验证了反演方法应用于实测数据的可靠性。
由于从反演出的三维模型图上看出异常块体主要位于第五层320~400 m范围,所以给出360 m的横截面切片图。从图 6(c)发现,几块明显的异常块体沿坐标点(2 000,0)呈arctan0.5=26.5°分布,异常块体分布的沿线距离长度在2.3 km左右。如图 6(d)所示,可将横截面显示出的异常块体划分为4块,反演出异常块体的中心分别位于(505 280, 6 470 670)、(505 610, 6 470 170)、(505 950, 6 469 500)和(506 170, 6 469 110)点处,异常体大小从300 m×300 m至300 m×500 m不等。由于异常块体中心深度位于360 m左右,距离地表较近,且该地区的地势平坦,在之前的地质研究中并未发现断层的存在,因反演的初值m均为0,并未加任何约束,反演出异常体的密度值最大可达到2.41×103 kg/m3。而该地区矿藏具有埋藏浅的特点,故推断是由于矿体的存在导致该地区存在明显的异常体。对比之前重力反演的结果[14],发现中心块体周围还存在着明显的异常块体。
利用改进算法,可反演出异常体的位置、大小、密度异常值,从而判定密度异常体的分布特征。
-
本文利用改进后的预条件共轭梯度算法对地下密度异常体进行三维密度结构反演,通过构建模型重力梯度场验证了方法的有效性和可靠性。并利用该算法对澳大利亚西澳洲Kauring地区实测重力梯度场进行了反演,给出了测区内密度异常体的分布信息。国内航空重力梯度测量还未大范围开展, 希望通过本文对国外实测航空重力梯度数据的处理与解算, 为今后国内重力梯度试验场数据的应用提供经验和基础;该方法在重力梯度矿产资源调查中具有广泛的应用前景,可以得到更加完整的密度异常体分布信息,从而进行更加深入的地质结构分析和解释工作。
Inversion of Three-Dimensional Density Structure Using Airborne Gradiometry Data in Kauring Test Site
-
摘要: 相对于传统的重力测量手段,重力梯度测量能够以更高的灵敏度和分辨能力反映出地下密度异常体的结构特征。由于拉格朗日经验参数在实测数据反演中存在不确定性,对预条件共轭梯度反演算法加以改进,利用L曲线的拐点值代替原反演算法中的拉格朗日经验参数作为正则化参数;为改善反演中存在的病态性问题并减弱核函数的快速衰减,将地下模型改进为不等间隔模型;为改善反演中解的非唯一性,利用重力梯度的5个独立分量进行联合反演;通过对澳大利亚Kauring试验场航空重力梯度张量进行联合反演,得到该地区异常体的三维密度分布,将重力梯度联合反演结果与之前的重力反演结果对比分析,发现在中心异常体附近沿线还分布着多个异常块体。结果表明,改进后的算法能够有效地利用实测重力梯度数据反演出密度异常体的分布信息。Abstract: Due to the uncertainty of the Lagrange empirical parameter, selecting empirical of parameters for diverse observed data sets introduces uncertainty into the results, which weakens the applicability of the inversion method. By using the turning point of the L curve to replace the Lagrange empirical parameter as the regularization parameter, the algorithm focusing on the preconditioned conjugate gradient algorithm has been improved. The underground models have been converted to models with unequally spaced aiming to solve ill conditioned problem as to well as weaken kernel function attenuation. In order to take full advantage of the gravity gradient multiple components, the method of joint five independent measured components of tensor gradient gravity data has been taken with the purpose of meliorating the non uniqueness of inversion results. The effectiveness and reliability of the improved method are validated by the statistical analysis of multiple sets of synthetic models. For the application of the field data, analysis result shows that the improved calculation method is effectively applicable to the inversion of measured gravity gradient data, through inversion of airborne gradiometry data on Australian Kauring test site, we obtained 3D distribution of underground density anomalies. According to the previous results of gravity data inversion, this paper verifies the effectiveness of the algorithm, and discovers more anomaly blocks besides the central anomaly blocks. Our results show that the improved algorithm using field measurements can inverse the distribution of density anomalies, and the inversion results provide more detailed and reliable pattern information for the density anomaly.
-
表 1 模型统计信息
Table 1. Statistics of Model Information
密度异常体 尺寸(L×W×H)/m3 密度差/(kg·m-3) 质心平面/m 质心位置/m 红色异常体A 200×200×200 1 000 (500, 500) 200 红色异常体B 200×200×200 800 (-125, 0) 200 蓝色异常体C 200×200×200 -800 (375, 0) 200 -
[1] Everton B. The Use of the GOCE Mission Data for Characterizations and Implications on the Density Structure of the Sedimentary Basins of Amazon and Solimoes, Brazil[D].Trieste: Università Degli Studi Di Trieste, 2012 [2] 曾华霖.重力场与重力勘探[M].北京:地质出版社, 2005:68-71 Zeng Hualin. Gravitational Field and Gravitational Prospecting[M]. Beijing:Geological Publishing House, 2005:68-71 [3] Zhang Changyou, Mushayandbvu M F, Reid A, et al. Euler Deconvolution of Gravity Tensor Gradient Data[J].Geophysics, 2000, 65(2):512-520 doi: 10.1190/1.1444745 [4] Zhdanov M S, Ellis R, Mukherjee S. Three-Dimensional Regularized Focusing Inversion of Gravity Gradient Tensor Component Data[J].Geophysics, 2004, 69(4):925-937 doi: 10.1190/1.1778236 [5] Martinez C, Li Y, Krahenbuhl R, et al. 3D Inversion of Airborne Gravity Gradiometry Data in Mineral Exploration:A Case Study in the Quadrilatero Ferrifero, Brazil[J].Geophysics, 2012, 78(1):B1-B11 [6] 蒋甫玉, 高丽坤, 黄磷云.油气模型的重力梯度张量研究[J].吉林大学学报(地球科学版), 2011, 41(2):545-551 http://d.old.wanfangdata.com.cn/Periodical/cckjdxxb201102031 Jiang Fuyu, Gao Likun, Huang Linyun. Study on Gravity Gradient Tensor of Oil-Gas Model[J].Journal of Jilin University(Earth Science Edition), 2011, 41(2):545-551 http://d.old.wanfangdata.com.cn/Periodical/cckjdxxb201102031 [7] 孙中科, 朱自强, 鲁光银, 等.基于粒子群算法的重力张量反演研究[J].物探化探计算技术, 2013, 35(2):128-133 doi: 10.3969/j.issn.1001-1749.2013.02.02 Sun Zhongke, Zhu Ziqiang, Lu Guangyin, et al. Research on Gravity Tensor Inversion Based on Particle Swarm Optimization[J].Computing Techniques for Geophysical and Geochemical Exploration, 2013, 35(2):128-133 doi: 10.3969/j.issn.1001-1749.2013.02.02 [8] Hansen P C. Analysis of Discrete Ill-posed Problems by Means of the L-curve[J].SIAM Review, 1992, 34:561-580 doi: 10.1137/1034115 [9] Constable S C, Parker R L, Constable C G. Occam's Inversion:A Practical Algorithm for Generating Smooth Models from Electromagnetic Sounding Data[J]. Geophysics, 1987, 52(3):289-300 doi: 10.1190/1.1442303 [10] Li Yaoguo, Oldenburg D W. 3-D Inversion of Gravity Data[J]. Geophysics, 1998, 63:103-119 http://d.old.wanfangdata.com.cn/OAPaper/oai_doaj-articles_a36c0a2315d6fd28e329717dfe4b3131 [11] Pilkington M. 3D Magnetic Data-Space Inversion with Sparseness Constraints[J].Geophysics, 2009, 74(1):L7-L15 doi: 10.1190/1.3026538 [12] 姚姚.地球物理反演基本理论与应用方法[M].武汉:中国地质大学出版社, 2002:57-59 Yao Yao. Basic Theory and Application of Geophysical Inversion Methods[M]. Wuhan:China University of Geosciences Press, 2002:57-59 [13] 管志宁.地磁场与磁力勘查[M].北京:地质出版社, 2005:199-204 Guan Zhining. Magnetic Field and Magnetic Exploration[M]. Beijing:Geological Publishing House, 2005:199-204 [14] Martinez C, Li Y G. Understanding Gravity Gra-diometry Processing and Interpretation Through the Kauring Test Site Data[C]. 22nd ASEG International Geophysical Conference and Exhibition, Brisbane, Australia, 2012 [15] 刘金钊, 柳林涛, 梁星晖, 等.重力梯度特征向量和多尺度分析法在密度异常深度探测中的应用[J].武汉大学学报·信息科学版, 2016, 41(3):323-330 http://ch.whu.edu.cn/CN/abstract/abstract4570.shtml Liu Jinzhao, Liu Lintao, Liang Xinghui, et al. Application of Density Anomaly Depth Detection Using Gravity Gradient Eigenvectors and Multiscale Analysis Approach[J]. Geomatics and Information Science of Wuhan University, 2016, 41(3):323-330 http://ch.whu.edu.cn/CN/abstract/abstract4570.shtml -