1 引言

近年来，传统的电磁方法在能源、矿产资源、水资源、环境地质及工程地质勘察中都得到了广泛应用且发挥了重要作用(底青云等，2006；Xue et al., 2014).随着国民经济的发展，需要探测埋深更深、精度更高的目标体时，这些方法的应用受到了一定的限制，亟待开发新的技术来提高探测精度与探测深度.

近年来，一些新技术新方法被陆续提出来，英国爱丁堡大学的Wright (2002)提出了多道瞬变电磁法的概念和探测油气目标体的处理技术，何继善(2010)提出了可提高分辨率及加大勘探深度的广域电磁法及相关技术.

在这些新方法中，都引入了伪随机编码发射技术，其全新的处理方式及应用效果引起各个研究机构的关注.近几年，中国科学院、中南大学、中国地质大学、国土资源部以及相关单位对这种发射技术(赵璧如等，2006；罗先中等，2014；罗延钟等，2015)及后续的处理工作都进行了不同程度的研究.汤井田和罗维斌(2008)分析了逆重复m序列的波形特征、相关函数及频谱特性，讨论了对各种干扰的压制方法及测量精度，薛国强等(2015)分析了多通道瞬变电磁法(MTEM)伪随机编码发射的响应特性和方法的关键技术，齐彦福等(2015)利用方波响应移位叠加和电流导数与阶跃响应褶积两种方法实现了理论m序列和实际发射波形的全时正演模拟，武欣等(2015)研究了m序列伪随机编码源电磁响应的精细辨识问题.

学者多把研究重点放在相关辨识及对噪声的压制上，对相关辨识的方法理论做了充分的探讨，但在模拟电磁信号方面却着墨较少.然而，采用伪随机编码源的资料处理工作大多是从电磁信号出发的，模拟伪随机编码源激发下的电磁信号不但可以为资料处理环节提供理论数据，而且可以为检测资料处理的质量提供中间结果，因此，模拟电磁信号工作是必不可少的.

本文根据获得接收信号的物理过程来实现伪随机编码源激发下的电磁信号合成.首先简述伪随机源激发下的电磁信号合成基本理论，然后根据合成公式中的大地脉冲响应、伪随机编码源、仪器带限影响及噪声逐一进行研究，最后给出信号合成结果，并与一个野外实测信号作对比，来验证信号合成的正确性.

2 基本理论

伪随机源激发下的电磁信号合成工作是指给定一个m序列伪随机编码源，在编码源的激发下，模拟大地及发射、接收系统和噪声的共同响应，继而合成一个仿真的野外观测信号的过程.

在野外工作中，发射机通过两个接地电极向地下供入伪随机交变电流信号，在离发射电极一定距离处用不极化电极排列接收相邻电极间的电位差，根据这种发射接收装置，发射机发射的电磁信号在大地中传播，然后通过电极与导线、接收系统接收.在勘探工作中，大地是探测目标，经过大地的信号为有用信号.对有用信号产生影响的因素主要为噪声，除此之外，电极附近的极化效应、接收器的频率响应限制(即带限)、电线的感抗等都对有用信号或多或少产生影响.在这些影响因素中，除噪声外，每一部分的影响都相当于在有用信号上加载了一个滤波器，最后接收到的信号是对有用信号经过多次滤波后的综合信号.因此，模拟电磁信号工作应当考虑各种因素的影响，然而，现实中有些因素的影响不易模拟或相对较小，本文只考虑噪声及仪器带限影响，并且将发射和接收系统的带限影响综合考虑.

接收信号的合成公式可以用(1)式表示(Ziolkowski et al., 2007)：

(1)

其中，t表示接收时间，v(t)表示接收机在接收点接收的信号，i(t)表示发射系统向大地发射的电流，g(t)表示经由大地滤波后在接收点处的大地脉冲响应，r(t)表示仪器频率响应.n(t)为接收点处的噪声.

从公式(1)可以看出，只要获得了模型地电结构的大地脉冲响应g(t)、伪随机发射源信号i(t)、仪器带限响应r(t)及噪声n(t)，就可合成出仿真的电磁信号v(t).下面对公式(1)中各项的实现过程逐一阐述.

3 大地脉冲响应

理想情况下，向地下输入脉冲电源，在不计及各种影响因素时，接收机测量的输出即大地脉冲响应.但在野外工作中，一般采用方波作为输入.在数值模拟时，先在频率域获得不同频率的电磁响应，然后通过余弦变换转换到时间域，得到大地阶跃响应，再通过计算阶跃响应的时间导数来得到大地脉冲响应.

3.1 频率域响应

本文用均匀半空间地电模型来模拟地表电场信号的合成过程.根据文献(Wright，2004)，只有与源同方向的电场可提供有用信息，其他方向的场不能提供更多的信息，因此本文只用与源同方向的电场来合成所需信号.假设源的方向为x方向，则频率域中同方向的电场可写为(朴化荣，1990)：

(2)

其中，P_E为电偶极距，ρ₁为均匀半空间的视电阻率，r为收发距，φ为接收点与发射中点的连线和发射偶极的夹角，k₁为波数，k₁²=-iωμ₀σ₁，σ₁为均匀半空间的电导率，μ₀为导磁率，ω为角频率，i表示纯虚数.

3.2 频时变换

用于频率域到时间域变换的方法较多，如正弦变换或余弦变换(Anderson, 1982；Newman et al., 1986; Li and Constable, 2010)、折线法(Li et al., 2011)，这三种方法同属于傅里叶逆变换.在这三种方法中，正弦变换与余弦变换需要计算较宽频带内的频率响应，其频带宽度与时间点有关；折线法相对来说计算速度快，但相对粗糙.除此之外，还有基于拉普拉斯逆变换的系列算法(Li and Farquharson, 2015)可以实现信号的频时域转换，这类算法的优势在于可以用较少的滤波系数得到较精确的结果，但其对计算机精度的要求太高(4倍精度).本文最终选择了余弦变换作为把信号从频率域转换到时间域的方法.

余弦变换的表达式为(朴化荣，1990)：

(3)

其中，f(t)为时间域的场，F(ω)为频率域的场.

参考Anderson (1982)的方法，选用函数对

(4)

来求余弦变换的滤波权系数，其中，a为大于0的实数.利用滤波系数可将(3)式写为：

(5)

其中C_i为余弦变换的滤波权系数，A_i-x为移动的横坐标，N₁和N₂的大小由积分时间自动调节，本文所用的余弦变换子程序为公开的代码(方文藻等，1993).

3.3 大地脉冲响应

殷长春等(2013)利用阶跃响应与电流的时间导数褶积代替脉冲响应与电流褶积的方法成功实现了任意波形全时响应的计算，避免了记录时间t→0时大地脉冲奇异值的出现，但只考虑了电流变化感应出的电磁场，并没有考虑电流自身产生的直流场，所以还需再将直流场加进感应的场中.本文采用Strack (1992)的提法，直接计算阶跃响应的时间导数来得到大地脉冲响应.t→0处的值用第一个阶跃值与该处时间的比值来代替，通过将起始时间设计得较小(如0.01 ms)来减小t→0时的误差.

设计一个电阻率是100 Ωm的均匀半空间模型，假设源的方向为x方向，源的长度为200 m，计算点位于源的延长线上且离源中心点500 m处，电场的频率域响应曲线如图 1a所示.图 1a显示了从10^-7 Hz到10⁷ Hz的计算结果，这个频带宽度可满足1 s内所有时间点上进行余弦变换时所需的频率范围.从图 1a可以看出，在高频和在低频一定范围内，电场曲线保持水平，不受频率的影响.通过余弦变换将频率域结果转换到时间域，得到如图 1b所示的阶跃响应，将其与解析公式计算结果对比，发现二者的吻合程度非常高，验证了余弦变换过程的正确性.采用3.3节中提到的阶跃响应的时间导数来计算脉冲响应，计算结果如图 1c所示，对图中所示的x方向电场来说，曲线中出现了一个峰值，这个峰值对后续的处理非常重要，因不是本文的研究内容，所以不进行阐述.这里为了突出峰值的形态，对曲线进行了截断.所用的时间域解析公式来自文献(朴化荣，1990).

4 m序列伪随机源 4.1 m序列伪随机编码

伪随机序列是由移位寄存器产生的，根据文献(林可祥和汪一飞，1978)，可通过查表的方式得到不同级移位寄存器产生的伪随机序列的本原多项式F(x).通过F(x)与序列多项式G(x)的关系式G(x)=1/F(x)，用多项式的长除算法可计算出相应的序列多项式G(x)，从G(x)中可直接获得m序列伪随机编码.

序列多项式G(x)最终可写为下面的表达形式：

(6)

其中⊕表示模二加法运算符，x为移位运算符，其上标的数字表示移位的位数，例如：a_kx²表示将a_k位移两位，a_n所组成的序列{a_n}=a₀, a₁, a₂, a₃, …即为m序列伪随机编码.(5)式中的N为伪随机码一个周期的码元个数，N=2^r-1，r为移位寄存器级数.

图 2是4级和6级的m序列伪随机编码一个周期的波形，它是通过文献(林可祥和汪一飞，1978)本原多项式列表中4级和6级第一个本原多项式得到的m序列的编码波形，从图中看出其随机性随着移位寄存器级数的增多而愈发明显.

4.2 m序列伪随机编码源的频谱

对m序列伪随机编码应用快速傅氏变换，获得了其作为源时相应的频谱.通过对比不同级数的频谱幅值，可了解不同级数的伪随机信号所包含的频率成分及幅频特性.本文给出了4、6、8、10、12级伪随机编码源的频谱，其中与伪随机编码相关的参数见表 1.为了对比不同源谱的宽度，将所有伪随机编码周期与采样率都定为0.25 s及16384 Hz.伪随机编码源的频谱如图 3所示.

从图 3可以看出，不同级移位寄存器得到的频谱形态不同，级数越高，主瓣的宽度越宽，当移位寄存器的级数r=12时，从0 Hz到8192 Hz范围内的频谱丰富，包络线的形态接近水平，接近于脉冲源的频谱，由此源激发得到的大地响应会接近于脉冲源激发得到的响应.图 3中从10级到12级伪随机编码源的频谱外包络线不光滑，这是采样率不够导致的.若提高采样率，则能得到外包络线较光滑、幅值较均匀的谱图.从图 3还可以看出，随着寄存器级数的增加，主瓣的幅值也随之下降，说明信号的能量降低.对地球物理勘探来说，信号能量降低会导致电磁信号的穿透深度变浅.换句话说，虽然主瓣频率成分提高，分辨率也提高，但信号强度因能量分散而降低，导致探测深度浅.若想解决这个问题，需在加大发射功率的同时增加垂直叠加的时间，而且在后续处理采用类地震水平叠加的技术.

5 其他影响因素 5.1 仪器频带限制的影响

发射机和接收器自身的频带宽度是有限的，作为理论研究，可以用有一定带宽的低通滤波器来模拟仪器的频带宽度(即带限问题).

低通滤波器的设计方法有多种，不同方法得到的滤波器的旁瓣幅值大小以及从主瓣到旁瓣的过渡带宽度不同.理想滤波器要求过渡带尽可能窄，旁瓣幅值尽可能小，以减小能量的泄露，使能量尽可能多地保留在主瓣内.本文选用海明窗函数法作为滤波器的设计方法，因为海明窗可将99.963%的能量集中在窗谱的主瓣内，旁瓣的峰值小于主瓣峰值的1%(程佩青，2007)，符合理想滤波器的要求.

海明窗函数为：

(7)

其中，n为滤波点数，ω_c为截止频率，α为滤波器的中心，，N为物理滤波器的阶数.

借鉴频率域电磁仪器的频带宽度，设计一个低通滤波器，其通带为8192 Hz，阻带为9600 Hz.理想的低通滤波器的频带如图 4a所示，在通带内，信号可以通过，通带外是阻带，阻带内的信号全部被滤除.理想的低通滤波器在时间域内如图 4c中的黑实线所示(实际为无限信号，为了和加窗后的信号相比，只显示了其中的一部分).图 4b是海明窗在频率域的表现形式，用海明窗对图 4c中的信号进行截断，截断后的信号如图 4c中的离散点所示，将加窗后的时间域滤波器变换到频率域，如图 4d所示，相对于理想低通滤波器的频带(图 4a)，图 4d出现了通带与阻带之间的过渡带，但较窄，阻带出现了小幅振荡(因其相对于主瓣来说，能量较小，振荡不明显，若放大该部分，可以看到振荡现象)，说明所设计的滤波器过渡带较窄，阻带能量泄漏很小，是一个性能优良的滤波器.

5.2 噪声的加入

当用电磁法在野外工作时，常遇到的噪声类型有两种，一是白噪声，二是50 Hz的工频噪声.

(1) 白噪声

本文用正态分布随机数来模拟噪声：

(8)

其中，n₁为所求的正态分布随机数，μ为正态分布的均值，本文为观测资料，σ²为正态分布的方差，在本文中σ为观测值与噪声水平之积，N_i为0到1之间均匀分布的随机数，m为随机数的个数.

(2) 50 Hz噪声

用一个频率是50 Hz的正弦波来模拟工频噪声：

(9)

其中n₂为噪声值，A为50 Hz工频噪声的振幅.

6 伪随机源激发下的电磁信号合成 6.1 信号合成

通过(1)式将以上内容进行组合，便可分别得到理想接收信号、加带限的信号和带噪声的信号.

以4阶m序列伪随机编码源作为激发源(见图 5a)，将源和均匀半空间大地脉冲响应(见图 1c)做褶积，得到电场的理想信号，如图 5b所示.将理想信号与低通滤波器(仪器带限)做褶积，便可得到考虑了带限影响的电场信号，如图 5c所示.从这些图可以看出，不加带限前，信号比较平滑，加了带限后，信号出现了一些毛刺状的干扰.

对加了带限的信号分别加上了5%及10%的白噪声，如图 6a和6b所示，噪声叠加在原来的信号之上，引起信号的畸变，当噪声达到10%时，形成的干扰已非常明显.图 6c为幅值是4倍于有用信号的50 Hz干扰信号，干扰信号和有用信号叠加后的信号如图 6d所示，可以看出，有用信号湮没在了50 Hz干扰信号中.

6.2 合成信号与实测信号的对比

中国科学院项目组2015年在河北张北地区进行了野外试验.用一个长度是300 m的接地电偶极发射一个12阶的伪随机编码源，采用轴向装置，用一系列极距是60 m的电偶极子一字排开，在离发射源不同距离处的接收端接收耦合了大地响应及各种影响的信号.为了和野外信号对比，在室内也进行了模拟.根据测区已知资料，在室内用电阻率是100 Ωm的均匀半空间模型来近似地下介质，发射源也用12级移位寄存器生成的m序列伪随机编码序列，编码的部分形态如图 7a所示.采用了和野外一致的码元长度及采样率(码元长度为1/1024 s，采样率为16 kHz，发射一个周期的信号大约需要4 s)，因信号太长，为了突出细节，我们仅显示了其中的一小部分.值得注意的是，室内的编码序列和室外的编码序列是不一致的，因为对于任何级的移位寄存器来说，可以产生很多种形态编码，级数越多，能产生的编码形态也越多，所以并未刻意和野外用相同的编码形式.模拟时，分别在偏移距由近及远的位置接收信号，300 m和2700 m处的脉冲响应分别如图 7b和图 7c所示，离源较近的接收点脉冲响应的信号很强，并且和离源较远时接收到的脉冲响应形态有差异，这是因为离源近，信号接近饱和，离源越远，得到的脉冲响应越平缓.

由于本文所依据的理论是电偶极子理论，当发射偶极与接收偶极子距离较近时，对发射偶极源进行了离散叠加.对300 m的收发距，将偶极源划分成一系列5 m长的小偶极子，离散后的小偶极子是300 m收发距的1/60，满足偶极子理论，同时大大减小了线间的感应耦合.

我们计算了部分点的接收信号，并和野外接收信号做了对比，发现二者的相似程度较高，这里只展示偏移距是300 m时的合成信号，如图 7c所示，通过和野外300 m处接收信号的屏幕截图 7e对比，可以看出信号的相似度较高.二者不尽相同的原因有三个：一是所用的编码形式不同，二是地下介质的电阻率不同，这两个在前面也有所描述，还有一个原因，就是在野外实测中所受的干扰更多，如电极附近的极化效应，接收电线的感抗，而在数值模拟中只是模拟了较为典型的干扰.即便如此，二者较高的相似度说明了本文信号合成过程是正确的.

7 结论

本文根据伪随机编码源信号的传播物理过程对接收信号进行了合成，通过研究，发现大地脉冲响应离源较近时会出现饱和现象，是一条单调衰减曲线，随着测点远离源，峰值逐渐出现.

本文通过计算阶跃响应的时间导数来得到大地脉冲响应，t→0处的值用第一个阶跃值与该处时间的比值来代替，通过将起始时间设计得较小(如0.01 ms)来减小t→0时的误差，通过数值模拟发现，这种近似方法可行.

当用移位寄存器得到m序列伪随机编码时，移位寄存器级数越多，其频谱主瓣的宽度越宽，当移位寄存器的级数r=12时，从0~8192 Hz范围内的各个频率点上的幅值外包络线形态接近水平，接近于脉冲响应的频谱，由此源激发而得到大地响应会接近于脉冲源激发得到的响应，从而有较高的电性结构分辨能力.然而随着寄存器级数的增加，主瓣的各个频率信号的幅值也随之下降，影响了勘探深度的提高，为此需通过增加源的发射时间(便于加大垂直叠加次数)和增加拟地震水平叠加次数来弥补因增加频率成分而导致的信号强度损失，以达到既提高分辨率，又提高探测深度的双重目标.

由本文的方法得到的大地脉冲响应与理论值吻合，合成信号与野外信号非常接近，因此，本文方法所产生的中间结果及理论信号合理可靠，可用在后续的数据处理中，并可用于检验后续数据处理方法的可靠性和有效性.