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

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