当前位置：文档之家 › 矩形棱柱体重力梯度张量异常正演计算公式

矩形棱柱体重力梯度张量异常正演计算公式

  • 格式：pdf
  • 大小：1.26 MB
  • 文档页数：6

下载文档原格式

  / 6

合集下载

重磁(梯度)张量数据边界识别方法研究

重磁(梯度)张量数据边界识别方法研究

重磁(梯度)张量数据边界识别方法研究赵建宇;明彦伯;孙成城;张志东;刘晓甲【摘要】The boundary identification is one of the regular tasks of gravity and magnetic data interpretation,but gravity and magnetic raw anomaly data does not well correspond to the boundary of geological body.The maximum value of the horizontal derivative of gravity and magnetic anomalies are corresponding to the zero value of vertical derivative,which is used mostly to complete the assignment of gravity and magnetic boundary identification.We find that false boundary is caused by using current boundary identification method.This article summarizes the current boundary identification method and corrects the problem of excessive boundaries to present new boundary identification method.This new method can effectively remove excessive boundaries and reduces the disturbance of noise,which provides a new thinking way for gravity and magnetic data boundary identification method.%边界识别是重力数据解释的常规任务之一,但地质体的边界不能很好地与重力原始异常数据对应,而与重磁异常水平导数极大值、垂直导数零值相对应,因此大多利用该性质完成重磁边界识别任务.研究发现,现有边界识别方法会出现虚假边界,针对现有边界识别方法进行总结,并针对多余边界问题进行改进,提出新的边界识别,通过模型试验和实际数据证明,该方法可有效地去除多余边界,且降低了噪声的干扰,为重磁数据边界识别方法提供新的思路.【期刊名称】《物探化探计算技术》【年(卷),期】2017(039)006【总页数】7页(P748-754)【关键词】边界识别;重磁(梯度)张量数据;重磁勘探【作者】赵建宇;明彦伯;孙成城;张志东;刘晓甲【作者单位】吉林大学地球探测科学与技术学院,长春130021;吉林大学地球探测科学与技术学院,长春130021;吉林大学地球探测科学与技术学院,长春130021;吉林大学地球探测科学与技术学院,长春130021;吉林大学地球探测科学与技术学院,长春130021【正文语种】中文【中图分类】P631.2重磁勘探是一种可以有效地圈定异常，划分构造的物理勘探方法，由于它的经济、快速、范围广等特点而得到广泛应用。

航空重力梯度测量技术研究

航空重力梯度测量技术研究

航空重力梯度测量技术研究航空重力梯度测量作为二十世纪末发展起来的尖端技术,随着测量系统和处理解释方法的逐步完善,在固体矿产和油气资源勘查中发挥着日益重要的作用,并因其快速、高效和高空间分辨率等特点而备受青睐。

航空重力梯度测量技术是目前国际研究热点和难点,成熟的商业勘探技术为美国Lockheed Martin公司垄断,我国在该领域起步较晚,基础相对薄弱。

2006年开始,国家863计划开始支持航空重力梯度关键技术研究,经过十多年的努力,国内多个研究团队在多项关键技术上取得了重大突破,并在“十二五”期间实现了实验室静基座条件下重力梯度效应的测量,加快了该项技术的实用化进程。

本文围绕突破航空重力梯度测量系统研制关键技术及测量结果实际应用开展研究。

首先,系统调研了国外航空重力梯度测量技术研发历程、应用现状和研究动态,详细剖析了旋转加速度计航空重力梯度仪的测量原理和设计思想,梳理了关键技术难点及解决方案,跟踪了系统完善过程中的各项技术改进,为航空重力梯度测量系统自主研制和持续改进提供了参考和借鉴。

立足国内基础,制定了基于石英挠性加速度计部分重力梯度张量测量系统总体研制方案。

突破多项关键技术,研制完成的重力梯度仪用高分辨率加速度计样机分辨率优于1×10<sup>-8</sup>g,重力梯度敏感器实验室测量精度优于70E,重力梯度稳定平台满足载荷要求,性能指标通过飞行测试。

完成航空重力梯度测量系统集成、减震和温控方案设计,为“十三五”航空重力梯度测量系统飞行试验和实用化奠定了基础。

针对航空应用和在研航空重力梯度测量系统特点,优选Y-12飞机平台,开展了典型航空地球物理勘探条件下的飞机振动、姿态、气压、温度和湿度等环境状态参数测量及研究,详细分析了飞机底板振动的频率特征,揭示了振动信号的周期分布及振动周期与螺旋桨转速基频之间倍频关系的基本规律,总结了不同飞行状态下飞机侧滚、侧滑姿态角的变化特点及变化范围,分析了机舱内气压、温度和湿度随飞行过程的变化情况,为航空重力梯度测量系统量程、结构、减震、温控和气密设计及后续改进完善提供了参考和依据。

重力异常正反演问题

重力异常正反演问题

正演问题的定义： 根据巳知的、具有剩余质量的地质体的形状、产状和剩余密度 分布，通过理沦计算，研究它们所引起的异常及其各阶导数异 常的数值大小、空间分布和变化规律。 反演问题的定义： （1）由观测上重力异常的分布，在给定物体边界位置函数的条 件下，求解物体的密度分布函数；（物性反演） （2）由观测面上重力异常分布，在给定物体密度函数的条件下， 求解物体的边界位置的数值；（几何反演） （3）由观测面上重力异常分布。在给定特殊约束（如设物体密 度均匀、形态规则)条件下，求解物体密度参数和几何参数。 给定的函数和特殊约束称为反演问题的定解条件。
用解析公式计算出每个小长方
最后，将所有长方体的重力异
常值累加，以求得整个地质体在 计算点的异常值。
体在计算点所产生的重力异常值。
点元法 “点元”法所取的各个点元的体积可以相同，也可不同。各 个点元的物性可以相同，也可不同。通常是将勘探剖面之间的 地质体用适当的长方体或立方体来近似，确定出各个点元的角 点坐标，即可计算出该点元的三重积分值。 对于一个点元而言，其计算公式如下：
(j-1)
。由（4-6）式可知，由ρ
(j-1)
产生的重力场频谱为 F[△g(j-1)]为
n n F [ ( j 1) ( DH L ( r ) DH u ( r ))]
F g ( j 1) 2G
n 1



(－ k ) n 1 n!
e
k zc
（4-9）
而已知场△g（r0,z0）的频谱 F[△g]也可由（4-6）式来表示、将 F[△g]与 F[△g(j-1)]相减并经整 理后可得
2．面元法
用一组垂直于z轴的平面
或者垂直于X轴、y轴的平 面切割地质体，地质体与平 面相交形成一系列的裁面。

重力正演、反演

重力正演、反演

2)当σ>o时，极大值一侧对应着上升盘，极小 值一侧对应着下降盘，在极小值十分清晰且大 干极大值的绝对值时，属正断层类型，反之则 属逆断层类型。
二度铅垂柱体 对于沿水平方向延伸较长而横截面近于矩形的 矿脉，可以当成二度铅垂柱体来研究。在正演 它的异常时，坐标系及有关参数的选取见图,用 (x+α)与(x一α)分别代替铅垂台阶各公式中的 x，并将结果相减，即获得这一形体的重力异 常及各阶导数异常的公式:
当柱体的下底 H→+∞ 时，便可获得底部无限延 伸的铅垂脉的相应公式Δg→∞
( x − a) 2 + h 2 V xz = Gσ ln ( x + a) 2 + h 2 h h 2ah V zz = 2Gσ (tg −1 − tg −1 ) = 2Gσtg −1 2 x−a x+a x + h2 − a2 ⎡ ⎤ x+a x−a 2a ( a 2 + h 2 − x 2 ) V zzz = 2Gσ ⎢ = 2Gσ 2 − 2 2 2 2 ⎥ ( x + a) + h ⎦ ( x + a 2 + h 2 ) 2 − 4a 2 x 2 ⎣ ( x + a) + h
GM GMD = 2 2 nD ( x1 / n + D 2 ) 3 / 2
x 1/n = ± D n 2 / 3 − 1
取n=2，得x1/2=0.766D(X正半轴)和x’1/2=－0.766 D (X负半轴)，说明异常半极值点的横坐标为球心 深的0.766倍
4、当D不变，使M加大m倍时，异常也同样加大
[( x + a ) 2 + H 2 ][( x − a ) 2 + h 2 ] V xz = Gσ ln [( x + a ) 2 + h 2 ][( x − a ) 2 + H 2 ] H h H h ) − tg −1 − tg −1 + tg −1 V zz = 2Gσ (tg −1 x+a x+a x−a x−a ⎡ ⎤ x+a x+a x−a x−a − + − V zzz = 2Gσ ⎢ ⎥ 2 2 ( x + a) 2 + H 2 ( x − a) 2 + h 2 ( x − a) 2 + H 2 ⎦ ⎣ ( x + a) + h

两种新的长方体重力异常正演公式及其理论推导

两种新的长方体重力异常正演公式及其理论推导

+(
z) arctan ( z)R | 2 | 2 | 2 ( x)( y) 1 1 1
g(x, y, z) = G ||| ( x) ln{( y) + R} + ( y) ln{( x) + R}
( z) arctan ( x)( y) | 2 | 2 | 2 ( z)R 1 1 1
g(x, y, z) = G ||| ( x) ln{( y) + R} + ( y) ln{( x) + R}

两种新的长方体重力异常正演公式及其理论推导
骆遥 1, 2
1 中国科学院地质与地球物理研究所，北京（100029） 2 中国科学院研究生院，北京（100049） E-mail：geo@
摘 要： 在前人推导长方体重力场、磁场正演理论表达式工作的基础上，重新对长方体重 力场正演理论表达式进行理论推导，提出了两种全新的长方体重力异常正演公式形式，并给 出了全部的理论推导过程，对比模型正演计算结果表明，新导出长方体重力场正演理论表达 式的正确。 关键词：长方体，重力场，正演，积分 中图分类号：P631
线数据单位为 g.u.
Fig2. The cubic model gravity contour map
5. 结论
综合前人对长方体重力场正演理论表达式的推导过程，并借鉴长方体磁场及其梯度场理 论表达式的推导，推导出了两种新的长方体重力场正演理论表达式（11）式和（13）式，对 比模型正演计算结果表明，新导出长方体重力场正演理论表达式是完全正确的。
似积分的推导[15~17]，对 2 的推导有：
2= 2( 1
z)2 ( {(
y) {( x) +
R}2

重力数据处理解释方法

重力数据处理解释方法

2. 简单规则形体的异常特征及应用
2)水平圆柱体的重力异常及特征应用
2Gh l 1 当l2 l1 l时,g 2 1/ 2 2 2 2 2 x h x l h 3 x 2 2l 2 3h 2 2Ghxl Wxz 2 2 2 2 2 2 3/ 2 x h x l h Wzz
l2 1/ 2 2

l1
2. 简单规则形体的异常特征及应用
2)水平圆柱体的重力异常及特征应用
Wxz 3G Wzz 3G
l2 l2
l1
Ghx 3x 2 3h d 2 2 2 5/ 2 2 2 2 2 2 2 3/ 2 x h x h x h xh

d
Wxz 3G Wzz 3G
l2
l2
l1
2 2 x h
xh
2 5/ 2

d ,
l1
2 2 x h
2 x2
2 5/ 2

d
G 6.67 108 cm3 / g s 2 , R 2 , dm d
2. 简单规则形体的异常特征及应用
2
G
M
2 1/ 2
,Wxy Wyz 0
g 2W xh Wxz 3GM 2 x zx ( x h 2 )5 / 2 g 2W 2h 2 x 2 W GM 2 z zz ( x h 2 )5 / 2 Wzzz 2 g 2W 2h 2 3 x 2 3 3GM 2 2 z z ( x h 2 )7 / 2
0 M t M (1 )
2. 简单规则形体的异常特征及应用

重力及重力梯度张量地形效应方法

第2章重力及重力梯度张量地形效应方法2.1基本原理重力勘探是指根据地球内部物质密度分布不均匀引起的重力变化来研究地质构造并进行矿产勘探[70]。

重力勘探的流程如下图2-1所示：先对野外测量得到的观测数据进行数据处理，得到重力异常，而后通过正演得到密度差，进而进行地质解释。

而本文所做的重力地形效应就是数据处理中的一部分。

图2-1重力勘探流程图重力测量数据的处理流程如下图2-2所示：获得了野外重力测量值后，先进行仪器零点校正，再进行高度校正和正常场校正，得到自由空气重力异常；再进一步对自由空气重力异常进行地形校正和中间层校正，最终得到布格重力异常。

而重力测量数据处理流程的地形改正环节中最重要就是起伏地形条件下重力效应的计算。

图2-2重力测量数据的处理流程图总地形效应是重力势及其在某一点上由于大地水准面上所有地形质量而产生的导数(图2-3)。

剩余地形效应是由于地形质量相对于通过点的水平表面，或相对于地形的其他近似。

也就是说，如果我们用这个水平表面(或与地形相似的其他近似)来构造一个地形，已经包含和尚未包含的附加质量就构成了产生剩余地形效应的剩余质量。

对于地球表面上的观测点P，布格板是位于大地水准面和经过该点的水平面之间的无限平面，两者近似为水平面。

地形改正是考虑(即去除)残留地形的引力作用，即实际地形质量与点P 布格板的偏差(如图2-3)。

对于一般的重力测量应用(重力或重力梯度测量)，我们将地形改正定义为适当的重力效应，即去除计算点平面以上的质量，并添加平面以下的质量，以创建布格板。

质量的去除对应于密度为负号的效应；而在下面质量的增加则对应于为正号的效应。

因此，地形改正可以解释为由平面以上的负密度和平面以下的正密度而产生的引力效应。

图2-3布格板地形的地形效应(白色区)和残余地形(红色区)的地形改正2.1.1重力地形效应理论基础某起伏地形的剖面如图2-4所示，表示了观测点A周围起伏地形的影响。

当地形相对平整时，比A点高的地形质量对A点产生的垂直方向的分力会使A点的重力值减小；比A点低的地形，由于缺少物质，也会使A点的重力值减小。

重力异常正反演问题


设：
（D.1）
所以：
按式（D.1）的形式累加起来，最后只需要求一次反正切函 数，这样处理后，计算速度提高一倍以上。
1.2.2 直立“线元”法
某工区物探、勘探工作布置示意图
正演问题的定义： 根据巳知的、具有剩余质量的地质体的形状、产状和剩余密度 分布，通过理沦计算，研究它们所引起的异常及其各阶导数异 常的数值大小、空间分布和变化规律。 反演问题的定义： （1）由观测上重力异常的分布，在给定物体边界位置函数的条 件下，求解物体的密度分布函数；（物性反演） （2）由观测面上重力异常分布，在给定物体密度函数的条件下， 求解物体的边界位置的数值；（几何反演） （3）由观测面上重力异常分布。在给定特殊约束（如设物体密 度均匀、形态规则)条件下，求解物体密度参数和几何参数。 给定的函数和特殊约束称为反演问题的定解条件。
什么是正问题与反问题？
反问题：m=G-1d
观测数据d
地质模型 m
正问题：d=Gm
（一）规则形体的正、反演问题
为了简化，假设地质形体孤立存在，密度均匀，地 面水平，所取剖面为中心剖面
规则形体：球体、水平圆柱体、垂直台阶、脉状体……
1、球体
规 则 形 体 的 正 、 反 演 问 题
近似于等轴状地质体，如盐丘、矿巢、溶洞等
lim g 0; g max
g max h02 m G 2 ; x1/ 2 0.766h0 ; m h0 G
2、水平圆柱体 2、水平圆柱体（线质量）
规 则 形 体 的 正 、 反 演 问 题
小柱体元在P(x,0,0)点产生的重力异常为
g G
h0 dy
(x y h )
用一组垂直于y轴的平面
和一组垂直于X轴的平面分 别切割地质体，则任意两 个平面的交线包合在地质 体之内的部分形成一个线 元。 用解析式计算每一个线 元在计算点产生的重力异 常作用值。 对所有钱元的作用值依 次进行X方向和Y方向的数 值积分，便得到整个地质 体在计算点所产生的重力 异常值。

16重力勘探-重力异常正反演解析


△gz
△g
FHale Waihona Puke rh1R0
测量垂直梯度原理 gz h2 h1
g h2 g h1
g ( z z ) g ( z )) g h1 g h2 z h2 h1
△g △g(x+△x) △g(x-△x)
△gx
A（x,0,0）

△g
F r
h
1
R
0
2 x1
2 3
2 3
2 3
h2
n 2 3
（n 1 ）
x1 h
n
（n 1）
2 3
• （3）反演剩余质量
m g max G 2 h
• （4）反演半径
g max h 2 m G
g max h 2 m G
4 3 m R 3
3 m R( ) 4
• 重力异常的正反演（正反演问题的关系：异常场源（地
形状
质因素产生的剩余质量）和重力异常之间的对应关系（互相关系）包 括数量上关系。 ）
大小
异常场源 位置 产状 深度 物性
根据数学物理方法： 万有引力 重力异常的 推断：定性或者定量 △g
分布规律 形态特征 幅度大小
A F △g
1）正问题是反问题的基础； 2）反问题强烈依赖于正问题。
2 7 2 2
②均匀的水平圆柱体（二维水平柱状体）
• • • • • •
在实际的地质现象中，如长轴背斜、向斜等， 可以近似看成水平圆柱体来讨论。 水平柱状体：向两端无限延伸 半径：R 埋深：h 延伸方向：y 剩余密度
定义：线密度
S R 2 dd
hx g x 4G 2 ( x h2 )2 h2 x2 g z 2G 2 ( x h2 )2 h 2 3x 2 g zz 4Gh 2 ( x h 2 )3

基于多种位场数据及3D反演结果的综合解释方法


关键词
重力棱柱体模 型
Ｐ ５ ５ ４ ；
重力梯度数据
文献标志码 Ａ
重力数据 ３ Ｄ反演
位场数据综合解释
中图法分类号
重 力 位 沿 着 垂 直 方 向 的 一 阶导 数— — 重 力 数 据在 位 场 勘 探 观 测 数 据 解 释 中起 到 重 要 的作 用 。
测数 据 的解 释 中。基 于 各 类 重 力 数 据 处 理 及 解 释
ｄ
（ ２ ）
式（ ２ ） 中ｒ 是从 点 Ｐ到体 积 元 的距 离 ， Ｇ是 万有
ห้องสมุดไป่ตู้
３ １ 期
田 维， 等： 基 于多种位场数据及 ３ Ｄ反演结果 的综合解 释方法
引力常 量 。因为 重 力位 是 处 处 连 续 和有 限 的 ， 所 以 重 力也 是处 处连 续 和有 限 的。 设 坐标 系 中 轴 垂直 地面 指 向地 下 ， 和 Ｙ 互 相
的积分上 下 限代入 式 （ ４ ） ， 相 应 的重力 一 阶导数
如下 ：
ｒ
垂直 ， 且ｘ ｏ ｙ 平面垂直 轴 ； 地面观测点坐标 Ｐ的坐
标为 （ ， ｙ ， ｚ ） 。 早期 的重 力 仪 测 量 的是 重 力位 垂 直 分量 ， 由上述 可 知 ， 现 用 ｇ表 示 重 力 大小 ， 为 了简化
工作中， 并将虎林盆地划分为 ８ 个构造单元 。
工作中¨ Ｉ ２ Ｊ 。近几年随着 国外位场数据观测仪器 的 快速发展 ， 位场数据的二 阶导数已可 以通过仪器观
测 获得 ； 因此 ， 发 展 基 于 位 场 数 据 梯 度 张 量 的解 释
１ 方法原理
１ ． １ 棱柱 体 的重 力数据计 算方 法
  1. 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

（ ξ ｉ － ｘ） （ ζ ｋ － ｚ） （ η ｊ － ｙ） ｒ ｉｊｋ
，
㊀ ㊀ 收稿日期： ２０１５⁃０６⁃０８
㊀ ㊀ 基金项目： 国家高技术研究发展计划（２０１１ＡＡ０６０５０１） ㊁中国地质调查局地质矿产调查评价项目（１２１２０１１０８７０１３）
㊃ １２１８㊃
Ｕ ｚｚ ＝ － Ｇ ０ ρð ð ð ｓ㊃ａｒｃｔａｎ
∬ （ ｘᶄ ʏʏ
∬ （ ｘᶄ
２
，带入式（７） 得：
３／ ２
２
１ ｄｙᶄｄｚᶄ㊂ ＋ ｙᶄ２ ＋ ｚᶄ２ ） ３ ／ ２
＋ ｙᶄ２ ＋ ｚᶄ ）
ｄｙᶄｄｚᶄ ＝ （８）
㊀ ㊀ 由于对角分量的公式相似， 简单交换积分次序 ｙᶄｚᶄ ， ｘᶄｒ ｘᶄｚᶄ ， ｙᶄｒ （１３）
－ Ｇ ０ ρａｒｃｔａｎ
ｙᶄｚᶄ ㊂ ｘᶄｒ
Ｖ
Ｕ ｙｙ ＝ Ｇ ０
以下从基本公式出发推导矩形棱柱体解析公
Ｕ ｙｚ ＝ Ｇ ０ Ｕ ｚｚ ＝ Ｇ ０
∭ρ ３ｙᶄ ｒ
２ Ｖ Ｖ
万有引力定律，在笛卡尔坐标系中体密度为 ρ （ ξ，η，
重力梯度是重力位的二阶导数，写成张量形式为 ∂２ Ｕ é ê ê ∂ｘ ２ Ｕ ｘｚ ù ê ú ∂２ Ｕ Ｕ ｙｚ ú ＝ ê ú ê ∂ｙ ∂ｘ Ｕ ｚｚ ú û ê ２ ê∂ Ｕ ê ë ∂ｚ ∂ｘ ∂２ Ｕ ∂ｘ ∂ｙ
∭ ρ（ ξ，ｒη，ζ） ｄξｄηｄζ，
（３）
∂２ Ｕ ù ú ∂ｘ ∂ｚ ú Ｕ Ｕ ｘｙ é ú ê ｘｘ ∂２ Ｕ ∂２ Ｕ ú ê Ｕ ｙｘ Ｕ ｙｙ ㊂ （４） ∂ｙ ２ ∂ｙ ∂ｚ ú ê ê ú ë Ｕ ｚｘ Ｕ ｚｙ ∂２ Ｕ ∂２ Ｕ ú ∂ｚ ∂ｙ ∂ｚ ２ ú û ㊀ ㊀ 对于图 １ 所示密度均匀矩形棱柱体而言， 将式 （３） 代入式（４） ，并令 ρ （ ξ， η， ζ ） ＝ ρ㊁ ξ － ｘ ＝ ｘᶄ㊁ η － ｙ ＝ ６ 个重力梯度张量分量计算公式： ３ｘᶄ２ － ｒ ２ Ｕ ｘｘ ＝ Ｇ ０ ρ ｄｘᶄｄｙᶄｄｚᶄ， ｒ５ Ｖ Ｕ ｘｙ ＝ Ｇ ０ ｙᶄ㊁ζ －ｚ ＝ ｚᶄ，则 ｄξ ＝ ｄｘᶄ㊁ｄη ＝ ｄｙᶄ㊁ｄζ ＝ ｄｚᶄ， 可得常用的
㊀ ㊀ 上述积分方法， 在推导过程中前两步积分都是 严格成立的，但第三步对 ｚᶄ 积分过程并不严谨㊂ 首 ＝ ｃ２ ＝ ０， 此时 ， 不满足变量替换条件 ， 由此得出的 先，在上半无源空间，矩形棱柱体顶面角点在计算平 面上的投影点存在ｘᶄ㊁ｙᶄ同时为零的情况，即：ｘᶄ２ ＋ｙᶄ２ 积分结果无意义，为 解析奇点 ； 其次， 在积分倒数 第二步分子分母同除 ｘᶄ２ 亦不严谨，因为在矩形棱柱 体顶面角点在计算平面上的投影点存在 ｘᶄ ＝ ０ 的 情况㊂
３ｚᶄ２ － ｒ ２ ｄｘᶄｄｙᶄｄｚᶄ㊂ （５） ｒ５ Ｖ ㊀ ㊀ 根据重磁位场泊松公式， 式 （ ５） 与三分量磁异 常的核函数完全一致 ［９］ ㊂ 实际应用中， 仅需计算矩 形棱柱体上半无源空间的重力梯度张量㊂ 同时， 由 于 Ｕ ｘｘ ㊁Ｕ ｙｙ ㊁Ｕ ｚｚ 及 Ｕ ｘｙ ㊁Ｕ ｘｚ ㊁Ｕ ｙｚ 分别在形式上具有相似 特征，只有推导出其中任意两个，然后再进行简单变 ρ 量替换便可获得其他几个张量分量的积分结果㊂ 对 Ｕ ｘｘ 采用分步积分和变量替换方法先计算其 三重不定积分： ３ｘᶄ２ － ｒ ２ Ｕ ｘｘ ＝ Ｇ ０ ρ ｄｘᶄｄｙᶄｄｚᶄ ＝ ｒ５ Ｖ Ｇ０ ρ ２ｘᶄ ｙᶄ ｚᶄ ∭ （ ｘᶄ ＋ ｙᶄ ＋ ｚᶄ ）
应用于重力梯度正演数值模拟之中㊂ 笔者研究发 常计算公式 ［１－５］ 均源自 Ｄｏｎａｌｄ Ｐｌｏｕｆｆ ［６］ ，该公式虽然 观测平面进行特殊剖分㊂ 因此，通过详细推导，剖析
图 １㊀
１㊀ 问题的提出
矩形棱柱体与外部场点相互关系示意
研究的文献 ［１－５］ 中，如图 １ 所示， 矩形棱柱体不含积 分限的解析公式为 （ η － ｙ） （ ζ － ｚ） Ｕ ｘｘ ＝ － Ｇ ０ ρａｒｃｔａｎ ， （ ξ － ｘ） ｒ Ｕ ｘｙ ＝ Ｇ ０ ρｌｎ [ （ ζ － ｚ） ＋ ｒ ] ， Ｕ ｘｚ ＝ Ｇ ０ ρｌｎ [ （ η － ｙ） ＋ ｒ ] （ ξ － ｘ） （ ζ － ｚ） ， （ η － ｙ） ｒ Ｕ ｙｚ ＝ Ｇ ０ ρｌｎ [ （ ξ － ｘ） ＋ ｒ ] ， （ ξ － ｘ） （ η － ｙ） Ｕ ｚｚ ＝ － Ｇ ０ ρａｒｃｔａｎ ， （ ζ － ｚ） ｒ Ｕ ｙｙ ＝ － Ｇ ０ ρａｒｃｔａｎ
ｉ＝１ ｊ＝１ ｋ＝１ ２ ２
Ｕ ｙｚ ＝ Ｇ ０ ρð ð ð ｓ㊃ｌｎ [ （ ξ ｉ － ｘ） ＋ ｒ ｉｊｋ ] ，
２ ｉ＝１ ｊ＝１ ｋ＝１ ２ ２
２
２
２
（ ξ ｉ － ｘ） （ η ｊ － ｙ） （ ζ ｋ － ｚ） ｒ ｉｊｋ
２
，
数值模拟中，当观测平面网格剖分点与矩形棱柱体 四个角点在观测平面投影重合时， Ｕ ｘｘ 和 Ｕ ｙｙ 在该点 无计算值， 为解析奇点， 如图 ２ 所示㊂ 在实际应用 中，为了避免奇点出现，通常在观测平面采用特殊的 剖分方法，避开棱柱体角点在观测平面的投影点㊂
并利用类似变量替换便可得到其他两个分量公式， 因此按照上述积分方法得到对角分量解析表达式为 Ｕ ｘｘ ＝ － Ｇ ０ ρａｒｃｔａｎ Ｕ ｙｙ ＝ － Ｇ ０ ρａｒｃｔａｎ Ｕ ｚｚ ＝ － Ｇ ０ ρａｒｃｔａｎ
π ＜ｔ ２
１ é ù Ｕ ｘｘ ＝ － Ｇ ０ ρｘᶄ ê ê （ ｘᶄ２ ＋ ｙᶄ２ ＋ ｚᶄ２ ） ３ ／ ２ ｄｙᶄ ú ú ｄｚᶄ ＝ ë û １ é ù － Ｇ ０ ρｘᶄ ê ㊃ｂｓｅｃ ２ｔ ｄｔ ú ｄｚᶄ ＝ ３ ３ ê ú ë ｂ ｓｅｃ ｔ û
∬ʏ
可得 （７） Ｕ ｘｘ ＝ － Ｇ ０ ρａｒｃｔａｎ ｙᶄｚᶄ ｘᶄ ｘᶄ ＋ ｙᶄ２ ＋ ｚᶄ２
２
∬
＝ （１２）
ｘᶄ ｘᶄ ｘᶄ２ ＋ａ ２
２
㊀ ㊀ 同理，令 ｘᶄ２ ＋ｚᶄ２ ＝ ｂ ２（ ｂ ＞０） ，ｙᶄ ＝ ｂ ㊃ｔａｎ ｔ，（ － ＜ π ） ，则 ｄｙᶄ ＝ ｂ㊃ｓｅｃ ２ｔ ｄｔ，带入式（８） 得 ２
式（１） 的离散计算公式为
２ ２ ２
，
（１）
Ｕ ｙｙ ＝ － Ｇ ０ ρð ð ð ｓ㊃ａｒｃｔａｎ
ｉ＝１ ｊ＝１ ｋ＝１
Ｕ ｘｚ ＝ Ｇ ０ ρð ð ð ｓ㊃ｌｎ [ （ η ｊ － ｙ） ＋ ｒ ｉｊｋ ] ，
２ ｉ＝１ ｊ＝１ ｋ＝１ ２ ２
Ｕ ｘｙ ＝ Ｇ ０ ρð ð ð ｓ㊃ｌｎ [ （ ζ ｋ － ｚ） ＋ ｒ ｉｊｋ ] ，
舒晴１，２ ，朱晓颖２ ，周坚鑫２ ，高维２ ，尹航２
存在的原因，并基于前人提出的长方体 ΔＴ 场及其梯度场无解析奇点理论表达式，建立了矩形棱柱体重力梯度张量 中图分类号： Ｐ６３２㊀ ㊀ ㊀ 文献标识码： Ａ㊀ ㊀ ㊀ 文章编号： １０００－８９１８（ ２０１５） ０６－１２１７－０６
㊀ ㊀ 矩形棱柱体作为最常用的模型构建单元被广泛 现，目前被广泛引用的矩形棱柱体重力梯度张量异 简洁㊁对称，但文献中仅给出了结果， 没有详细地推 导过程，且存在解析奇点 ［７］ ， 在实际应用中需要对 解析奇点存在原因， 并寻找新途径构建新的无解析 奇点公式显得十分必要㊂
矩形棱柱体重力梯度张量异常正演计算公式
（ １． 吉林大学 地球探测科学与技术学院， 吉林 长春 ㊀ １３００２１；２． 中国国土资源航空物探遥感中心， 北京㊀ １０００８３）
摘 要：矩形棱柱体重力梯度张量异常常用的正演计算公式含有解析奇点㊂ 结合理论推导过程详细分析了 奇点 无解析 奇点 公式，通过模型计算结果对比和理论分析证明了其正确性㊂ 关键词：重力梯度张量异常；矩形棱柱体；解析奇点
在众多重力梯度张量分量正演数值模拟和反演
其 中， ｒ ＝
１０ －１１ ｍ ３ ／ （ ｋｇ㊃ｓ２ ） ，为万有引力常数㊂ Ｕ ｘｘ ＝ － Ｇ ０ ρð ð ð ｓ㊃ａｒｃｔａｎ
ｉ＝１ ｊ＝１ ｋ＝１ ２ ２ ２ ｉ＝１ ｊ＝１ ｋ＝１ ２ ２ ２
（ ξ －ｘ） ２ ＋（ η －ｙ） ２ ＋（ ζ －ｚ） ２ ， Ｇ ０ ＝ ６． ６７２ ˑ （ η ｊ － ｙ） （ ζ ｋ － ｚ） （ ξ ｉ － ｘ） ｒ ｉｊｋ
２
∭
ｄｘᶄｄｙᶄｄｚᶄ， ∭ρ ３ｙᶄｚᶄ ｒ
５
５
－ ｒ２
ｄｘᶄｄｙᶄｄｚᶄ，
∭
－
２
－
２
∭
Ｕ ｘｚ ＝ Ｇ ０
ｄｘᶄｄｙᶄｄｚᶄ， ∭ρ ３ｘᶄｙᶄ ｒ ｄｘᶄｄｙᶄｄｚᶄ， ∭ρ ３ｘᶄｚᶄ ｒ
Ｖ ５ Ｖ ５
㊀ ㊀ 详细推导过程如下：首先对 ｘᶄ进行积分，在上半 无源空间 ｚᶄ ＝ ζ －ｚ ＞０，利用变量替换，令 ｙᶄ２ ＋ｚᶄ２ ＝ ａ ２（ ａ π π ＞０） ，ｘᶄ ＝ ａ㊃ｔａｎ ｔ，（ － ＜ ｔ ＜ ） ， 则 ｄｘᶄ ＝ ａ ㊃ｓｅｃ ２ｔ ｄｔ， ２ ２ 代入得 ２ｘᶄ２ － ｙᶄ２ － ｚᶄ２ é ù ｄｘᶄ ú Ｕ ｘｘ ＝ Ｇ ０ ρ ê ｄｙᶄｄｚᶄ ＝ ２ ２ ２ ５／ ２ ê ú ＋ ＋ （ ｘᶄ ｙᶄ ｚᶄ ） ë û ２ａ ２ ｔａｎ ２ ｔ － ａ ２ é ù ２ Ｇ０ ρ ê ê （ ａ ２ ｔａｎ ２ｔ ＋ ａ ２ ） ５ ／ ２ ㊃ａ ｓｅｃ ｔｄｔ ú ú ｄｙᶄｄｚᶄ ＝ ë û
ｘᶄｙᶄ ㊂ ｚ得到 ｓｉｎ ｔ ＝ Ｕ ｘｘ ＝ － Ｇ ０ ρｘᶄｙᶄ
ʏ
ｓｉｎ ｔ ｂ２ ｙᶄ
ｄｚᶄ，
（９）
（１０） ㊀ ㊀ 参照前两步积分方法，利用变量替换，令ｘᶄ２ ＋ ｙᶄ２ ＝ ｃ２ ，当 ｘᶄ，ｙᶄ不同时为零时，ｃ ＞０，ｚᶄ ＝ ｃｔａｎ ｔ，（ － π ） ，则 ｄｚᶄ ＝ ｃ㊃ｓｅｃ ２ｔｄｔ，代入式（１０） 得 ２
Ｖ
２
２
２
５／ ２
ｄｘᶄｄｙᶄｄｚᶄ㊂
（６）