长方体场及其梯度场无解析关键点的理论公式

长方体场及其梯度场无解析关键点的理论公式

1膜分离器磁场梯度为了提高航空磁体研究的应用范围，以航磁梯度为主的航磁大质量分析和反演解释方法已成为国内外航磁机研究的发展趋势和方向。例如，在过去的十多年里，传统的空间梯度分析信号法（即总梯度模型法）和欧洲-美国的反演模拟方法越来越受到重视。在中国，1993年，中国人民大学（北京）对磁体进行了系统的理论研究。为了协助国内航磁干燥梯度数据的处理以及总梯度模型和优化方法的应用，中国航磁干燥动态模拟中心在1992年、1994年、1998年和2003年底在东海进行了多次试验测量。为了协助国内航磁干燥梯度数据的计算，以及采用国内航磁干燥梯度模型和优化反演方法的改进项目的研究，林场及其梯度模型的正演计算非常重要。体和长方体是正演计算中最常用的两种三级体模型。当使用文献中现有的长方体模型t及其梯度场的理论方程进行正演计算时，这些算法问题被忽视。为了不全取消去特定长方体模型的所有数据，某些特定长方体模型的高级建模数据被广泛收集。因此，需要进行深入的分析分析和研究“点”产生的原因，并提供相应的解决方案。这具有重要的理论和实际意义。2原理论的概念2.1e[2e]-z2+-y1-z2+-2/-z2+-2-z2+-2-y2+-2/-z2-2+-2-22+-2-22+-2-y2-22+-2-22+-2-22+-2e2-2+-2-22+-2e2-2+-2e2-2+-2e2-2+-2e2+-2e2-22+-2e2-22+-2e2-22+-2e2-22+-2e2+-2e2-22+-2e2-22+-2e2-22+-2e2+-2e2-22+-2e2-22+-2e2+-2e2+-2e2+-2e2+-2e2+-2e2+-2e2+-2e2+-2e2+-2e2+-2e2+-2e2+-2e2+-2e2+-2e2+-2e2+-2e2+-2e2+-2e2+-2e2+-2e2+-2e2+-2e2+-2e2+-2e2+-2e2+-2e2长方体ΔΤ场及其梯度场1)的正演计算理论表达式为ΔΤ(x,y,z)=-μ04π{Κ1arctan(ξ-x)(η-y)r(ζ-z)+Κ2arctan(η-y)(ζ-z)r(ξ-x)+Κ3arctan(ξ-x)(ζ-z)r(η-y)+Κ4ln[r-(ξ-x)]+Κ5ln[r-(η-y)]+Κ6ln[r-(ζ-z)]}|x0+a/2x0-a/2|y0+b/2y0-b/2|2z0-hh‚(1)∂ΔΤ(x,y,z)∂x=-μ04π{Κ1(η-y)(ζ-z)[(ξ-x)2-r2]r[r2(ζ-z)2+(ξ-x)2(η-y)2]+Κ2(η-y)(ζ-z)[(ξ-x)2+r2]r[r2(ξ-x)2+(η-y)2(ζ-z)2]+Κ3(η-y)(ζ-z)[(ξ-x)2-r2]r[r2(η-y)2+(ξ-x)2(ζ-z)2]+Κ41r+Κ5(ξ-x)r[(η-y)-r]+Κ6(ξ-x)r[(ζ-z)-r]}|x0+a/2x0-a/2|y0+b/2y0-b/2|2z0-hh‚(2)∂ΔΤ(x,y,z)∂y=-μ04π{Κ1(ξ-x)(ζ-z)[(η-y)2-r2]r[r2(ζ-z)2+(ξ-x)2(η-y)2]+Κ2(ξ-x)(ζ-z)[(η-y)2-r2]r[r2(ξ-x)2+(η-y)2(ζ-z)2]+Κ3(ξ-x)(ζ-z)[(η-y)2+r2]r[r2(η-y)2+(ξ-x)2(ζ-z)2]+Κ4(η-y)r[(ξ-x)-r]+Κ51r+Κ6(η-y)r[(ζ-z)-r]}|x0+a/2x0-a/2|y0+b/2y0-b/2|2z0-hh‚(3)∂ΔΤ(x,y,z)∂z=-μ04π{Κ1(ξ-x)(η-y)[(ζ-z)2+r2]r[r2(ζ-z)2+(ξ-x)2(η-y)2]+Κ2(ξ-x)(η-y)[(ζ-z)2-r2]r[r2(ξ-x)2+(η-y)2(ζ-z)2]+Κ3(ξ-x)(η-y)[(ζ-z)2-r2]r[r2(η-y)2+(ξ-x)2(ζ-z)2]+Κ4(ζ-z)r[(ξ-x)-r]+Κ5(ζ-z)r[(η-y)-r]+Κ61r}|x0+a/2x0-a/2|y0+b/2y0-b/2|2z0-hh‚(4)2.2关于三重体积分积分理论上,上半无源空间(即ζ-z>0、r>0)位场是不存在“奇点”的,式(2)～(4)由式(1)推导而出,文中将详细给出式(1)导出的积分过程.设(ξ,η,ζ)为长方体内体积微分元的坐标,而(x,y,z)为无源区计算点的坐标,为简便起见,省略公式前面因采用国际单位制时的系数μ0/(4π),由泊松公式可得该体积微分元产生的磁位和ΔΤ场:dU=-J⋅Δ(1r)⋅dξ⋅dη⋅dζ=-J∂∂t(1r)⋅dξ⋅dη⋅dζ‚dΔΤ=-∂∂t0(dU)=J∂2∂t∂t0(1r)dξdηdζ,式中:∂∂t0=l∂∂x+m∂∂y+n∂∂z,∂∂t=L∂∂x+Μ∂∂y+Ν∂∂z,l、m、n及L、M、N分别为地磁场和总磁化强度J的单位矢量t0、t的方向余弦。dΔΤ=J[(l∂∂x+m∂∂y+n∂∂z)(L∂∂x+Μ∂∂y+Ν∂∂z)(1r)]dξdηdζ=J{lL∂∂x[∂∂x(1r)]+mΜ∂∂y[∂∂y(1r)]+nΝ∂∂z[∂∂z(1r)]+(lΜ+mL)∂∂x[∂∂y(1r)]+(lΝ+nL)∂∂x[∂∂z(1r)]+(mΝ+nΜ)∂∂y[∂∂z(1r)]}dξdηdζ‚由于对任意f(ξ-x,η-y,ζ-z)类型的函数有∂f∂ξ=-∂f∂x,∂f∂η=-∂f∂y,∂f∂ζ=-∂f∂z,并且∂∂x(1r)=ξ-xr3‚∂∂y(1r)=η-yr3‚∂∂z(1r)=ζ-zr3,则上式积分后可以写成:ΔΤ(x,y,z)=J∫∫∫V{-lL∂∂ξ(ξ-xr3)-mΜ∂∂η(η-yr3)-nΝ∂∂ζ(ζ-zr3)+(lΜ+mL)∂∂ξ[∂∂η(1r)]+(lΝ+nL)∂∂ξ[∂∂ζ(1r)]+(mΝ+nΜ)∂∂η[∂∂ζ(1r)]}dξdηdζ‚(5)令Q1=∫∫∫V∂∂ζ(ζ-zr3)dξdηdζ,在求解该三重体积分的过程中将用到积分公式∫dv(v2+A2)3/2=vA2√v2+A2+C‚(A、C为常数,并且A≠0),并用到关系式r2=[(ξ-x)2+(η-y)2+(ζ-z)2],体积分的积分限范围为ξ1=x0-a/2,ξ2=x0+a/2,η1=y0-b/2,η2=y0+b/2,ζ1=h,ζ2=2z0-h,则有Q1=∫η2η1∫ξ2ξ1(ζ-zr3)dξdη|ζ2ζ1=∫η2η1∫ξ2ξ1(ζ-z)d(ξ-x)[(ξ-x)2+(η-y)2+(ζ-z)2]3/2dη|ζ2ζ1=∫η2η1(ζ-z)(ξ-x)[(η-y)2+(ζ-z)2]rdη|ξ2ξ1|ζ2ζ1=∫η2η1(ζ-z)(ξ-x)r[(ξ-x)2+(ζ-z)2][(η-y)2+(ζ-z)2][(ξ-x)2+(ζ-z)2]dη|ξ2ξ1|ζ2ζ1=∫η2η1(ζ-z)(ξ-x)r[r2-(η-y)2]r2(ζ-z)2+(ξ-x)2(η-y)2dη|ξ2ξ1|ζ2ζ1=∫η2η1r(ζ-z)(ξ-x)-(ζ-z)(ξ-x)(η-y)2/rr2(ζ-z)2+(ξ-x)2(η-y)2dη|ξ2ξ1|ζ2ζ1=∫η2η1(ξ-x)r(ζ-z)-(ξ-x)(η-y)2(ζ-z)/rr2(ζ-z)21+(ξ-x)2(η-y)2r2(ζ-z)2dη|ξ2ξ1|ζ2ζ1=∫η2η1d[(ξ-x)(η-y)r(ζ-z)]1+[(ξ-x)(η-y)r(ζ-z)]2|ξ2ξ1|ζ2ζ1=arctan[(ξ-x)(η-y)r(ζ-z)]|ξ2ξ1|η2η1|ζ2ζ1‚(6)如文献中所述,可推出下面积分∫∫∫V∂∂ξ(ξ-xr3)dξdηdζ=arctan[(η-y)(ζ-z)r(ξ-x)]|ξ2ξ1|η2η1|ζ2ζ1,(7)∫∫∫V∂∂η(η-yr3)dξdηdζ=arctan[(ξ-x)(ζ-z)r(η-y)]|ξ2ξ1|η2η1|ζ2ζ1,(8)(5)式中另外几项积分的推导过程中将用到积分公式:∫dv√v2+A2=-ln(√v2+A2-v)+C,(A、C为常数),则有:∫∫∫V∂∂ξ[∂∂η(1r)]dξdηdζ=∫ζ2ζ1∫η2η1∂∂η(1r)dηdζ|ξ2ξ1=∫ζ2ζ11rdζ|ξ2ξ1|η2η1=∫ζ2ζ1d(ζ-z)√(ζ-z)2+(ξ-x)2+(η-y)2|ξ2ξ1|η2η1=-ln[r-(ζ-z)]|ξ2ξ1|η2η1|ζ2ζ1‚(9)同理,可导出:∫∫∫V∂∂ξ[∂∂ζ(1r)]dξdηdζ=-ln[r-(η-y)]|ξ2ξ1|η2η1|ζ2ζ1,(10)∫∫∫V∂∂η[∂∂ζ(1r)]dξdηdζ=-ln[r-(ξ-x)]|ξ1ξ2|η1η2|ζ1ζ2,(11)将(6)～(11)式结果代入(5)式则可得到表达式(1).上面推导过程看似无懈可击,实际上问题出在(7)和(8)式的导出过程中及(9)式的导出结果上.以(7)式为例,在类似(6)式的推导过程中,在其第四步需要对分子和分母同乘以[(ξ-x)2+(η-y)2],在其第七步则对分子和分母同除以r2(ξ-x)2,但是当长方体顶面的四个角点在计算平面上的正投影点与网格点重合时(即x=ξ、y=η),即使在ζ-z>0、r>0的上半无源空间,由于这些网格点处这两个式子的值为零,因此这两步运算及导出的表达式在这四个特殊的网格点(即所谓的“奇点”)上是无意义和不成立的,虽然在其他网格点上是有意义和成立的.(8)式的导出过程也存在类似的问题,虽然(6)式的导出过程及结果完全正确,但类推出(7)～(8)式的过程和结果均存在所谓“奇点”问题,另外在“奇点”处,r-(ζ-z)值为零也导致(9)式的对数表达式无意义,这些就是(1)式的ΔΤ场及由它导出的(2)～(4)式的梯度场存在解析“奇点”问题的原因所在.3-z2+-z2+-z2+-z2+-z2+-z2+-z2+2+-2+-z2+-2+-z2+-z2+-2+-z2+-2+-z2+-2+-z2+-2+-z2+-2+-z2+-2+-z2+-2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-z2+-2+-z2+-z2+-z2+-z2+-z2+-2+-2+-2+-z2+-2+-z2+-z2+-为了解决长方体ΔΤ场及其梯度场正演计算原理论表达式的“奇点”问题,必须重新推导(7)～(8)式左端的积分,以得到详细的正确推导过程和无“奇点”的解析表达式.Q2=∫∫∫V∂∂ξ(ξ-xr3)dξdηdζ=∫ζ1ζ2∫η1η2(ξ-xr3)dηdζ|ξ1ξ2=∫ζ1ζ2∫η1η2(ξ-x)d(η-y)[(η-y)2+(ξ-x)2+(ζ-z)2]3/2dζ|ξ1ξ2=∫ζ1ζ2(ζ-x)(η-y)[(ξ-x)2+(ζ-z)2]rdζ|ξ1ξ2|η1η2=∫ζ1ζ2(ξ-x)(η-y)/r[r+(ζ-z)]2[(ξ-x)2+(ζ-z)2][r+(ζ-z)]2dζ|ξ1ξ2|η1η2‚为简化,令u2=[(ξ-x)2+(ζ-z)2],则Q2=∫ζ1ζ2(ξ-x)(η-y)/r[r2+2r(ζ-z)+(ζ-z)2]u2[r2+2r(ζ-z)+(ζ-z)2]dζ|ξ1ξ2|η1η2=∫ζ1ζ2(ξ-x)(η-y)/r[r2+2r(ζ-z)+(ζ-z)2]u2[(ξ-x)2+(ζ-z)2+(η-y)2+2r(ζ-z)+(ζ-z)2]dζ|ξ1ξ2|η1η2=∫ζ1ζ2(ξ-x)(η-y)/r[r2+2r(ζ-z)+(ζ-z)2]u2[u2+2r(ζ-z)]+u2[(ζ-z)2+(η-y)2]dζ|ξ1ξ2|η1η2‚上式中u2[(ζ-z)2+(η-y)2]=[(ξ-x)2+(ζ-z)2][(ζ-z)2+(η-y)2]=[(ξ-x)2+(ζ-z)2+(η-y)2](ζ-z)2+(ξ-x)2(η-y)2=r2(ζ-z)2+(ξ-x)2(η-y)2,将结果代入可得Q2=∫ζ1ζ2(ξ-x)(η-y)[r+2(ζ-z)+(ζ-z)2/r]u2[u2+2r(ζ-z)]+r2(ζ-z)2+(ξ-x)2(η-y)2dζ|ξ1ξ2|η1η2=∫ζ1ζ2(ξ-x)(η-y)[2(ζ-z)+r+(ζ-z)2/r][u2+r(ζ-z)]2+(ξ-x)2(η-y)2dζ|ξ1ξ2|η1η2=∫ζ1ζ2(ξ-x)(η-y)[2(ζ-z)+r+(ζ-z)2/r][(ξ-x)2+(ζ-z)2+r(ζ-z)]2+(ξ-x)2(η-y)2dζ|ξ1ξ2|η1η2=∫ζ1ζ2(ξ-x)(η-y)[2(ζ-z)+r+(ζ-z)2/r][(ξ-x)2+(ζ-z)2+r(ζ-z)]21+(ξ-x)2(η-y)2[(ξ-x)2+(ζ-z)2+r(ζ-z)]2dζ|ξ1ξ2|η1η2=∫ζ1ζ2(ξ-x)(η-y)[2(ζ-z)+r+(ζ-z)2/r][(ξ-x)2+(ζ-z)2+r(ζ-z)]21+[(ξ-x)(η-y)(ξ-x)2+(ζ-z)2+r(ζ-z)]2dζ|ξ1ξ2|η1η2=-∫ζ1ζ2d[(ξ-x)(η-y)(ξ-x)2+(ζ-z)2+r(ζ-z)]1+[(ξ-x)(η-y)(ξ-x)2+(ζ-z)2+r(ζ-z)]2|ξ1ξ2|η1η2=-arctan[(ξ-x)(η-y)(ξ-x)2+(ζ-z)2+r(ζ-z)]|ξ1ξ2|η1η2|ζ1ζ2‚(12)与上面类似的推导过程可推出∫∫∫V∂∂η(η-yr3)dξdηdζ=-arctan[(ξ-x)(η-y)(η-y)2+(ζ-z)2+r(ζ-z)]|ξ1ξ2|η1η2|ζ1ζ2,(13)为了避免(9)式可能出现的“奇点”问题,利用积分公式:∫dvv2+A2=ln(v2+A2+v)+C,可导出(9)～(11)式的另一种无解析“奇点”的表达式∫∫∫V∂∂ξ[∂∂η(1r)]dξdηdζ=ln[r+(ζ-z)]│ξ1ξ2│η1η2│ζ1ζ2,(14)∫∫∫V∂∂ξ[∂∂ζ(1r)]dξdηdζ=ln[r+(η-y)]│ξ1ξ2│η1η2│ζ1ζ2,(15)∫∫∫V∂∂η[∂∂ζ(1r)]dξdηdζ=ln[r+(ξ-x)]│ξ1ξ2│η1η2│ζ1ζ2,(16)将(6)、(12)～(16)式的结果代入(5)式中可获得长方体ΔT场无解所“奇点”的正演计算理论表达式(17).该表达式与Bhattacharyya给出的向下无限延深直立四方棱柱体公式有相似之处,但形式更为简单和整齐;另外Bhattacharyya给出的推导过程不详细.在式(17)基础上本文导出长方体ΔT梯度场无解析“奇点”的正演计算理论表达式(18)～(20)为ΔΤ(x,y,z)=J{k1ln[r+(ξ-x)]+k2ln[r+(η-y)]+k3ln[r+(ζ-z)]+k4arctan(ξ-x)(η-y)(ξ-x)2+r(ζ-z)+(ζ-z)2+k5arctan(ξ-x)(η-y)(η-y)2+r(ζ-z)+(ζ-z)2+k6arctan(ξ-x)(η-y)r(ζ-z)}|x0-a/2x0+a/2|y0-b/2y0+b/2|h2z0-h,(17)∂ΔΤ(x,y,z)∂x=J{-k11r-k2(ξ-x)r[(r+(η-y)]-k3(ξ-x)r[r+(ζ-z)]+k4(η-y)[(ξ-x)2-r(ζ-z)]r[(ξ-x)2+(ζ-z)2][r+(ζ-z)]-k5(η-y)r[r+(ζ-z)]-k6(η-y)(ζ-z)r[(ξ-x)2+(ζ-z)2]}|x0-a/2x0+a/2|y0-b/2y0+b/2|h2z0-h,(18)∂ΔΤ(x,y,z)∂y=J{-k1(η-y)r[r+(ξ-x)]-k21r-k3(η-y)r[r+(ζ-z)]-k4(ξ-x)r[r+(ζ-z)]+k5(ξ-x)[(η-y)2-r(ζ-z)]r[(η-y)2+(ζ-z)2][r+(ζ-z)]-k6(ξ-x)(ζ-z)r[(η-y)2+(ζ-z)2]}|x0-a/2x0+a/2|y0-b/2y0+b/2|h2z0-h,(19)∂ΔΤ(x,y,z)∂z=J{-k1(ζ-z)r[r+(ξ-x)]-k2(ζ-z)r[r+(η-y)]-k31r+k4(ξ-x)(η-y)r[(ξ-x)2+(ζ-z)2]+k5(ξ-x)(η-y)r[(η-y)2+(ζ-z)2]+k6(ξ-x)(η-y)[r2+(ζ-z)2]r[(ξ-x)2+(ζ-z)2][(η-y)2+(ζ-z)2]}|x0-a/2x0+a/2|y0-b/2y0+b/2|h2z0-h,(20)其中k1=mN+nM,k2=lN+nL,k3=lM+mL,k4=lL,k5=mM,k6=-nN,l=cosI0cosA0,m=cosI0sinA0,n=sinI0,L=cosIcosA,M=cosIsinA,N=sinI,I0、A0与I、A分别为地磁场与总磁化强度方向的倾角和相对X轴方向偏角。采用与上面类似的积分过程和方法也可导出长方体重力场的正演计算理论表达式.表达式(1)～(4)、(17)～(20)均是长方体的长、宽边分别与X、Y轴平行情况下的ΔΤ场及其梯度场正演计算公式,通过坐标的旋转、平移变换,令x=(x′-x0)cosφ+(y′-y0)sinφ,y=(y′-y0)cosφ-(x′-x0)sinφ,借助式(17)～(20)还可计算中心位置为(x0,y0,z0),走向与X′轴夹角为φ的顶面水平的长方体ΔΤ场及其梯度场数据.但在计算旋转后长方体的水平梯度时,要考虑两个水平求导方向的变化.需要指出的是ΔΤ场是磁异常场在地磁场方向上的分量场,如果用l、m、n表示任意磁分量方向,则式(17)～(20)实际上也是任意磁分量场(例如Hax、Hay、Za场等)及其梯度场正演计算的理论通式.4计算结果比较为了检验首次推出的无解析“奇点”的长方体ΔΤ场及其梯度场正演计算理论表达式的正确性,对表达式(1)～(4)、(17)～(20)的模型计算结果进行了比较.图2,3是斜磁化长方体有“奇点”和无“奇点”ΔΤ场及∂ΔΤ∂y场等值线图.比较结果表明,除4个解析“奇点”外,由式(1)～(4)和(17)～(20)计算的其他网格点磁场及梯度场数据(精确到两位小数情况下)完全一致,这充分说明新导出的长方体无解析“奇点”ΔΤ场及其梯度场正演计算理论表达式(17)～(20)是完全正确的.5投影点位置上的正演传统5.1文献中原有的长方体ΔT场及其梯度场正演计算理论表达式是上半无源空间存在场值无法计算的解析“奇点”的局限公式,对某些特定长方体模型无法获得完整的正演场数据.5.2新导出的长方体ΔT场及其梯度场无解析“奇点”正演计算理论表达式是上半无源空间不存在场值无法计算的解析“奇点”的普遍适用公式,对所