1 频率域比值算法原理
1.1 均匀半空间
均匀大地表面单位水平电偶极子(偶极长度为1 m，发射电流为1 A)，其轴向电场的阶跃响应(Wait,1960)为
(1)
其中，μ是大地磁导率，一般取真空磁导率4π×10-7，单位为H/m；r为偏移距，单位为m；ρ是背景电阻率，单位为Ωm，erf为误差函数.(1)式对变量t求导，可以得到大地脉冲响应：
(2)
进一步对(2)式求导，令导数为零，求得峰值时刻：
(3)
其中，tpeak，r表示偏移距为r时大地脉冲响应的峰值时刻，与偏移距的二次方成正比，与背景电阻率成反比.
利用(3)式定义视电阻率(Ziolkowski et al.,2007)，假设均匀大地模型的脉冲响应与记录点计算的大地脉冲响应具有相等的峰值时刻，则把假设的均匀大地模型的背景电阻率定义该记录点的视电阻率，计算公式为
(4)
其中，ρa为视电阻率，r为偏移距.野外观测时，视电阻率记录点位置为收发电极中点的连线的中点，与传统偶极-偶极电阻率测深的记录点位置一致，该记录点的偏移距为收发电偶极子中点连线的长度.对于每个记录点，先计算完整大地脉冲响应曲线，然后求取峰值时刻，利用(4)式计算电阻率.
但是，上述方法需要求出完整的大地脉冲响应，而频率域比值法可以利用单个频点值计算峰值时刻，选择较低频率、信噪比较高的频点计算视电阻率.具体原理如下：
利用峰值时刻定义归一化时间参数：
(5)
大地脉冲响应(2)式可以改写为
(6)
其中A为常数因子.如果定义标准大地脉冲响应
是峰值时刻为1 s与峰值为1的大地脉冲响应：
(7)
结合式(5)、(6)、(7)，任意偏移距及任意背景电阻率下一般均匀大地脉冲响应(2)式可以进一步写成如下形式
(8)
从(8)式可以看出，在不同偏移距与不同背景电阻率条件下，均匀半空间的大地脉冲响应形态是一致的，只需要通过伸缩变换和尺度变换就可以与标准大地脉冲响应互相转换.
例如，当均匀半空间的背景电阻率值为20 Ωm，偏移距为2500 m时，大地脉冲响应的峰值时刻为39.2699 ms.首先对该脉冲响应幅值进行归一化尺度变换，即所有幅值除以其峰值使得归一化后峰值变为1，然后对横坐标进行伸缩变换，即时间除以峰值时刻39.2699 ms，从而在相对时间尺度下其峰值时刻变成1.在图1a中，背景电阻率为20 Ωm的均匀大地、偏移距为2500 m的一般大地脉冲响应经过尺度与伸缩变换后如黑色实线所示，黑色十字虚线表示标准大地脉冲响应，二者完全重合，说明一般大地脉冲响应可以通过尺度与伸缩变换转换为标准大地脉冲响应.尺度与伸缩变换的过程如图1b所示，黑色虚线表示一般均匀大地脉冲响应，黑色双向虚箭头表示其幅值与时间范围，首先对其幅值进行归一化尺度变换(除以其峰值)，得到黑色圆点所示的脉冲响应，其幅值范围为0～1，用黑色双向实箭头表示；然后，对黑色圆点所示脉冲响应进行伸缩变换(除以峰值时刻)使其峰值时刻为1，得到黑色实线所示标准大地脉冲响应.当然，这一过程是可逆的，标准大地脉冲响应也可以通过尺度变换(乘以一般大地脉冲响应峰值)与伸缩变换(乘以一般大地脉冲响应峰值时刻)得到均匀半空间下具有任意背景电阻率任意偏移距的一般大地脉冲响应.
图1 (a)均匀大地偏移距2500 m脉冲响应通过尺度与伸缩变换转换为标准大地脉冲响应；(b) 一般大地脉冲响应转换成标准大地脉冲响示意图
Fig.1 (a) Impulse response of homogeneous earth with offset 2500 m is transformed into standard impulse response via time and amplitude scaling; (b) Sketch of how to transform a general impulse response into the standard impulse response
从频率域来讨论这一过程，令标准大地脉冲响应的傅里叶变换为
(9)
由傅里叶变换性质可知一般大地脉冲响应的傅里叶变换为
(10)
由(10)式可知，与时间域类似，一般大地脉冲响应频谱(本文指傅里叶变换)，可以由标准大地脉冲响应的频谱通过归一化尺度变换与伸缩变换获得，伸缩变换因子即为峰值时刻.
除特殊说明外，下文频谱均指归一化频谱，单位为dB，定义如下：
(11)
其中，G(ω)表示幅度谱，G(ω0)表示归一化尺度变换参考频点的幅度谱，在本文中采用幅值最大的频点作为归一化参考频点.为了表述方便，下文称一般大地脉冲响应归一化频谱为一般谱，称标准大地脉冲响应归一化频谱为标准谱，转换或者变换除特殊说明外均指对归一化的频谱进行伸缩变换.
由上文论证可知，一般谱可由标准谱变换得到，而且这个过程是可逆的.因此，均匀半空间情况下同一发射不同偏移距的两个记录点的脉冲响应频谱之间可以通过伸缩互相转换，只是伸缩因子不是某一个峰值时刻，而是二者的比值.仍以背景电阻率为20 Ωm均匀大地模型为例进行说明，图2a—c分别是偏移为500 m频谱、标准谱、偏移距为1000 m频谱的示意图.偏移距为500 m的频谱变换到偏移距为1000 m的频谱有两种方法.第一种方法，先把偏移距500 m的频谱(图2a)进行伸缩，即每个频点的频率乘以偏移距500 m大地脉冲响应峰值时刻，得到标准谱(图2b)，伸缩后偏移距500 m的频谱与标准谱完全重合(图2d)；得到标准谱后，把标准谱进行伸缩，即每个频点的频率除以偏移距1000 m大地脉冲响应峰值时刻，则可由标准谱得到偏移距为1000 m的频谱(图2c)，伸缩后标准谱与偏移距为1000 m频谱完全重合(图2e)；通过上述两步，偏移距500 m频谱转换为偏移距1000 m频谱.第二种方法，直接由偏移距为500 m的频谱(图2a)通过伸缩得到偏移距1000 m的频谱(图2c)，偏移距500 m频谱每个频点的频率乘以偏移距为500 m大地脉冲响应峰值时刻与偏移距为1000 m大地脉冲响应峰值时刻的比值，与偏移距1000 m的频谱完全重合(图2f).偏移距500 m频谱直接变换为偏移1000 m频谱的伸缩因子为两个峰值时刻的比值.事实上，均匀半空间情况下，任意两个不同偏移距的一般谱可以通过伸缩直接互相转换，伸缩因子为二者峰值时刻的比值
图2 (a)偏移距500 m频谱；(b) 标准谱；(c) 偏移距1000 m频谱; (d) 偏移距500 m频谱伸缩变换为标准谱；(e) 标准谱伸缩变换为偏移距1000 m频谱；(f) 偏移距500 m频谱直接伸缩变换为偏移距1000 m频谱
Fig.2 (a) Spectrum at offset 500m; (b) Standard spectrum; (c) Spectrum at offset 1000 m; (d) Spectrum at offset 500 m is transformed into standard spectrum using time scaling; (e) Standard spectrum is transformed into spectrum at offset 1000 m; (f) Spectrum at offset 500 m is directly transformed into spectrum at offset 1000 m
(12)
其中，下标表示α12是频谱1变换到频谱2时频谱1每个频点的频率乘以的伸缩因子.例如从偏移距500 m频谱直接变换到偏移距1000 m频谱，则每个频点的频率乘以的伸缩因子为tpeak，500/tpeak，1000=1.5708(ms)/6.2832(ms) =0.25.当然，这个转换也是可逆的.
均匀半空下两个一般谱可以相互转换，伸缩因子为二者峰值时刻的比值.反之，如果已知二者频谱变换的伸缩因子，就可以求得二者峰值时刻关系.通过乘以伸缩因子α12，频谱1转换成频谱2，相等幅值对应的频率f1乘以伸缩因子后与f2相等，有
(13)
仍然以背景电阻率为20 Ωm的均匀半空间为例，如图3a所示，实线与虚线分别表示偏移距为500 m与1000 m的频谱.任取一个幅值，例如-1 dB，偏移距为500 m频谱与偏移为1000 m频谱对应的频率分别为18.3209 Hz与4.5802 Hz，根据(13)式，则可以得到伸缩因子为0.25，这与由峰值时刻计算得到伸缩因子一致.这样，由500 m偏移距与1000 m偏移距频谱计算得到伸缩因子，利用伸缩因子与偏移距500 m大地脉冲响应峰值时刻可以计算偏移距1000 m大地脉冲响峰值时刻.
可以选择多个幅值，得到多个伸缩因子求平均值.如图3a中圆圈所示，分别在偏移距500 m与偏移距1000 m的频谱从-1 dB到-7 dB每1 dB取一个频点，求得两条谱线对应频点.具体数值列在表1中，第一列是抽取的频点幅值，用A表示；第二列是偏移距500 m频谱对应某一幅值的频率，用f500表示，下标表示偏移距为500 m；第三列为偏移距为1000 m频谱与500 m频谱相等幅值对应的频率，用f1000表示；第四列为对应的频率比值f1000/f500，仍用α12表示；第五列是频率比值α与理论值之间相对百分误差.可以看出，偏移距为500 m与1000 m频谱相等幅值对应的频率比值是一个固定值(如表1与图3b所示)，这个定值就是两个频谱互相转换的伸缩因子，再次说明两个频谱可以通过伸缩互相转换.
综上所述，总结频率域比值算法，计算流程如图4所示，首先，计算参考频谱1与一般谱2，其中参考频谱1峰值时刻已知，对于均匀半空间，参考频谱可取标准谱，对于其他非均匀半空间模型，参考谱也可能是另一个一般谱(1.2节).然后，按幅值Ai分别在频谱1与频谱2抽取第i个频点，计算二者的频率比值(即伸缩因子)：
图3 (a) 偏移距500 m与1000 m频谱，选择频点幅值为-1 dB到-7 dB；(b) 频率比值
Fig.3 Spectrums at offset 500 m and 1000 m，with amplitudes from -1 dB to -7 dB selected;(b) Ratios of selected frequencies
表1 一般谱等幅频点及频率比值
Table 1 Frequencies of general spectrums with the same amplitudes and ratios of these frequencies
图4 频率域比值法计算视电阻率流程图
Fig.4 Flowchart of calculating apparent resistivity with frequency-domain ratio method
(14)
其中，
与式(13)意义类似，但是表示由频谱2(一般谱)变换到频谱1(参考谱)的伸缩因子，与分别表示幅值为Ai谱1与谱2的频率；求得第i个伸缩因子后，抽取下一个频点，直至N个，计算N个伸缩因子的平均
(15)
根据(12)式求得频谱2的峰值时刻
tpeak，r2=α21tpeak，r1,
(16)
根据(4)式即可求得一般谱2的视电阻率.
1.2 非均匀半空间
在1.1节论述了均匀半空间的频率域比值算法理论基础与算法实现.频率域比值算法要求不同的一般谱(包括标准谱)之间可以相互转换，转换的伸缩因子为二者幅值相等频点对应的频率的比值.上文在均匀半空间中严格证明了这点.对于非均匀半空间，一般而言，一般谱之间不能通过简单的伸缩而互相转换.由(14)—(15)式可知，两个一般谱等幅频点的频率比值决定了二者变换的伸缩因子.在均匀半空间中，伸缩因子为常数；对于非均匀半空间，伸缩因子不再是一个常数，而且与实际模型密切相关.但是，对于一维层状介质模型，伸缩变换因子虽然不是一个常数，但是在一定频率范围内变化不大，则仍然可以认为近似不变.陆上MTEM勘探，测点之间的点距一般不超过100 m，如果大地横向地电性质变化不是很剧烈，则可以近似为一维介质，相邻近记录点的一般谱具有相似特征，可以假设二者能够通过伸缩进行转换.若该假设成立，对于一维模型，如果已知一个一般谱及其峰值时刻，则可以通过频率域比值法求得邻近记录点一般谱对应的峰值时刻.本文通过正演模型与野外实测数据来验证这一假设.
2 算例
2.1 均匀半空间模型
为了验证方法的可行性，首先考虑均匀半空间模型.设均匀半空间的背景电阻率值为100 Ωm，偏移距为2500 m，其理论大地脉冲响应峰值时刻为7.8540 ms.输入电流为10阶的m序列，码元频率为512 Hz，电流大小1 A，采样率为16384 Hz(图5a)；通过输入电流与理论大地脉冲响应的卷积合成输出电压响应(图5b).应用频率域比值法，首先计算参考谱，由于是均匀半空间，选择标准谱作为参考谱；然后，由输出电压响应与输入电流的频谱计算偏移距为2500 m的一般谱，与标准谱进行对比，抽取相等幅值频点.如图5c所示，实线表示标准谱，虚线表示偏移距2500 m的一般谱.分别抽取幅值为-3、-6、-10、-15、-20、-30 dB的六组频点，得到标准谱频点和一般谱频点，每组频点的比值可以计算一般谱的峰值时刻，从而算出视电阻率值.抽取的六个频点统计在表2中，由于标准谱峰值时刻为1，因此数值上伸缩因子等于一般谱峰值时刻，为了表示方便，α的单位为10-3，对应峰值时刻单位为ms，抽取的6组频点计算的峰值时刻与理论值7.8540 ms的相对百分误差小于0.03%(表2第5列).由峰值时刻计算的视电阻率值ρa(Ωm)及其与理论值100 Ωm的百分误差分别列在表2的第6列与第7列，可以看出，视电阻率误差小于0.03%.
图5 均匀半空间： (a) PRBS发射电流；(b) 电压响应；(c) 偏移距2500 m频谱与标准谱
Fig.5 Homogenous haft space: (a) PRBS input current; (b) Voltage response;(c) Spectrum at offset 2500 m and standard spectrum
表2 一般谱与标准谱等幅频点及其频率比值
Table 2 Frequencies of general spectrum and standard spectrum with the same amplitudes and ratios of these frequencies
2.2 三层油储模型
图6 三层油储模型(根据Ziolkowski等，2007)
Fig.6 Three-layer model with hydrocarbon reservoir (Replotted from Ziolkowski et al.，2007)
为了验证频率域比值算法适应于一维介质，研究一个三层油储模型，如图6，该模型油层厚度为25 m，顶界面埋深为500 m，电阻率为500 Ωm，背景电阻率为20 Ωm.正演模拟使用empymod(Werthmüller,2017)软件包，发射源为水平接地电偶极子源，接收点偏移距从500 m到12000 m，接收点距为100 m.典型脉冲响应频谱如图7所示，当偏移距为500 m时，三层油储模型响应频谱(实线)与均匀半空间背景频谱非常接近(虚线)；当偏移距为2400 m，可以看到三层油储模型响应频谱(实线)与均匀半空间背景频谱(虚线)有显著差异，油层反映明显；当偏移距为5400 m，三层油储模型响应的频谱(实线)仍然与均匀半空间背景频谱(虚线)有差异，但是差异开始变小.偏移距为500 m的一般谱可认为与20 Ωm均匀半空间背景一般谱相同，从而可以根据(3)式计算其峰值时刻为1.5708 ms.为了计算偏移距600 m的峰值时刻，选择500 m频谱作为参考谱，按照图4所示算法，挑选的频点分布及伸缩因子如图8所示，具体数值见表3.可以看出，正如猜想的一样，二者的频谱在一段频率范围内(500 m频谱频率为19.9566～71.5346 Hz，600 m频谱频率为13.8917～49.5059 Hz)伸缩因子近似为常数，随着频率进一步增高，伸缩因子开始缓慢增大.由表3中第3列可知，伸缩因子的平均均值为1.4578，标准差为0.0307.根据(16)式得到600 m偏移距峰值时刻为500 m偏移距峰值时刻1.5708 ms与平均伸缩因子1.4578之积，为2.2900 ms.根据(4)式求得其视电阻率值为20.0 Ωm，与背景真实电阻率相同，这是由于偏移距较小，油层响应可以忽略.依照上述方法类推，可以由近及远得到不同偏移距视电阻率值，结果如图9所示.图中黑线虚线表示均匀半空间模型的背景视电阻率，与真实电阻率一致；图中的实线是通过计算完整的大地脉冲响应求得的峰值时刻视电阻率，点实线为频率域比值法求得的视电阻率.二者在小偏移距时都接近均匀半空间视真实电阻率，随着偏移距的增加，油层的影响开始体现出来，视电阻率增大，随着偏移距的进一步增加，油层影响又变小，视电阻率变小，频率域比值法计算的视电阻率在偏移距较大时完全趋近真实基岩电阻率.尽管视电阻率数值与真实电阻率之间存在差异，但是作为初步了解地下电性结构的手段，两种方法计算的电阻率都能够反映了真实电阻率低高低的三层特点.
图7 三层油储模型典型响应频谱
Fig.7 Typical spectrum responses of three-layer model with hydrocarbon repository
图8 (a)偏移距500 m与600 m频谱，选择频点幅值由-1 dB到-7 dB；(b)频率比值
Fig.8 (a) Spectrums at offset 500 m and 600 m with amplitude from -1 dB to -7 dB selected;(b) Ratios of selected frequencies
图9 三层油储模型视电阻率
Fig.9 Apparent resistivity of three-layer model with hydrocarbon repository
表3 偏移距500 m与600 m等幅频点及频率比值
Table 3 Frequencies of spectrums of offset 500 m and offset 600 m with the same amplitudes and ratios of these frequencies
2.3 野外实测数据
为了进一步验证本文频率域比值方法的可行性，对野外数据进行了处理.发射偶极位于0～240 m，接收偶极分别位于1320～1380 m与1380～1440 m，偏移距为发射与接收偶极中点距离，分别为1230 m与1290 m，两个测点相距60 m.发射电流为±22 A、12阶、基频4096 Hz的m序列.发射电流实际波形如图10a所示；两个测点的电压响应分别如图10b与10c所示；采用Wiener滤波法计算的脉冲响应如图10d所示，实线表示偏移距为1230 m的脉冲响应，虚线表示偏移距为1290 m的脉冲响应；二者的频谱如图10e所示；选择偏移距为1230 m的频谱作为参考谱，通过频率域比值法计算的伸缩因子如图10f所示.频率域比值法抽取的频点列在表4中，伸缩因子平均值为1.1896，均方差为0.0041，如果把其均方差与平均值的商定义为误差，可以得到伸缩因子误差为0.35%.
根据Wiener滤波法，偏移距为1230 m脉冲响应的峰值时刻为0.7105 ms，视电阻率为267.60 Ωm；偏移距为1290 m脉冲响应的峰值时刻为0.8395 ms，视电阻率为249.09 Ωm.应用频率域比值法，把偏移距为1230 m的频谱作为参考谱，平均伸缩因子为1.1896，根据(16)式计算偏移距为1290 m的频谱对应峰值时刻为0.7105 ms×1.1896=0.8451 ms，与Wiener滤波计算的峰值时刻0.8395 ms相差0.67%.频率域比值计算的视电阻率为247.45 Ωm，与Wiener滤波计算的视电阻率249.09 Ωm一致.
表4 偏移距1230 m与1290 m等幅频点及频率比值
Table 4 Frequencies of spectrums of offset 1230 m and offset 1290 m with the same amplitudes and ratios of these frequencies
3 结论
本文根据多通道瞬变电磁法峰值时刻视电阻的物理意义给出了一种频率域比值方法计算视电阻率.频率域比值法利用低频段几个孤立频点，估计大地脉冲响应的峰值时刻，计算视电阻率，这样可以不用计算整个大地脉冲响应就给出峰值时刻视电阻率.视电阻率可以对地电结构作初步解释.通过均匀半空间与三层油储模型的模拟数据以及野外实测数据说明了该方法的可行性.
图10 (a) PRBS 输入电流；(b) 偏移距1230 m处电压响应；(c) 偏移距1290 m处电压响应；(d) 归一化的Wiener滤波计算的大地脉冲响应；(e)偏移距1230 m与偏移距1290 m频谱；(f)频率比值
Fig.10 (a) PRBS input current; (b) Voltage response at offset 1230m; (c) Voltage response at offset 1290 m；(d) Normalized impulse response calculated by using Wiener filtering; (e) Spectrums at offset 1230 m and offset 1290 m; (f) Ratios of frequencies
参考谱的选择与模型有关，对于均匀半空间，或者短偏移，或者背景电阻率相对高阻的情况，可以选择标准谱作为参考谱.对于一维模型或者横向地电性质变化小的介质，可以选择邻近记录点的一般谱作为参考谱，因此，参考谱是变化的.
