-宣扬地学成果，传播勘查技术方法-
点击上方“覆盖区找矿”，关注更精彩！
伪随机波形瞬变电磁全时响应数值方法及应用
1 中国矿业大学深地工程智能建造与健康运维全国重点实验室
2 中国矿业大学资源与地球科学学院
作者简介：鲁凯亮，博士，助理研究员，研究时间域电磁正演和数据解释。
通讯作者：岳建华，博士生导师，教授，研究煤田地球物理勘探与煤矿防治水。
传统瞬变电磁法以方波（梯形波）为主，缺乏高频谐波成分，难以精细表征地质异常体。提出在供电期间发射伪随机波形，通过改善场源分辨率，提高瞬变电磁探测能力。给出了伪随机波形全时域电磁正演方法。通过多个阶跃响应的移位线性代数运算，能够计算一维伪随机电磁响应；对于三维数值模拟，基于伪随机波形特征，提出了源项解耦和位移逆Krylov子空间技术，求解全时电磁场仅需一次LU矩阵分解和几十次矩阵回代，并与后退欧拉法的结果对比，验证了本文方法的准确性。通过理论分析、数值模拟和实测数据处理相结合的研究方式，对伪随机波形激励的二次场特征和探测能力进行了详细研究，研究结果表明：伪随机波形相比方波具有更丰富的高频谐波，有效改善了时间域电磁法的探测分辨力，但由于低频谐波能量占比较低，导致晚期电磁场衰减较快；伪随机波形的二次场表达式中包含正负两项，过零点时刻设置不当会导致晚期二次场变号，通过减小心与t+i（i为奇数）之间的时间间隔，可解决这一问题；通过层状模型对比了不同脉宽伪随机波形的二次场，发现窄脉宽波形具有更高的探测分辨力；通过三维地电模型，进一步验证了伪随机波形的探测能力；最后基于实测数据处理，表明伪随机波形激励的二次场能够更加精细地表征地下介质的电阻率特征。基金项目：国家自然科学基金重点资助项目(42230811)；江苏省重大科技示范资助项目(BE2020644)
0 引言
1 伪随机波形时间域电磁数值方法
1.1 伪随机波形一维数值方法
1.2 伪随机波形三维数值方法
1.3 伪随机波形全时段结果验证
2 伪随机波形的频谱特征和电磁响应特征
2.1 伪随机波形的频谱特征
2.2 伪随机波形二次场变号原因分析
2.3 均匀半空间电磁响应特征
2.4 低阻层状模型电磁响应特征
3 三维模型响应特征
3.1 三维低阻矿脉
3.2 三个低阻块状体
4 野外试验
5 结论与展望
时间域电磁法具有多种装置形式，可以满足不同工程应用场景，在地下水探测、有色矿产资源勘探、煤田电磁法勘探以及地质灾害调查等领域得到了广泛应用。但目前时间域电磁法的发射波形主要以方波（梯形波）为主，在野外工作中也不乏三角波、半正弦波以及其他波形的应用，总体来讲，这类波形普遍缺乏高频谐波成分，难以满足不良地质体精细表征的重大需求。本文从发射波形入手，通过补充高频谐波成分，改善时间域电磁法的探测分辨力。为了提高探测分辨力，众多学者将On-time的发射波形作为研究切入点。Liu于1998年研究了不同发射波形对机载电磁系统响应的影响，给出了常用的半正弦波、方波、三角波和梯形波的最佳脉宽和探测时间之间的关系。陈曙东等针对这四种发射波形，推导了回线源在自由空间中的电磁场表达式，研究了对电磁探测的影响，并指出方波与梯形波适用于地面电磁法，半正弦波和三角波适用于航空电磁法。殷长春等详细研究了半正弦波和梯形波对地质目标体的探测能力，根据实际勘探要求，可以通过设计最佳的航空电磁参数，进而提高工作效率。根据赵越等人的研究，这类发射波形缺乏高频谐波成分，导致时间域电磁探测分辨率难以得到根本改善。李文翰等在城市地下空间成像方法中指出，微分脉冲（由两个脉宽相同的连续正负方波组成）可以减少低频谐波干扰，提高探测分辨率。本文根据这一思路，将频率域电磁法中广泛应用的伪随机波形，扩展至时间域电磁法。针对伪随机复杂波形的数值模拟方法，可将其归类为考虑发射波形的时间域电磁正演。对于一维正演，殷长春等利用系统阶跃响应和电流波形导数进行褶积，计算了全时段航空电磁响应。岳鑫使用阶跃响应移位相加减，计算了微分脉冲的二次场，本文将其扩展至伪随机波形一维正演。对于三维正演，齐彦福等对全时段波形采用后退欧拉法，使用多个时间步长的离散策略，实现了航空电磁全波形三维正演。但是这种计算方法涉及多次矩阵分解，计算成本较高。Zhou等将磁场控制方程写成偏微分方程形式，使用位移逆Krylov子空间方法直接求解全时段电磁响应。但该方法在计算On-time电磁响应时，每个时刻均需重新构建Krylov子空间，耗时较长。本文基于源项解耦和位移逆Krylov子空间技术，提出基于四面体网格的三维快速正演方法，计算全时段电磁场只需要一次矩阵分解和几十次矩阵回代，特别适合伪随机一类复杂波形的正演模拟。本文内容主要分为以下三个部分：首先给出了伪随机波形时间域电磁一维和三维正演方法，研究了伪随机波形的频谱特征与电磁响应特征；然后通过层状模型和三维地电模型，验证了伪随机波形具有较高的探测分辨力；通过野外试验，表明伪随机波形具有较大的发展潜力。最后对本文内容进行总结，并对存在的问题以及下一步研究工作进行说明。伪随机波形如图1a所示，可看作由多个连续的正负方波（图1b）组成，本节内容首先给出对应的数值模拟方法。
图1 波形示意图（a）伪随机波形；（b）方波
对于伪随机波形的一维数值模拟，本文给出两种求解策略。第一种方法使用阶跃响应与发射波形导数进行褶积，本文做简要提及：上式中，B(t)是伪随机波形的磁感应强度，B(t)和Bs(t)分别是脉冲函数与下阶跃函数的磁感应强度，I(t)是伪随机发射电流。其中阶跃波的电磁响应可由Key的开源代码计算。对于第二种方法，本文提出伪随机波形可以看作由多个不同脉宽的正负方波组成，每个方波可以等效为两个阶跃波相减。这样即可通过多个阶跃响应的移位线性代数运算得到伪随机波形的电磁场响应。对于由n（n大于1且为正整数）个不同正负方波组成的伪随机波形，将每个方波的过零点时刻分别记作t0、t1和tn（图1a）。以磁感应强度B(t)为例，第i个方波（图1b）的磁感应强度为：对于伪随机波形三维数值模拟，本文采用非结构四面体矢量有限元方法。通过观察伪随机波形可知，每个方波的持续时间不同，若整体采用较小的时间步长离散，势必造成时间步过多；若采用不同时间步长离散，则需要多次系数矩阵分解(因为系数矩阵是时间步长的矩阵函数)。为了解决这一难题，本文基于源项解耦和位移逆Krylov子空间技术，将伪随机波形解耦为时变电流函数和空间分布常数的乘积，只需要1次矩阵分解和数十次矩阵回代，特别适合研究伪随机一类复杂波形的电磁场特征。由Maxwell方程组可知，电场的双旋度方程如下：其中E为电场强度，m和s分别是介质的磁导率和电导率，t是时间，Js是外加电流密度。对方程(8)使用伽辽金法，可得：由于源项可以表示为固定空间分布的矢量场和随时间变化函数的解耦形式，因此向量g(t)可以写成下面的形式，其中u为常向量。
通过Krylov子空间基向量矩阵Vm和投影矩阵Hm满足的关系Hm=VTmAVm，本文对方程（11）两边同时乘以矩阵VTm，并由VmVTm=I可得： 求得电场强度后，可通过Maxwell方程组求取衰减电压和磁感应强度。
本文采用脉宽为1ms、2ms、3ms和4ms的正负方波组成伪随机波形，设置接地导线源长度100m，发射源端点的坐标分别是(-50m，50m，0m)和(50m，50m，0m)，接收点坐标为(0m，0m，0m)。空气层和均匀半空间地层的电阻率分别设为106Ω.m和100Ω.m。图2是均匀半空间计算结果，其中图2a和图2d是伪随机发射波形；图2b和图2e分别是衰减电压z分量和磁感应强度z分量，实线、圆圈和星分别表示一维、本文方法和后退欧拉法的计算结果；图2c和图2f分别是衰减电压z分量和磁感应强度z分量的相对误差曲线，从图中可看出，除了发射电流突变点外，本文方法和后退欧拉法的相对误差整体较小，证明了其准确性。此外相比于传统的后退欧拉法，本文方法的计算效率可提高10倍以上，如表1所示。 图2 均匀半空间结果
（a），（d）波形；（b），（e）衰减电压与磁感应强度；（c），（f）相对误差
表1 两种方法对伪随机波形的运行时间
目前时间域电磁法的发射波形主要以方波(梯形波)为主，缺乏高频谐波成分，难以满足日益增长的精细勘探需求，本节内容研究伪随机波形的频谱特征与电磁响应特征。伪随机波形较为复杂，研究其频谱特性以及波形控制对时间域电磁探测具有重要意义。本节以四个不同脉宽的正负方波组成伪随机波形作为研究对象。为了便于对比，本文同时使用脉宽10ms的方波作为对照组，如图3所示。设置方波占空比为1:1，从图3b可以看出，大约在100Hz时方波的频谱振幅衰减为零，因此本文设置100Hz为截至频率，小于100Hz的谐波视为低频成分，大于100Hz的谐波视为高频成分。本文利用高频谐波与横坐标所围面积占总面积的百分比，可计算方波频谱的高频谐波能量占比约为68.9%。首先从简单波形入手，研究两个正负5ms方波组成的伪随机波形，如图4a所示，相对于方波频谱（图3b），图4a对应频谱的高频谐波能量占比约为78.9%（图4b）。继续减小每个方波脉宽，并保持每个正负方波的脉宽一致。当方波脉宽减小至2.5ms时（图4c），低频谐波能量占比减少，高频谐波能量占比进一步增强，约为94.0%（图4d）。随着每个方波脉宽进一步减小至1.25ms和0.625ms（图4e和图4g），对应的高频谐波能量所占比重分别是97.4%和98.7%。图4 相同脉宽方波组成的波形及频谱图
（a）、（c）、（e）和（g）中方波数目依次是2、4、8、16；（b）、（d）、（f）和（h）为不同波形对应的频谱
通过研究大量的伪随机波形（图5a、图5c、图5e和图5g），发现相比于传统发射波形，伪随机波形可以有效增加高频段的谐波能量（图5b、图5d、图5f和图5h）。同时由于低频谐波能量占比较低，导致晚期电磁信号衰减过快，信噪比较低，加大了对深部空间的探测难度。因此对于伪随机波形需要配备大功率发射机，同时减小电性源的接地电阻，提高辐射电流强度。图5 伪随机波形及频谱图
（a）、（c）、（e）和（g）不同伪随机波形；（b）、（d）、（f）和（h）对应频谱图
目前瞬变电磁数据处理方法主要针对关断电流后的电磁数据，因此研究关断后的电磁场特征尤为重要。在大量的数值模拟中，本文发现晚期结果时常出现变号现象。通过观察伪随机波形电磁场Bn(t)的表达式可知：公式中除了第一项B(t)，其余项均含系数(-1)i，因此在某些情况下电磁场会发生变号。将公式（7）展开发现，其中B(t-ti)-B(t-ti+1)等项恒大于零（i为正偶数），-B(t-ti)+B(t-ti+1)等项恒小于零（i为正奇数），为了保证Bn(t)正常衰减，则需要-B(t-ti)+B(t-ti+1)等项尽可能小，即使时间ti与ti+1之间的时间间隔尽量小(i为正奇数)。本文设置如下四组伪随机波形进行论证，令伪随机波形的过零点时刻为t1至t4。波形1中t1与t2持续时间为9ms，波形2中t1与t2持续时间为1ms；波形3中t1与t2、t3与t4持续时间为4ms，波形4中t1与t2、t3与t4持续时间为1ms，如表2所示。
表2 四种波形的过零点时刻
图6是四种不同波形的计算结果，图中实线表示正值，虚线表示负值。图6a至图6d对应的发射波形为波形1至波形4。从图6可以看出，波形1中t1与t2的持续时间和波形3中tl与t2、t3与t4的持续时间均较长，导致晚期电磁信号出现负值。而波形2与波形4的过零点时刻满足前文分析，所以在整个时间段上二次场没有变号现象(图6b和图6d)。图6 四种不同波形计算结果
（a）波形1；（b）波形2；（c）波形3；（d）波形4
本节研究伪随机波形均匀半空间的电磁响应特征，收发装置和物性参数与前文一致。设置方波和三种伪随机波形，其中方波的持续时间为10ms，波形1与波形2每个方波持续时间为5ms和2.5ms，波形3每个方波持续时间依次为3ms、2ms、1ms和4ms，如表3所示。表3 四种波形的过零点时刻
图7为表3中四种波形的计算结果。其中图7a和图7b分别是衰减电压z分量和磁感应强度z分量。从图中可以看出，当发射波形的供电时间固定时(此处供电时间均为10ms)，离散的方波数目越多，晚期信号衰减越快，这是因为低频谐波的能量占比随离散方波数目的增多而降低(图8a至图8c)，这与前文分析较为一致。由于波形2和波形3的谐波特征接近，所以波形2和波形3的衰减曲线没有较大差异。图7 四种不同波形计算结果
（a）衰减电压z分量；（b）磁感应强度z分量
图8 四种不同波形频谱
（a）方波；（b）波形1；（c）波形2；（d）波形3
目前针对特定地电模型的最优伪随机波形设计较为困难，因此在接下来的内容中，本文以两个均匀分布的正负方波组成伪随机波形进行研究。当均匀地层中含有异常体时，感应电动势曲线的“Undershoot”和“Overshoot”现象会造成假异常，因此本文采用磁感应强度z分量反映伪随机波形二次场的相对异常。本文设置图9所示的层状模型，均匀半空间的电阻率为100 Ω.m，其中含有两个低阻薄层。其中第一个低阻薄层的厚度为10m，埋深30m，电阻率为5Ω.m。第二个低阻薄层的厚度为10m，埋深130m，电阻率为1 Ω.m。电性源沿y轴布设，长度200m，发射源的两个端点坐标分别为(-100m，100m，0m)和(100m，100m，0m)。接收点坐标为(0m，0m，0m)。图9 层状模型
本文设置伪随机波形和方波的供电时间分别为0.5ms、1.5ms、3ms和10ms，并对比两种波形对低阻地层的分辨能力。本文定义相对异常如下：其中fhalfspace表示同装置下，均匀半空间磁感应强度z分量，f3d表示含有异常体(层)的磁感应强度z分量。
图10是不同波形的计算结果，其中图10a和图10b分别是激励方波时，均匀半空间(100Ω.m)和层状模型(图9)的磁感应强度z分量曲线，图10c和图10d分别是激励伪随机波形时，均匀半空间和层状模型的磁感应强度z分量曲线，图10e和图10f是不同方波和伪随机波形的相对异常曲线。从图10e可以看出，当方波的脉宽为0.5ms时，两个低阻薄层的相对异常较为明显，当方波的脉宽进一步增大时，浅部低阻薄层引起的电磁异常迅速减弱。从图10f可以看出，当伪随机波形的脉宽为0.5ms时，两个低阻薄层的电磁异常均比图10e明显。与方波的电磁异常规律类似，随着伪随机波形脉宽的增大，高频谐波能量占比迅速减少，从而降低了浅部薄层的分辨能力。尽管小脉宽的波形具有更强分辨力，但是其低频谐波能量占比更低，晚期电磁场衰减更快，因此在探测分辨力和晚期信噪比之间要取得平衡。
图10 不同波形结果
(a)方波均匀半空间响应；(b)方波层状响应；(c)伪随机均匀半空间响应；(d)伪随机层状响应；(e)方波相对异常；(f)伪随机相对异常
本文设置图11所示的三维低阻矿脉模型。模型参数设置如下，空气电阻率为106Ω.m，均匀半空间电阻率为200 Ω.m。矿脉顶部埋深60m，沿x方向的长度为600m，沿y方向延伸400m，沿z方向厚度为200m，矿脉的电阻率为2 Ω.m。剖分最小网格长度为8m，网格扩大因子为1.3。接地导线源长度1km，发射源两个端点坐标分别为(-500m，100m，0m)和(500m，100m，0m)。测线长度800m，偏移距为100m，测线两端坐标分别是(-400m，0m，0皿)和(400m，0m，0m)，点距20m，共41个测点，方波与伪随机波形的供电时间均为1ms。图11 三维低阻矿脉模型
图12是方波和伪随机波形在xoz平面的磁感应强度z分量相对异常等值线图，对应的关断后时间依次为0.05ms、0.12ms和1ms，图12a至图12c是方波的相对异常，图12d至图12f是伪随机波形的相对异常。从图上整体可以看出，磁感应强度z分量的相对异常呈现先增大后减小的特征。通过对比可以看出，图12a和图12d的相对异常都只能反映低阻异常体的浅部形态，且两者差异较小。随着时间推移，当关断后0.12ms时，图12b和图12e可以更进一步表征异常体的形态。并且进一步发现，图12e的相对异常明显较强。当关断后1ms时，电磁场进一步向下传播，通过图12c和图12f可以发现，此时的异常等值线已经无法清晰刻画异常体的浅部特征，只能较为清楚地反映异常体的深部特征。相应地，图12f中的相对异常也明显强于图12c。图12 不同时刻磁感应强度z分量相对异常等值线图
（a）-（c）方波；（d）-（f）伪随机波形
图13是方波和伪随机波形地表测线数据的相对异常，从图中可以看出，两种波形的相对异常形态大致相同，但是伪随机波形的相对异常更加明显。通过该模型，可以证明伪随机波形在时间域电磁勘探中具有较强分辨率。图13 地表测线磁感应强度z分量相对异常
（a）方波；（b）伪随机波形
本文设置图14所示的三个低阻块状体模型。模型参数设置如下，空气电阻率和半空间电阻率分别是106 Ω.m和100Ω.m。最小网格长度为10m，网格扩大因子为1.3。接地导线源置于地表，长度1km，发射源两个端点坐标分别为(-500m，100m，0m)和(500m，100m，0m)。三个低阻异常体的顶部埋深分别为30m、100m和180m，异常体的电阻率分别为3Ω.m、2Ω.m和1Ω.m；几何尺寸分别为40m×40m×10m、80mx80mx10m和120mx120mx10m。在地表设置长度为400m的测线，起始点与终止点的坐标分别为(-200m，0m，0m)和(200m，0m，0m)，点距10m，共41个测点。方波与伪随机波形的供电时间均为0.05ms，研究两种波形对该模型的分辨能力。图15是不同时刻xoz平面磁感应强度z分量相对异常等值线图，其中图15a至图15c是方波的相对异常，图15d至图15f是伪随机波形的相对异常。观测时刻分别是关断电流后0.1ms、1ms和3ms。图15a和图15d是关断电流后0.1ms的相对异常，从图上可以看出，两种波形对三个块状体都有反映，但是伪随机波形的异常信息更加丰富。图15b和图15e是关断电流后1ms的相对异常，此时方波与伪随机波形的相对异常明显增强，伪随机波形对第二个和第三个低阻异常体的反映更为明显。图15c和图15f是关断电流后3ms的相对异常，随着电磁场向下扩散，方波与伪随机波形几乎都无法分辨第一个低阻异常体，伪随机波形对第三个异常体有较为明显的反映，而方波的异常信息相对微弱。图15 不同时刻磁感应强度z分量相对异常等值线图
(a)-(c)方波；(d)-(f)伪随机波形
图16是该模型地表测线数据的相对异常，可以看出，方波的相对异常较弱，几乎无法分辨第二个异常体；而伪随机波形的相对异常不仅更强，还对第二个异常体有更好的分辨能力。图16 地表测线磁感应强度z分量相对异常
（a）方波；（b）伪随机波形
本次试验位于山东省某矿区，采用无人机地空电磁法，以接地长导线作为发射源，在一次场间歇期，以无人机载平台搭载传感器在空中接收电磁场。试验采用的装备主要包括发电系统、信号发射系统和无人机采集信号系统。将发射系统与接收系统的时间同步，编写程序(已知二进制存储文件的数据格式)提取伪随机波形激励的二次场。图17 野外工区示意图
试验矿区地势平坦，区内第四纪地层分布广泛，深部以焦家断裂带主裂面为界，上部主要为马连庄序列栾家寨岩体，下部为各含矿物蚀变岩带。工区地面标高位于11.2m至18.6m，图17为野外工区示意图。根据试验需求及现场施工条件，设置地面发射源长度为1400m，发射电流25A，整个试验期间发射电流状态稳定。布置6条平行于发射源的无人机飞行测线，飞行高度约为10m，最小偏移距为350m，每条测线长度500m，测线间距50m，距发射源由近及远依次是测线L1至测线L6。图18为野外铺设电极板与电磁场采集照片，图19为发射波形示意图，方波与伪随机波形的供电时间均为10ms。 图18 铺设电极板与电磁场采集
图19 发射波形
图20是方波和伪随机波形在测线L3中心处测点的衰减曲线。早期阶段，两种波形的衰减曲线几乎重合；1ms前后，由于伪随机波形高频谐波成分更丰富，使得两条曲线有较为明显的差异；晚期阶段伪随机波形低频谐波能量占比较低，导致电磁响应衰减较快。图21是使用全域视电阻率方法计算的结果，图21a和图21b分别是方波和伪随机波形的视电阻率图。根据钻孔地质资料：①地下250.1m至420.3m，水系较为发育，表现为相对低阻特征；②地下663.9m至1063.0m，岩石裂隙发育，裂隙面见有绢云母、绿泥石等蚀变矿物、灰白色构造角砾岩等，导电性较差，电阻率相对较高；③地下1089.0m至1404.5m，主体为花岗碎裂岩，岩石裂隙发育且呈现弱富水性，电阻率相对较低。视电阻率整体随深度的变化特征与已知地质资料吻合较好。通过对比可知，方波的视电阻率整体呈“团状"分布，伪随机波形的结果可以更为准确地表征地下250.1m至420.3m的低阻层位置信息（图21中黑色虚线），能够突出地下电性介质的局部特征。图21 测线L3视电阻率
本文研究了伪随机发射波形在时间域电磁探测中的数值方法与应用，提出了适用于伪随机波形的1D和3D时间域电磁数值模拟方法，并给出1D数值解的表达式。通过研究有以下几点结论：（1）相较于传统波形，伪随机波形具有更丰富的高频谐波，可以提高浅层地质目标体的探测分辨力；（2）伪随机波形的低频谐波能量占比较低，因此在野外作中必须配备大功率发射机，减小接地电阻，提高发射电流强度；(3)数值结果表明，伪随机波形脉宽越窄，高频谐波越丰富，电磁异常越明显，同时导致晚期电磁场衰减更快，因此在探测分辨力与晚期信噪比之间要取得平衡；(4)伪随机波形的过零点时刻设置不当会使晚期电磁场发生变号，本文研究了背后机理并针对这一现象提出了解决方法；(5)通过数值模拟和野外试验，证明了伪随机波形在时间域电磁探测中具有较大发展潜力。
展望如下：(1)伪随机波形对深部地层的探测分辨力有待进一步研究；(2)研究叠加次数对实际工作中伪随机电磁场的影响；(3)对于野外实际工作，如何设计针对特定地质构造的最优伪随机波形需要重点研究；(4)多辐射源伪随机时间域电磁探测技术和数据解释方法有待深入研究；(5)伪随机电磁场的三维反演以及伪地震子波成像也是后续研究的重点方向。原文来源：鲁凯亮，岳建华，樊亚楠，习丹阳．伪随机波形瞬变电磁全时响应数值方法及应用[J/OL].煤炭学报.
https://doi.Org/10.13225/j.cnki.jccs.2023.1755
导读评论和排版整理等：《覆盖区找矿》公众号.
推荐读者下载、阅读和引用原文！
------------往期精彩回顾-------------
------关注“覆盖区找矿”，拥有更多新方法------![]()
