一种多辐射场源瞬变电磁法多分量全域视电阻率定义方法

摘要

本发明公开了一种多辐射场源瞬变电磁法多分量全域视电阻率定义方法，基于多辐射场源地空瞬变电磁法理论体系，建立了相应的多分量全域视电阻率定义方法。系统分析了利用磁场强度进行多辐射场源地空系统全域视电阻率定义的优点，针对磁场强度的各个分量，提出了各自的全域视电阻率算法，实现了多分量、全时域、全空域视电阻率计算，并分析了偏移距对全域视电阻率的影响。多辐射源瞬变电磁地空系统不仅可以加强信号强度，削弱随机干扰，通过调整源之间的相对位置，还可更好的分辨地下异常体的位置。通过对所设计模型的处理，证实了多辐射场源地空系统多分量全域视电阻率算法的有效性，也验证了多辐射场源地空系统加大勘探深度、提高分辨率的优越性。

说明书

技术领域

本发明涉及一种定义方法，涉及一种多辐射场源瞬变电磁法多分量全域视电阻率定义方 法。

背景技术

地空瞬变电磁法已成功应用在地热调查、火山结构、地下水盐渍化及地下水监测等领域， 由于电性源地空系统仍属时间域电磁法的范畴，且时间域电磁法特别适合深部、隐伏矿床勘 查，因此该法在矿产勘探中也具备较大的发展潜力。目前我国的矿业形势十分严峻，一方面 老矿山的储量日趋枯竭，另一方面新发现的矿床又越来越少，因此国家明确提出了危机矿山 和第二深度空间探矿的专项计划，这两项专项计划针对的都是大深度问题，都要求矿产勘查 向深部发展，目前已知的电性源地空系统的最大探测深度可达800m，若采用多辐射场源， 还可进一步提高该法的探测深度，将瞬变电磁法引入到深部探测。另一方面矿产勘查工作迫 切需要向地形复杂的山区等尚未被详细探查的区域发展，由于地空瞬变电磁法在空中进行数 据采集，可实现阵列式观测及三维信息采集，相比较地面瞬变电磁法在这些地区具有更大的 可操作性。

到目前为止，电性源地空系统仍采用单辐射源工作模式，数据解释工作大都仍基于地面 LOTEM，李貅等结合TEM三维拟地震解释技术，建立了地空逆合成孔径成像的三维解释理 论体系。关于电性源地空系统的大部分研究工作都集中在数据校正以及设备研发方面，校正 工作包括有运动校正、消除当地噪声，基线漂移校正、去噪等，Verma等分析了飞行高度、 偏移距等参数对电性源地空系统响应特征的影响。目前得到成功应用的系统有两套，一是 Mogi主持研发的直升机平台的GREATEM系统，它采用磁力计直接测量磁感应强度(B(t))， 该系统已成功应用于地热调查及火山结构调查；另一套是吉林大学自主研制的飞艇搭载的电 性源发射系统，该系统采集的是电压响应，它与磁感应强度对时间的变化率成 正比，该系统在地下水盐渍化及地下水监测方面取得了良好的应用效果。由于对多辐射源地 空瞬变电磁法的研究工作尚未展开，解决了全时域、全空域的多分量全域视电阻率定义问题， 不仅可以提供一种定性的直观的初步解释方法，为三维解释技术以及反演提供更丰富的信 息，还可利用各个分量各自的特点更好的反映地下电性信息。

当遇到良导目标体时，由于会很快衰减到噪声水平之下，这会限制电性源地空 系统的探测深度；而B(t)作为关于t的积分，虽然对于浅部的分辨能力要相对弱一 些，但是由于B(t)衰减慢，在晚期时幅值仍较大，因此在晚期仍有较高的分辨率。另一方面， 通过计算分析我们发现，B(t)关于电阻率参数的函数形态要相对简单的多，而不仅 是多解函数，在某些时刻，甚至是无合适解、无解的，此外它还会带来诸如“overshoot”和 “undershoot”的假异常干扰，显然相对于利用B(t)进行多辐射源地空系统数据解 释具有其独特的优势。但是从目前来说，由于受到制作工艺和使用范围等方面的限制，在地 空系统设备中使用磁力计仍需解决航空电磁法中存在的诸如高温、带宽等问题。

发明内容

本发明的目的在于克服上述技术存在的缺陷，提供一种多辐射场源瞬变电磁法多分量全 域视电阻率定义方法，采用对进行积分求取B(t)的方法进行数据解释。鉴于多辐射 场源本身存在的诸多优势，如加强采集信号强度、减弱随机噪声、更全面反映地下异常体位 置信息等，若再结合B(t)进行多分量数据解释，必会更大提高地空系统的探测深度以及分辨 率。理论模型计算证实了利用B(t)进行多辐射源地空系统多分量全域视电阻率定义的有效 性，也验证了多辐射场源地空系统提高分辨率的优越性。

其具体技术方案为：

一种多辐射场源瞬变电磁法多分量全域视电阻率定义方法，包括以下步骤：

1)磁场强度的x和z分量

将磁场强度的x和z分量记为B_p(ρ，C，t)(p＝x，z)，其中C表示测点在空中的位置坐标参数。 给定初值在的邻域内对B_p(ρ，C，t)(p＝x，z)进行泰勒展开

$\left(\begin{array}{c}{B}_p(\rho ,C,t)=B_p({\rho }_{\tau }^{\left(0\right)},C,t)+{B_p}^{\prime }({\rho }_{\tau }^{\left(0\right)},C,t)(\rho -{\rho }_{\tau }^{\left(0\right)})+\frac{{B_p}^{\prime \prime }({\rho }_{\tau }^{\left(0\right)},C,t)}{2!}{(\rho -{\rho }_{\tau }^{\left(0\right)})}^2+...\\ +\frac{{B_p}^{\left(n\right)}({\rho }_{\tau }^{\left(0\right)},C,t)}{n!}{(\rho -{\rho }_{\tau }^{\left(0\right)})}^n+R_n\left(\rho \right)(p=x,z)\end{array}\right)---\left(7\right)$

由于该邻域非常小，在该邻域内用一小段直线对曲线进行近似，即保留(7)式的前两 项，有

$B_p(\rho ,C,t)\approx B_p({\rho }_{\tau }^{\left(0\right)},C,t)+{B_p}^{\prime }({\rho }_{\tau }^{\left(0\right)},C,t)(\rho -{\rho }_{\tau }^{\left(0\right)})(p=x,z)---\left(8\right)$

对上式进行变换

$\rho =\frac{B_p(\rho ,C,t)-B_p({\rho }_{\tau }^{\left(0\right)},C,t)}{{B_p}^{\prime }({\rho }_{\tau }^{\left(0\right)},C,t)}+{\rho }_{\tau }^{\left(0\right)}(p=x,z)---\left(9\right)$

将(9)式写成迭代的形式

${\rho }_{\tau }^{(i+1)}\approx \Delta {\rho }_{\tau }^{\left(i\right)}+{\rho }_{\tau }^{\left(i\right)}(i=\mathrm{0,1,2},...)---\left(10\right)$

其中

${\mathrm{\Delta \rho }}_{\tau }^{\left(i\right)}=\frac{B_p(\rho ,C,t)-B_p({\rho }_{\tau }^{(i-1)},C,t)}{{B_p}^{\prime }({\rho }_{\tau }^{(i-1)},C,t)}(p=x,z)---\left(11\right)$

(10)式的迭代终止条件为

$\left|\frac{B_p(\rho ,C,t)-B_p({\rho }_{\tau }^{\left(i\right)},C,t)}{B_p(\rho ,C,t)}\right|<ϵ(p=x,z)---\left(12\right)$

ε是给定的迭代终止误差限，B_p(ρ，C，t)(p＝x，z)是测得的磁场分量；

2)磁场强度的y分量

$\left(\begin{array}{c}{{\rho }_{\tau }}^{\left(k\right)}\left(i\right)=-\frac{1}{6}\{{{\rho }_{\tau }}^{(k-1)}(i+2)+{{\rho }_{\tau }}^{\left(k\right)}(i-2)-4[{{\rho }_{\tau }}^{(k-1)}(i+1)+{{\rho }_{\tau }}^{\left(k\right)}(i-1)\left]\right\}(i=1,...,M)\\ \left(\begin{array}{c}{\rho }_{\tau }(i-1)=2{\rho }_{\tau }\left(i\right)-{\rho }_{\tau }(i+1)(i=1)\\ {\rho }_{\tau }(i+1)=2{\rho }_{\tau }\left(i\right)-{\rho }_{\tau }(i-1)(i=M)\\ {\rho }_{\tau }(i-2)={\rho }_{\tau }(i+2)-2[{\rho }_{\tau }\left(i\right)-{\rho }_{\tau }(i-1)](i=1)\\ {\rho }_{\tau }(i+2)={\rho }_{\tau }(i-2)+2[{\rho }_{\tau }\left(i\right)-{\rho }_{\tau }(i-1)](i=M)\end{array}\right)\end{array}\right)---\left(13\right)$

其中k表示插值的次数，M表示参与插值的总个数，一般在需要插值的区域左右两端各 取2个点参与计算。

与现有技术相比，本发明的有益效果为：

本发明关于电阻率参数的函数形态更为简单、晚期幅值较大、不会带来类似于 “overshoot”和“undershoot”的假异常、削弱高频随机噪声等；基于反函数思想提出了利 用磁场强度进行多辐射源地空瞬变电磁法多分量全域视电阻率定义的算法，实现了时间上不 分早晚、距离上不分远近的多分量视电阻率计算，可以完整、光滑、渐变地反映出地下的电 性信息变化；并分析了偏移距变化对全域视电阻率的影响；建立起了一套完整的多辐射源地 空瞬变电磁法全域视电阻率算法，为多辐射源地空瞬变电磁法提供了直观、基础、定性解释 方法。模型试验验证了利用磁场强度的多辐射场源地空系统多分量全域视电阻率算法的有效 性，也证明了多辐射场源地空瞬变电磁方法提高分辨率的优越性。

附图说明

图1是坐标系示意图，其中，图1(a)地表电偶极子坐标系示意图，图1(b)多辐射源剖 分俯视图；

图2是均匀半空间模型双源地空系统及测点坐标俯视图，其中，图2(a)源所在坐标系俯 视图，图2(b)测点坐标俯视图；

图3是不同偏移距下磁场响应在不同时刻随均匀半空间模型电阻率变 化，其中，图3(a)偏移距500m，图3(b)偏移距500m，图3(c)偏移 距500m，图3(d)偏移距5000m，图3(e)偏移距5000m，图 3(f)偏移距5000m，

图4是不同偏移距下磁场响应B_p(t)(p＝x，y，z)在不同时刻随均匀半空间模型电阻率变化， 其中图4(a)偏移距500m，B_x(t)，图4(b)偏移距500m，B_y(t)，图4(c)偏移距500m，B_z(t)， 图4(d)偏移距5000m，B_x(t)，图4(e)偏移距5000m，B_y(t)，图4(f)偏移距5000m，B_z(t)；

图5是添加信噪比为100、50、20、10高斯白噪声的积分求解的B_Z(t)与滤波解 对比图(模型和计算参数：电阻率ρ₁＝100Ω·m，ρ₂＝10Ω·m，ρ₃＝100Ω·m，层厚h₁＝20m，h₂＝10m， 源的长度1000m，电流大小1A，飞行高度100m，偏移距500m坐标(300m，400m，-100m))，其 中，图5(a)偏移距500m测点坐标示意图(b)添加噪声后求得的B_Z(t)与滤波解对比图；

图6是磁场强度x和z分量全域视电阻率定义方法流程图；

图7是磁场强度y分量全域视电阻率定义方法流程图；

图8是进行最小曲率补足前后的全域视电阻率对比图(模型和计算参数：电阻率 ρ₁＝100Ω·m，ρ₂＝10Ω·m，ρ₃＝100Ω·m，层厚h₁＝20m，h₂＝10m，源1、2的长度均为1000m，两源 夹角300，电流大小均为50A，飞行高度100m，偏移距500m坐标(300m，400m，-100m))，其 中，图8(a)电性源及测点坐标俯视图，图8(b)补足前，图8(c)补足后；

图9是偏移距500m时用B_Z(t)和定义的视电阻率曲线(模型和计算参数：电阻率 ρ₁＝100Ω·m，ρ₂＝10Ω·m，ρ₃＝100Ω·m，层厚h₁＝20m，h₂＝10m，源1、2的长度均为1000m，两源 夹角300，电流大小均为50A，飞行高度100m，偏移距500m坐标(300m，400m，-100m))，其 中图9(a)电性源及测点坐标俯视图，图9(b)B_Z(t)所得的视电阻率曲线，图9(c)所得的视电阻率曲线；

图10是偏移距500m时用B_p(t)(p＝x，y，z)定义的全域视电阻率曲线(模型和计算参数：电 阻率ρ₁＝100Ω·m，ρ₂＝10Ω·m，ρ₃＝100Ω·m，层厚h₁＝20m，h₂＝10m，源1、2的长度均为1000m， 两源夹角300，电流大小均为50A，飞行高度100m，偏移距500m坐标(300m，400m，-100m)) 其中，图10(a)电性源及测点坐标俯视图，图10(b)B_x(t)所得的视电阻率曲线，图10 (c)B_y(t)所得的视电阻率曲线，图10(d)B_z(t)所得的视电阻率曲线；

图11偏移距变化对多辐射源地空系统多分量全域视电阻率曲线的影响(模型和计算参 数：电阻率ρ₁＝500Ω·m，ρ₂＝100Ω·m，ρ₃＝500Ω·m，层厚h₁＝800m，h₂＝800m，源1、2的长度均 为1000m，两源夹角300，电流大小均为50A，飞行高度100m，偏移距500m、1500m、3000m、 5000m对应的位置坐标分别为(300m，400m，-100m)、(1060m，1060m，-100m)、 (2121m，2121m，-100m)和(3000m，4000m，-100m))，其中，图11(a)电性源及测点坐标俯视图， 图11(b)偏移距对B_x(t)视电阻率曲线的影响，图11(c)偏移距对B_y(t)视电阻率曲线的影响，图 11(d)偏移距对B_z(t)视电阻率曲线的影响；

图12是模型示意图，其中图12(a)低阻体模型示意图，图12(b)低阻体尺寸参数示意 图；

图13是单辐射场源模型俯视图及偏移距10m(Line24)的视电阻率剖面图，其中，图 13(a)模型俯视图，图13(b)由B_x(t)得到的Line24视电阻率剖面图，图13(c)由B_y(t)得到的 Line24视电阻率剖面图，图13(d)由z(t)得到的Line24视电阻率剖面图；

图14是多辐射场源模型俯视图及偏移距10m(Line24)的视电阻率剖面图，其中，图 14(a)模型俯视图，图14(b)由B_x(t)得到的Line24视电阻率剖面图，图14(c)由B_y(t)得到 的Line24视电阻率剖面图，图14(d)由B_z(t)得到的Line24视电阻率剖面图。

具体实施方式

下面结合具体附图和实施例对本发明的技术方案作进一步详细地说明。

1理论与分析

1.1多辐射场源地空系统理论

图1(a)所示为均匀大地地表电偶极子所在坐标系示意图，AB表示位于地表的电偶极 子，电偶极子的中心即为坐标系原点O，z轴指向大地，整个坐标系满足右手坐标系准则， M表示空中任意位置的测点，它与O之间为R，M在地表的投影为P，P与O之间为r，即 本文中提及的偏移距。图1(b)所示为多辐射场源的坐标系俯视图，A_iB_i表示第i个电性源， 由于实际应用中铺设的电性源可能长达数公里，因此我们采用了对电性源进行剖分，采用有 限求和的方法进行计算，相应的多辐射源均匀大地地表电性源在空中产生的频率域总磁感应 强度B_p(ω)(p＝x，y，z)如(1)式所示。

其中P_Eij＝Ids_ij，I是电流大小，ds_ij表示第i个电性源剖分的第j个偶极子的长度， k₁是介质的波数，和r_ij分别表示对应第i个电性源剖分 的第j个偶极子的夹角和偏移距，n表示电性源的总个数，n_i表示第i个源的剖分个数。

当使用阶跃波作为发射波形时，时间域的场可以由频率域的场通过正、余弦变换得到， 见(2)式。

$\left(\begin{array}{c}{B}_p\left(t\right)=\frac{2}{\pi }{\int }_0^{\infty }\frac{\mathrm{Im}{B}_p\left(\omega \right)}{\omega }\cos \left(\mathrm{\omega t}\right)\mathrm{d\omega }(p=x,y,z)\\ \partial B_p\left(t\right)/\partial t=\frac{2}{\pi }{\int }_0^{\infty }\mathrm{Im}{B}_p\left(\omega \right)\sin \left(\mathrm{\omega t}\right)\mathrm{d\omega }(p=x,y,z)\end{array}\right)---\left(2\right)$

1.2多辐射源地空系统场特征分析

视电阻率定义是从均匀半空间的电阻率入手的，如果在均匀半空间模型中，相同的位置 参数下得到的响应与测得的响应相等的话，就把均匀半空间模型的电阻率视为该时间道的视 电阻率。但是由于多辐射源地空系统响应与电阻率参数之间存在复杂的函数关系，我们无法 直接得到一个用场表示的显式关系式来表达视电阻率。为了使分析不失一般性，我们采用了 两个相交的源进行分析，源所在的坐标系以及测点坐标的俯视图见图2(a)和(b)。

计算采用的参数如下：源1和源2的长度均为1000m，两个源之间的夹角为30°，电流 大小均为50A，飞行高度100m，偏移距500和5000m在X₁Y₁坐标系中对应的坐标分别为 (300m，400m，-100m)和(3000m，4000m，-100m)。图3和4所示分别为不同偏移距下对应图2的 多辐射源地空系统响应和B_p(t)(p＝x，y，z)随均匀半空间模型电阻率的变 化。从图3(a)～(c)中可以看出，在电阻率范围(0.001～10000Ω·m)内，小偏移距(500m)情况下， 基本在任意时刻都是关于均匀半空间模型电阻率的双值函数，在大偏移 距(5000m)情况下，某些时刻可视为单值函数，而在另外一些时刻仍是双值函数；从图4中 可以看出，B_x(t)和B_z(t)无论是在大偏移距还是小偏移距情况下，基本在任何时刻都是关于半 空间电阻率参数的单调函数，而B_y(t)的情况与类似，尽管如此，很多学 者已经证实，使用做视电阻率定义会带来类似于“overshoot”和“undershoot”的假异常， 而使用B(t)得到的视电阻率曲线能够平滑渐变的反映出地下的电性信息。因此对于多辐射源 地空系统瞬变电磁法，相比较使用B(t)做视电阻率定义具有更大的优势。

1.3 B(t)和的比较

对于时间域电磁方法，目前大部分仪器都是采用感应线圈直接测量但遇到良 导目标体时，会很快衰减到噪声水平之下，使得时间域电磁法的探测深度受限，而B(t) 由于衰减较慢，而且且作为对时间的积分，幅值也较大，这就使得B(t)在晚期仍能保 持较高的分辨率，更适合深部探测，利用B(t)做数据解释，由于积分过程相当于一个低通滤 波器，这种做法还可削弱随机的高频噪声。此外，虽然现在有些电性源地空系统瞬变电磁法 设备，如GREATEM系统，可以采用磁力计直接对B(t)进行测量，但跟航空电磁法类似，仍 存在较多的问题需要处理，如高温、带宽等。事实上，从目前来说，对于实测的我们可以采用积分的办法求取B(t)。(3)式所示为利用采用倒序积分的方法求取B(t)， t_i和t_i+1分别为第i个时间道积分窗口的上下限，N为最大测道数。

$B_p\left(t_i\right)={\int }_{t_i}^{t_{i+1}}\frac{\partial B_p\left(t\right)}{\partial t}\mathrm{dt}+{\int }_{t_{i+2}}^{t_{i+3}}\frac{\partial B_p\left(t\right)}{\partial t}\mathrm{dt}+···+{\int }_{t_N}^{t_{N+1}}\frac{\partial B_p\left(t\right)}{\partial t}\mathrm{dt}+{\int }_{t_{N+1}}^{\infty }\frac{\partial B_p\left(t\right)}{\partial t}\mathrm{dt}(p=x,y,z)---\left(3\right)$

令

$C_p={\int }_{t_{N+1}}^{\infty }\frac{\partial B_p\left(t\right)}{\partial t}\mathrm{dt}(p=x,y,z)---\left(4\right)$

将(3)式写成求和的形式，即

$B_p\left(t_i\right)=\underset{j=i}{\overset{N}{\Sigma }}\Delta t_j\frac{\partial B_p\left(t_j\right)}{\partial t}+C_p(p=x,y,z)---\left(5\right)$

其中Δt_j＝t_j+1-t_j，关于积分常数C的求取，由于在晚期，和t在双对数坐 标系下具有线性关系，且斜率小于-1，因此利用了最后几个时间道的采用线性回 归的方法求出$\ln \left(\frac{\partial B_p\left(t\right)}{\partial t}\right)=a_p\ln \left(t\right)+b_p(p=x,y,z)$中的a_p和b_p，则

$C_p=-e^{b_p}\frac{1}{a+1}{t}^{a_p+1}(p=x,y,z)---\left(6\right)$

我们采用正演得到的数据对上述积分方法进行了验证，计算参数如下：三层模型的电阻 率分别为ρ₁＝100Ω·m，ρ₂＝10Ω·m，ρ₂＝100Ω·m，层厚h₁＝20m，h₂＝10m，单源的长度为 1000m，电流大小1A，飞行高度100m，偏移距500m的测点对应的位置坐标为 (300m，400m，-100m)。图5(a)所示为偏移距500m测点坐标示意图，(b)为对添加 信噪比为100、50、20、10的高斯白噪声后积分得到的B_z(t)与滤波解对比图，从图中可以看 出，随着信噪比的下降，积分结果在晚期几道的偏差较大，这是因为我们并未对添加噪声后 的做圆滑处理，导致(6)式中对C_p的求取不够精确，而由于早晚期数据的量级相 差较大，所以对早期数据的影响较小。此外由于是积分过程，B_z(t)对高频随机噪声的压制效 果明显，所以我们仍能得到光滑的B_z(t)曲线。

2多分量多辐射源全域视电阻率定义方法

这里所说的全域视电阻率定义方法，是在视电阻率计算过程中同时考虑到位置坐标、时 间等各个参数，实现时间上不分早晚、距离上不分远近的视电阻率定义。在第2节中已经分 析过，B_p(t)(p＝x，z)可被视为关于均匀半空间电阻率参数的单调函数，其他分量的情况要复 杂一些，因此将其分成两部分进行讨论。

2.1磁场强度的x和z分量

将磁场强度的x和z分量记为B_p(ρ，C，t)(p＝x，z)，其中C表示测点在空中的位置坐标参数。 给定初值在的邻域内对B_p(ρ，C，t)(p＝x，z)进行泰勒展开

$\left(\begin{array}{c}{B}_p(\rho ,C,t)=B_p({\rho }_{\tau }^{\left(0\right)},C,t)+{B_p}^{\prime }({\rho }_{\tau }^{\left(0\right)},C,t)(\rho -{\rho }_{\tau }^{\left(0\right)})+\frac{{B_p}^{\prime \prime }({\rho }_{\tau }^{\left(0\right)},C,t)}{2!}{(\rho -{\rho }_{\tau }^{\left(0\right)})}^2+...\\ +\frac{{B_p}^{\left(n\right)}({\rho }_{\tau }^{\left(0\right)},C,t)}{n!}{(\rho -{\rho }_{\tau }^{\left(0\right)})}^n+R_n\left(\rho \right)(p=x,z)\end{array}\right)---\left(7\right)$

由于该邻域非常小，在该邻域内用一小段直线对曲线进行近似，即保留(7)式的前两 项，有

$B_p(\rho ,C,t)\approx B_p({\rho }_{\tau }^{\left(0\right)},C,t)+{B_p}^{\prime }({\rho }_{\tau }^{\left(0\right)},C,t)(\rho -{\rho }_{\tau }^{\left(0\right)})(p=x,z)---\left(8\right)$

对上式进行变换

$\rho =\frac{B_p(\rho ,C,t)-B_p({\rho }_{\tau }^{\left(0\right)},C,t)}{{B_p}^{\prime }({\rho }_{\tau }^{\left(0\right)},C,t)}+{\rho }_{\tau }^{\left(0\right)}(p=x,z)---\left(9\right)$

将(9)式写成迭代的形式

${\rho }_{\tau }^{(i+1)}\approx \Delta {\rho }_{\tau }^{\left(i\right)}+{\rho }_{\tau }^{\left(i\right)}(i=\mathrm{0,1,2},...)---\left(10\right)$

其中

${\mathrm{\Delta \rho }}_{\tau }^{\left(i\right)}=\frac{B_p(\rho ,C,t)-B_p({\rho }_{\tau }^{(i-1)},C,t)}{{B_p}^{\prime }({\rho }_{\tau }^{(i-1)},C,t)}(p=x,z)---\left(11\right)$

(10)式的迭代终止条件为

$\left|\frac{B_p(\rho ,C,t)-B_p({\rho }_{\tau }^{\left(i\right)},C,t)}{B_p(\rho ,C,t)}\right|<ϵ(p=x,z)---\left(12\right)$

ε是给定的迭代终止误差限，B_p(ρ，C，t)(p＝x，z)是测得的磁场分量，该迭代算法的流程图 见图6。

2.2磁场强度的y分量

大量的模型计算表明，对于B_y(t)，视电阻率定义问题可大致分为三种情况：单解问题、 双解问题和无解问题。在第2节中已经分析过，在早期和晚期时间道，B_y(t)视电阻率的求解 是一个单解问题，这种情况与磁场强度的x和y分量类似。对于双解问题，首先需要找到极 值点所在位置，然后再在极值点两侧分别求取这两个解。进一步的研究表明，如果想从B_y(t) 数据中得到光滑渐变完整的视电阻率曲线，还需要将双解问题进一步分解为以下三种情况： 当两个解相差较小时(如10％)，对这两个解取平均作为该时间道的解；当两个解相差较大 时(如300％)，选择与前一个时间道的解更接近的那个解作为该时间道的解；其他情况都视 为无合适解问题。当层状大地产生的磁场比均匀大地的磁场大时，就会出现无解问题，这种 情况在LOTEM中也存在。对于上述分析中提及的无合适解和无解问题，为了得到完整的视 电阻率曲线，本文采用最小曲率插值方法对空缺的数据进行补足，该法能够保证曲线按最小 曲率变化且保持曲线光滑，无约束点原位最小曲率差分迭代算法见(13)式。利用B_y(t)进行 视电阻率定义的流程图见图7。

$\left(\begin{array}{c}{{\rho }_{\tau }}^{\left(k\right)}\left(i\right)=-\frac{1}{6}\{{{\rho }_{\tau }}^{(k-1)}(i+2)+{{\rho }_{\tau }}^{\left(k\right)}(i-2)-4[{{\rho }_{\tau }}^{(k-1)}(i+1)+{{\rho }_{\tau }}^{\left(k\right)}(i-1)\left]\right\}(i=1,...,M)\\ \left(\begin{array}{c}{\rho }_{\tau }(i-1)=2{\rho }_{\tau }\left(i\right)-{\rho }_{\tau }(i+1)(i=1)\\ {\rho }_{\tau }(i+1)=2{\rho }_{\tau }\left(i\right)-{\rho }_{\tau }(i-1)(i=M)\\ {\rho }_{\tau }(i-2)={\rho }_{\tau }(i+2)-2[{\rho }_{\tau }\left(i\right)-{\rho }_{\tau }(i-1)](i=1)\\ {\rho }_{\tau }(i+2)={\rho }_{\tau }(i-2)+2[{\rho }_{\tau }\left(i\right)-{\rho }_{\tau }(i-1)](i=M)\end{array}\right)\end{array}\right)---\left(13\right)$

其中k表示插值的次数，M表示参与插值的总个数，一般在需要插值的区域左右两端各 取2个点参与计算。

图8所示为由求得的视电阻率曲线进行补足前后的对比图，模型参数如下：模 型的电阻率分别为ρ₁＝100Ω·m，ρ₂＝10Ω·m，ρ₂＝100Ω·m，层厚h₁＝20m，h₂＝10m，源1、2 的长度均为1000m，两源之间的夹角为30°，电流大小均为50A，飞行高度100m，偏移距 500m的测点对应的位置坐标为(300m，400m，-100m)。图8(b)中空白的那段曲线即无合适解 和无解情况，与(c)对比可以看出，经过最小曲率补足后，得到的视电阻率曲线能够光滑 渐变完整的反映出所设计模型的电性变化，而且由于该法可以保证曲线按照最小曲率光滑变 化，该法适合对B(t)定义的视电阻率曲线进行补足。

为了说明用进行多辐射源地空系统视电阻率定义的不足，将以数据 为例进行说明。图9所示为分别由B_z(t)和得到的视电阻率对比图，采用的模型同图 8，从图中可以看出，使用B_z(t)定义的全域视电阻率曲线可以光滑渐变完整的反映出模型的 电性信息，而使用得到的视电阻率曲线存在诸如B_y(t)进行视电阻率定义时的无解问 题，但由于利用进行视电阻率定义会带来类似于“overshoot”和“undershoot”的假异常， 最终导致补足等这类问题并不好解决。综合第2节的分析结果，显然对于多辐射源地空系统， 利用B(t)进行数据解释在很多方面都优于如视电阻率曲线更直观、晚期仍能保持较 大的幅值、压制随机高频噪声等。

3模型试验

3.1层状模型

设计了改变第二层电阻率的二层模型对第3节中提及的多分量全域视电阻率定义方 法进行验证，模型参数如下：源1、2的长度均为1000m，两源之间的夹角为30°，电流大小 均为50A，飞行高度100m，偏移距500m的测点对应的位置坐标为(300m，400m，-100m)，二 层模型第一层的电阻率ρ₁＝100Ω·m，层厚h₁＝20m，改变第二层的电阻率 ρ₂＝2，5，10，30，80，200，500，800Ω·m。使用B_p(t)(p＝x，y，z)定义的多分量全域视电阻率曲线见图 (10)，可见随着第二层电阻率的变化，视电阻率曲线呈现出有规律的变化，不仅能够在早 期和晚期逐渐趋于模型的第一层和最后一层电阻率，还能够完整渐变光滑的反映出模型的电 性信息变化。

为了说明偏移距对多辐射源地空系统全域视电阻率的影响，设计了如下模型：源1、2 的长度均为1000m，两源之间的夹角为30°，电流大小均为50A，飞行高度100m，偏移距 500m、1500m、3000m和5000m的测点对应的位置坐标分别为(300m，400m，-100m)， (1060m，1060m，-100m)，(2121m，2121m，-100m)和(3000m，4000m，-100m)，三层模型的电阻率 ρ₁＝500Ω·m，ρ₂＝100Ω·m，ρ₃＝500Ω·m，层厚h₁＝800m，h₂＝800m。从图(11)中可以看出， 全域视电阻率定义方法同时适用于大偏移距和小偏移距，曲线均能完整渐变光滑的反映出地 下的电性信息变化，解决了多辐射源地空系统视电阻率定义距离上不分远近的问题，且通过 比较可以看出偏移距变化对各分量视电阻率曲线的分异均较小。

3.2复杂模型

为了验证多辐射场源地空系统的分辨能力，设计了如下三维模型：在电阻率为100Ω·m的 均匀大地中，赋存一个埋深为100m，电阻率为10Ω·m、尺寸为60m*60m*10m的板状体，见 图(12)，计算采用的源的长度均为110m，电流大小50A，接收高度-80m。三维正演采用矢 量有限元方法完成，首先得到再通过2.3节所示的积分方法得到 B(t)(p＝x，y，z)，利用B(t)(p＝x，y，z)进行数据解释。图13和14的(a)(b)(c)(d)所示分别 为单辐射场源和多辐射场源情况下模型俯视图以及偏移距10m(记为Line24)测线的多分量 全域视电阻率剖面图，图中的虚线表示测线所在位置，粗线表示源所在的位置，图中采用的 色标均相同。通过对比可以看出，由于低阻异常体埋深较大，且规模较小，在单辐射源情况 下几乎分辨不出异常体的存在，但是通过调整源的位置，采用多辐射场源后，可以较为准确 地反映出异常体的存在。

以上所述，仅为本发明较佳的具体实施方式，本发明的保护范围不限于此，任何熟悉本 技术领域的技术人员在本发明披露的技术范围内，可显而易见地得到的技术方案的简单变化 或等效替换均落入本发明的保护范围内。