﻿ 多姿态具磁性未爆弹磁异常快速正演
 文章快速检索
 同济大学学报(自然科学版)  2018, Vol. 46 Issue (5): 687-694.  DOI: 10.11908/j.issn.0253-374x.2018.05.017 0

### 引用本文

WU Jiansheng, JIANG Wancong. Magnetic Anomaly Fast Forward-modeling of Magnetized Unexploded Ordnance with Different Shapes and Positions[J]. Journal of Tongji University (Natural Science), 2018, 46(5): 687-694. DOI: 10.11908/j.issn.0253-374x.2018.05.017

### 文章历史

1. 同济大学 海洋与地球科学学院, 上海 200092;
2. 同济大学 海洋地质国家重点实验室, 上海 200092

Magnetic Anomaly Fast Forward-modeling of Magnetized Unexploded Ordnance with Different Shapes and Positions
WU Jiansheng 1,2, JIANG Wancong 1,2
1. College of Ocean and Earth Science, Tongji University, Shanghai 200092, China;
2. State Key Laboratory of Marine Geology, Tongji University, Shanghai 200092, China
Abstract: Based on surface magnetic element integral theory, a fast high-precision magnetic forward method was proposed to solve the magnetic fast forward-modeling of magnetized unexploded ordnance with different shapes and positions, by using triangular surface element method and spline interpolation method.The method realized the fast interpolation forward modeling of finite-length cylinder and variable-diameter finite-length cylinder.It is shown that the error of the method is relatively small compared with that of non-interpolation method.
Key words: magnetic anomaly    forward modeling    surface element    spline interpolation

1 基于三角面元法的磁异常正演方法 1.1 三角面元法正演原理

 ${\sigma _{\rm{m}}} = {\mu _0}\left( {{\mathit{\boldsymbol{M}}_1} - {\mathit{\boldsymbol{M}}_2}} \right) \cdot \mathit{\boldsymbol{n}}$ (1)

 ${U_{\rm{m}}} = \frac{1}{{4{\rm{ \mathsf{ π} }}{\mu _0}}}\int_\mathit{\boldsymbol{V}} {\frac{{{\rho _{\rm{m}}}{\rm{d}}\mathit{\boldsymbol{V}}}}{\mathit{\boldsymbol{r}}}} + \frac{1}{{4{\rm{ \mathsf{ π} }}{\mu _0}}}\int_\mathit{\boldsymbol{S}} {\frac{{{\sigma _{\rm{m}}}{\rm{d}}\mathit{\boldsymbol{S}}}}{\mathit{\boldsymbol{r}}}}$ (2)

 $\mathit{\boldsymbol{T}} = - {\mu _0}\frac{{\partial {U_{\rm{m}}}}}{{\partial \mathit{\boldsymbol{t}}}} = \frac{{{\mu _0}}}{{4{\rm{ \mathsf{ π} }}}}\int_\mathit{\boldsymbol{S}} {{\sigma _{\rm{m}}}\frac{{\mathit{\boldsymbol{t}} \cdot \mathit{\boldsymbol{r}}}}{{{\mathit{\boldsymbol{r}}^3}}}{\rm{d}}\mathit{\boldsymbol{S}}}$ (3)

Um分别对xyz求导即可求得磁场水平x方向分量Hax, 水平y方向分量Hay和垂直方向分量Za的积分式, 如下所示:

 $\begin{array}{l} {H_{{\rm{a}}x}} = \frac{{{\mu _0}}}{{4{\rm{ \mathsf{ π} }}}}\int_\mathit{\boldsymbol{S}} {\frac{{{\mathit{\boldsymbol{M}}_{\rm{n}}}\cos \left( {\mathit{\boldsymbol{r}},x} \right)}}{{{\mathit{\boldsymbol{r}}^2}}}{\rm{d}}\mathit{\boldsymbol{S}}} = \frac{{{\mu _0}}}{{4{\rm{ \mathsf{ π} }}}}\int_\mathit{\boldsymbol{S}} {\frac{{{\mathit{\boldsymbol{M}}_{\rm{n}}}\left( {x - \xi } \right)}}{{{\mathit{\boldsymbol{r}}^3}}}{\rm{d}}\mathit{\boldsymbol{S}}} \\ {H_{{\rm{a}}y}} = \frac{{{\mu _0}}}{{4{\rm{ \mathsf{ π} }}}}\int_\mathit{\boldsymbol{S}} {\frac{{{\mathit{\boldsymbol{M}}_{\rm{n}}}\cos \left( {\mathit{\boldsymbol{r}},y} \right)}}{{{\mathit{\boldsymbol{r}}^2}}}{\rm{d}}\mathit{\boldsymbol{S}}} = \frac{{{\mu _0}}}{{4{\rm{ \mathsf{ π} }}}}\int_\mathit{\boldsymbol{S}} {\frac{{{\mathit{\boldsymbol{M}}_{\rm{n}}}\left( {y - \eta } \right)}}{{{\mathit{\boldsymbol{r}}^3}}}{\rm{d}}\mathit{\boldsymbol{S}}} \\ {Z_{\rm{a}}} = \frac{{{\mu _0}}}{{4{\rm{ \mathsf{ π} }}}}\int_\mathit{\boldsymbol{S}} {\frac{{{\mathit{\boldsymbol{M}}_{\rm{n}}}\cos \left( {\mathit{\boldsymbol{r}},z} \right)}}{{{\mathit{\boldsymbol{r}}^2}}}{\rm{d}}\mathit{\boldsymbol{S}}} = \frac{{{\mu _0}}}{{4{\rm{ \mathsf{ π} }}}}\int_\mathit{\boldsymbol{S}} {\frac{{{\mathit{\boldsymbol{M}}_{\rm{n}}}\left( {z - \zeta } \right)}}{{{\mathit{\boldsymbol{r}}^3}}}{\rm{d}}\mathit{\boldsymbol{S}}} \end{array}$ (4)

 $\Delta T = {H_{{\rm{a}}x}}\cos I\cos D + {H_{{\rm{a}}y}}\cos I\sin D + {Z_{\rm{a}}}\sin I$ (5)

1.2 以球体模型为基础的正演实验

1.2.1 球体模型的剖分方法

 $\begin{array}{l} \left[ {\begin{array}{*{20}{c}} {{D_{ij}}}\\ {{D_{i + 1,j}}}\\ {{D_{i + 1,j + 1}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{x_{ij}}}&{{y_{ij}}}&{{z_{ij}}}\\ {{x_{i + 1,j}}}&{{y_{i + 1,j}}}&{{z_{i + 1,j}}}\\ {{x_{i + 1,j + 1}}}&{{y_{i + 1,j + 1}}}&{{z_{i + 1,j + 1}}} \end{array}} \right]\\ \left[ {\begin{array}{*{20}{c}} {{D_{ij}}}\\ {{D_{i + 1,j + 1}}}\\ {{D_{i,j + 1}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{x_{ij}}}&{{y_{ij}}}&{{z_{ij}}}\\ {{x_{i + 1,j + 1}}}&{{y_{i + 1,j + 1}}}&{{z_{i + 1,j + 1}}}\\ {{x_{i,j + 1}}}&{{y_{i,j + 1}}}&{{z_{i,j + 1}}} \end{array}} \right] \end{array}$ (6)
 图 1 球体三角面元剖分示意图 Fig.1 Schematic diagram of sphere division in triangular surface element

 图 2 球体三角面元剖分投影示意图 Fig.2 Projection schematic diagram of sphere division in triangular surface element

 $\begin{array}{*{20}{c}} {{x_{ij}} = \sqrt {{R^2} - {{\left( {\left( {i - 1} \right)\Delta d} \right)}^2}} \cos \left( {\left( {j - 1} \right)\Delta \theta } \right)}\\ {{y_{ij}} = \sqrt {{R^2} - {{\left( {\left( {i - 1} \right)\Delta d} \right)}^2}} \sin \left( {\left( {j - 1} \right)\Delta \theta } \right)}\\ {{z_{ij}} = h + \left( {i - 1} \right)\Delta d} \end{array}$ (7)

 $\begin{array}{*{20}{c}} {{x_{ij}} = \frac{{\sqrt {{R^2} - {{\left( {\left( {i - 1} \right)\Delta d} \right)}^2}} }}{{\cos \left( {2{\rm{ \mathsf{ π} }}/2P} \right)}}\cos \left( {{\rm{ \mathsf{ π} }}/P + \left( {j - 1} \right)\Delta \theta } \right)}\\ {{y_{ij}} = \frac{{\sqrt {{R^2} - {{\left( {\left( {i - 1} \right)\Delta d} \right)}^2}} }}{{\cos \left( {2{\rm{ \mathsf{ π} }}/2P} \right)}}\sin \left( {{\rm{ \mathsf{ π} }}/P + \left( {j - 1} \right)\Delta \theta } \right)}\\ {{z_{ij}} = h + \left( {i - 1} \right)\Delta d} \end{array}$ (8)

1.2.2 球体模型的正演结果

 图 3 不同方法球体磁异常剖面曲线对比 Fig.3 Comparison of section curves of spherical magnetic anomaly among different methods
2 不同姿态UXO模型的正演方法

UXO大体上可以分为4类:集束弹药、普通陆战武器弹药、航空炸弹、地雷[10].UXO外壳多为铁磁性物质, 磁性参数与周围环境的差别较大, 在一定条件下能够被磁法勘探仪器所识别.然而, UXO外形差别较大.航空炸弹一般体形较大, 长度较长, 外形与圆柱体比较接近, 因此可使用均匀有限长的三度体圆柱体模拟航空炸弹; 其他几类弹药体形较小, 体长较短, 长度和埋深的比值很小, 理论和实践已经证明这几类弹药产生的磁场近似于球体所产生的磁场.常见的航空炸弹的外形往往不是标准的有限长圆柱体, 而是变径有限长圆柱体, 即沿长轴方向半径不同, 用单一有限长圆柱体的模型无法很好地近似, 但可以用不同半径的圆柱体元进行组合.

2.1 有限长圆柱体模型的剖分方法

 图 4 圆柱体三角面元剖分示意图 Fig.4 Schematic diagram of cylinder division in triangular surface element

2.2 全空间任意姿态有限长圆柱体磁异常计算方法

 $\begin{array}{l} {\mathit{\boldsymbol{T}}_\beta } = \left[ {\begin{array}{*{20}{c}} {\cos \beta }&{ - \sin \beta }&0\\ {\sin \beta }&{\cos \beta }&0\\ 0&0&1 \end{array}} \right]\\ {\mathit{\boldsymbol{T}}_\alpha } = \left[ {\begin{array}{*{20}{c}} {\cos \alpha }&0&{\sin \alpha }\\ 0&1&0\\ {\sin \alpha }&0&{\cos \alpha } \end{array}} \right] \end{array}$ (9)

 $\begin{array}{l} {\left( {{{x'}_{ij}},{{y'}_{ij}},{{z'}_{ij}}} \right)^{\rm{T}}} = \\ {\mathit{\boldsymbol{T}}_\beta } \times {\mathit{\boldsymbol{T}}_\alpha } \times {\left( {{x_{ij}},{y_{ij}},{z_{ij}}} \right)^{\rm{T}}} + {\left( {{x_1},{y_1},{z_1}} \right)^{\rm{T}}} \end{array}$ (10)

2.3 有限长圆柱体模型的正演结果

 图 5 不同走向和倾角下有限长圆柱体的形态及其产生的磁场 Fig.5 Shape of finite-length cylinder in different strikes and dip angles and the corresponding magnetic field

 图 6 不同方法垂直有限长圆柱体磁异常剖面曲线对比 Fig.6 Comparison of section curves of vertical finite-length cylindrical magnetic anomaly among different methods

3 基于三次样条的UXO磁异常快速正演方法

3.1 三次样条插值方法基本原理

Lagrange插值、Newton插值及Hermite插值等高次插值中会出现龙格(Runge)现象.分段插值虽然可以避免龙格现象, 但是分段线性插值在节点处不光滑。分段Hermite插值仅仅要求插值函数在节点上等于被插函数的一阶导数, 并且导数值往往并不已知.三次样条插值先由函数值确定导数值, 再由分段Hermite插值解决问题, 因而避免了龙格现象, 并且插值曲线光滑.

(1) S(l)∈C2[a, b], 即具有连续的一阶、二阶导数, Sk, v(lk)=Sk-1, v(lk), v=0, 1, 2, k=0, 1, …, n-1.

(2) S(l)在每一个小区间[lk, lk+1]上是次数≤3的多项式, 则称S(l)为关于剖分Δ的一个三次样条函数.

(3) 设给定函数表(lq, f(lq)), q=0, 1, …, n, 若三次样条函数S(l)还满足以下插值条件:

 $S\left( {{l_q}} \right) = f\left( {{l_q}} \right),\;\;\;q = 0,1, \cdots ,n$

(1) 第一种边界条件.S′(l0)、S′(ln)已知, 即S′(l0)=f′(l0), S′(ln)=f′(ln).

(2) 第二种边界条件.S″(l0)、S″(ln)已知, 即S″(l0)=f″(l0), S″(ln)=f″(ln), 若S″(l0)=S″(ln)=0, 则称为自然边界条件.

(3) 第三种边界条件(周期边界条件).f(l)为周期函数(此时有f(l0)=f(ln)), 要求S(l)亦是周期函数, 周期为(b-a), 即取Sν(l0+)=Sν(ln-), v=0, 1, 2, 称S(l)为周期样条函数, 此时不用给出边界值.

(4) 非扭结边界.使两端点处多项式三次项系数与这两端点的邻近点的三次项系数相等, 此时不用给出边界值.

3.2 不同模型的正演结果

 图 7 水平圆柱体磁异常三次样条插值结果 Fig.7 Results of horizontal cylindrical magnetic anomaly calculated with spline interpolation

 图 8 三角面元构造的变径纺锤形UXO模型 Fig.8 Fusiform UXO model with variable diameters made by triangular surface element

 图 9 变径有限长圆柱体磁异常三次样条插值结果 Fig.9 Results of variable diameter cylindrical magnetic anomaly calculated with spline interpolation
4 结语

 [1] ZALEVSKY Z, BREGMAN Y, SALOMONSKI N, et al. Resolution enhanced magnetic sensing system for wide coverage real time UXO detection[J]. Journal of Applied Geophysics, 2012, 84(4): 70 [2] CHURCHILL K M, LINK C, YOUMANS C C. A comparison of the finite-element method and analytical method for modeling unexploded ordnance using magnetometer[J]. IEEE Transactions on Geoscience and Remote Sensing, 2012, 50(7): 2720 DOI:10.1109/TGRS.2011.2174796 [3] SINEX D B. Advancing the state of the art of UXO discrimination for total-field magnetic data[D]. Golden: Colorado School of Mines, 2006. [4] SANCHEZ V, LI Y G, NABIGHIAN M N, et al. Numerical modeling of higher order magnetic moments in UXO discrimination[J]. IEEE Transactions on Geoscience and Remote Sensing, 2008, 46(9): 2568 DOI:10.1109/TGRS.2008.918090 [5] 葛健, 陆承达, 董浩斌, 等. 基于Overhauser传感器的近地表UXO磁梯度法探测技术[J]. 仪器仪表学报, 2015, 36(5): 961 GE Jian, LU Chengda, DONG Haobin, et al. The detection technology of near-surface UXO based on magnetic gradient method and Overhauser sensor[J]. Chinese Journal of Scientific Instrument, 2015, 36(5): 961 [6] 郭志馗, 陈超, 陶春辉, 等. 有限长圆柱体磁异常场全空间正演方法[J]. 地球物理学报, 2017, 60(4): 1557 GUO Zhikui, CHEN Chao, TAO Chunhui, et al. Forward modeling of magnetic anomaly of finite-length cylinder in 3D space[J]. Chinese Journal of Geophysics, 2017, 60(4): 1557 DOI:10.6038/cjg20170428 [7] 朱慧慧, 刘得军, 冯硕, 等. 基于磁偶极子构造法几何建模的铁管磁异常正演[J]. 计算机测量与控制, 2017, 25(3): 201 ZHU Huihui, LIU Dejun, FENG Shuo, et al. Magnetic dipole reconstruction geometry modeling for underground ferromagnetic pipe magnetic abnormal detection[J]. Computer Measurement and Control, 2017, 25(3): 201 [8] 吴文鹂, 管志宁, 高艳芳, 等. 重磁异常数据三维人机联作模拟[J]. 物探化探计算技术, 2005, 27(3): 227 WU Wenli, GUAN Zhining, GAO Yanfang, et al. Interactive man/computer modeling 3D body of gravity and magnetic anomaly data[J]. Computing Techniques for Geophysical and Geochemical Exploration, 2005, 27(3): 227 [9] 管志宁. 地磁场与磁力勘探[M]. 北京: 地质出版社, 2005 GUAN Zhining. Geomagnetic field and magnetic exploration[M]. Beijing: Geological Publishing House, 2005 [10] 李小康, 杨磊, 刘磊. 未爆弹药问题及地球物理解决方案综述[J]. 中国矿业, 2010, 19: 187 LI Xiaokang, YANG Lei, LIU Lei. Overview on unexploded ordnance problem and a solution:the geophysical scheme[J]. China Mining Magzine, 2010, 19: 187 [11] 关冶, 陆金甫. 数值方法[M]. 北京: 清华大学出版社, 2011 GUAN Ye, LU Jinfu. Numerical method[M]. Beijing: Tsinghua University House, 2011