在人文活动密集的城市环境下,工业及居民用电的工频干扰往往极为强烈,严重制约了电磁勘探在城市环境中的应用^[1⇓-3]。因此,深入研究和了解城市环境下背景噪声的特征变得尤为重要。目前,国内外学者针对工频噪声进行了诸多研究,研究方向多集中于对工频噪声进行去除,如Butler等^[4]提出了基于工频干扰时域形态相似性和基于正弦波形拟合的两种方法来去除工频噪声及其他周期性噪声;Tang等^[5]提出了基于频域稀疏分解的方法来压制大地电磁数据中的工频干扰;Larsen等^[6]提出了基于机器学习的方法来去除包含瞬变电磁法数据在内的多种物探数据中的工频干扰,均取得了较好的应用效果。然而在实际工程中,环境的多样性对方法的泛用性提出了巨大挑战,需要提出更加泛用且易行的方法来应对工频噪声的影响。

为解决这一问题,实现高质量有效信息的获取,急需开展提高勘探数据的质量和准确评估数据质量两个方面的研究。其一是应用高阶伪随机信号作为发射信号,可以在提高工作效率的同时延长每个测点的数据采集时长,得到更大的数据量和更强的抗干扰能力,其在城市环境下的多次应用均取得优良勘探效果,本次研究所用的就是来自此类信号在世界范围内的早期应用数据;其二是针对高阶伪随机信号的频谱特征,首次于实际工作中应用频谱包络评价算法,对数据中所有的潜在有效频率进行逐频率同标准的质量评估,将受到明显干扰的频率进行剔除,从而实现对工频干扰的规避,提取更多高质量的有效信息。

基于上述两个方面的共同作用,本文提出一种将高阶伪随机信号和频谱包络评价算法相结合的有效信息提取方法,进行了实际应用并取得了明显的效果,为城市环境下可控源电磁法(CSEM)探测提供了一种切实可行的数据处理方法,提高了在复杂电磁环境中电磁勘探数据的准确性和可靠性,为城市地质勘探领域的进一步研究和应用提供了丰富有效的电磁基础资料。

研究区内实际采集的背景电磁噪声频谱如图1所示,城市内存在由大量的电力设备开关、电极极化以及游散电流等产生的低频偏移噪声,因此在低频段上会造成明显的频谱振荡^[7-8],其造成的背景噪声能量较强,在时间域上具有明显的噪声形态。此外,由于城市内大量的车辆干扰和游散电流,会产生高频段的通信噪声以及宽频段的高斯噪声^[9-10],这些背景噪声对高频段有效信号会造成较大的影响。

城市环境下的工频噪声通常是诸多背景噪声中对电磁勘探影响最大的噪声,由于我国工业输电网络的频率为50 Hz,所以工业和居民输用电产生的干扰主要为50 Hz工频噪声及其高次谐波^[11],在图1a 的频谱上表现为50 Hz及其谐波位置的高能量尖峰。这些尖峰会给同频率的有效信号带来极强烈的干扰,故在选择信号频率时只能避开这些频率。值得注意的是,实际输电频率常有细微的变动,并非完全稳定于50 Hz,这意味着工频干扰的影响并不局限于50 Hz及其高次谐波。将50 Hz工频部分局部放大后(图1b),可以明显看到其影响范围是以50 Hz为中心向两侧进行辐射,并且在工频高次谐波频率处也同50 Hz处一样,向两侧辐射一定范围。随着用电情况的变化,这一影响范围也会有所不同,因此对邻近这些频率的有效信号频率也能产生干扰。同时,在部分高频段上,工频噪声也会造成显著的影响,从图1c中可以明显看出,在高频部分,工频噪声的高次谐波能量依旧很高。由于工频干扰的自身能量极高,即便是其高次谐波能量也足够大,而高频信号经过大地衰减,实际能接收到的信号能量往往很小,使得有效信号通常被工频噪声所淹没,因此,在高频段工频噪声依旧不容忽视。

综上所述,城市环境内的背景电磁噪声频谱主要由4部分组成:由游散电流、电极极化等原因导致的低频偏移噪声;由无屏蔽措施的输电电缆等导致的工频噪声;在全频段上普遍存在的高斯噪声;由各类原因导致的其他噪声。其中,工频噪声能量最强,且影响范围较广,是对电磁勘探信号综合影响最大的干扰类型。

伪随机信号在电磁勘探中的应用源于20世纪70年代^[12]。在其之前,主要应用的电磁勘探工作方法为变频法与奇次谐波法,前者具有工作效率低、精度难以提升的缺陷;后者具有相邻频点频差固定和谐波电流强度低等弱点^[13-14]。而伪随机信号能够克服上述两种方法的缺点,可以同时向地下发射多个频率,且振幅相近、频差合理。随后在此基础上,又发展出了2ⁿ序列伪随机信号^[15],2ⁿ 序列伪随机信号具有主频在频率轴上呈对数均匀分布的特征,且所有主频能量(振幅)基本一致,能够在实际勘探工作中实现更多频率的同时发射,因此在电磁法勘探中得到了广泛应用。以基频为1 Hz的13频2ⁿ 序列伪随机信号为例,其理论波形及频谱如图2所示,时域波形的特征为幅值的绝对值不变、正负号按照一定规律变化,频谱则具有由13个对数均匀分布的主频(1~4 096 Hz,图中以红色圈表示)。除主频以外,在2ⁿ 序列伪随机信号的频谱中还存在许多能量稍弱的脉冲,即主频的谐波。尽管这些谐波的能量较弱,但在条件合适的情况下,一样能够作为有效频率使用,通过提取2ⁿ 序列伪随机信号中的有效谐波成分,就能够在不增加任何野外工作量的情况下,进一步增加有效频率数目^[14,16]。

尽管2ⁿ序列伪随机信号可以实现多个频率的同时发射,大大提高了勘探工作的效率和精度,然而在实际工作中,仍然需要进行发射信号的更换,以7主频2ⁿ序列伪随机信号(7频波)为例,当需要测量的频率数量为40个时,就需要使用6组发射信号来分别进行测量。此时,使用伪随机信号进行勘探工作就面临与变频法相似的问题,即需要进行发射信号的更换。为了解决上述问题,杨洋等^[17]提出了高阶伪随机信号,相比传统2ⁿ序列伪随机信号,高阶信号的时域波形同样只包含正负幅值,但频谱更加密集(图3a),可以同时发射更多勘探频率,并保持了2ⁿ序列伪随机信号频差合理、主频能量一致的特点,能够更好地满足日益增长的勘探规模和精度需求。以基频为0.25 Hz的3阶伪随机信号为例,其理论频谱如图3b所示,其包含39个能量基本一致的主频(0.25~3 072 Hz,图中以红色圈表示)。

高阶伪随机信号与传统2ⁿ序列伪随机信号相比,其频谱更为密集,勘探主频率和谐波数量更多,这种密集的频谱为工频及其谐波频率附近提供了更丰富的潜在有效频率点,能够在单次测量中捕获更多的频率信息。对比39频高阶伪随机信号(图3)与13频伪随机信号(图2)可知,能明显观察到39频高阶伪随机信号在工频及谐波附近具有更丰富的主频和谐波频率能量分布,有效且显著地提高了工频附近的潜在可用信息量。

在实际勘探工作中,对所采集数据的质量进行评价是十分重要的一环,然而,对于如伪随机信号一般拥有多个主频且频差较大的发射信号而言,使用传统相对均方根误差(RMSE)方法对数据质量进行评价具有下列问题:①若对所有主频采用统一时窗长度进行评价,则由于时窗内低频所含周期数显著少于高频,致使低频信息质量被严重低估,且时窗长度的选取缺乏统一标准;②若对不同主频采用不同时窗长度进行评价,构建统一标准则变得更加困难;③谐波在实际工作中往往具有成为有效频率的潜力,然而由于上述的不稳定性,常常无法准确评估谐波的有效性。

因此,Zhou等^[18]提出了一种基于噪声连续性的数据质量评估方法,通过对采集数据进行特征包络分析,获取噪声在有效频率位置最大影响的估计值,以包络值替代真实值,并通过定义参数“噪声评价数”来代替真实噪信比,从而对数据质量进行更贴近真实情况的保守估计。其包络及评价流程如图4所示,主要包含数据预处理、包络算法、数据质量评估3个部分。

数据预处理:首先对频谱进行修正,将有效频率成分的幅值置为相邻两点平均值,接着通过Haar小波卷积来识别原始时域数据中的阶跃,并利用多项式分段拟合提取数据中的低频趋势成分,随后对剩余部分使用db5小波基进行a阶小波分解。Haar小波构建方式为

式中:L为Haar小波长度;i为自变量。阶跃检测方式为

式中:D为检测结果序列;T为向两侧延拓后的待检测时间序列;D中波峰的位置与T中阶跃的位置对应。

包络算法:对小波分解得到的震荡项进行基于Hilbert变换的解析包络和b次峰值包络,将包络结果与提取趋势项及分解趋势项叠加,得到修正频谱的最终包络线。解析包络获取方式为

式中:s为自变量;i为虚数单位;f(s)为原函数;fA(s)为f(s)的解析包络;H(f(s))为f(s)的Hilbert变换,变换方程为

式中:p.v.为柯西主值;h(t)=1πt为相变函数;τ为计算卷积使用的中间变量。该变化的实质是将f(s)与h(t)进行卷积,将f(s)的所有频率分量相位推迟90°。

数据质量评估:使用包络线在有效频率位置f的幅值与实际接收幅值作比,得到噪声评价数,即可评估噪声的最大影响:

式中:Vf为有效频率f处的包络幅值;Sf为f处的实际接收幅值;kf为f的噪声评价数,该参数可近似表征噪信比,值越小则说明噪声影响程度越小。

相比于传统RMSE评估方法,该方法通过对噪声进行估计,直接评价信号受噪声影响的程度,在信号能量显著弱于同频率噪声时能够得到更加准确的评估结果;同时,该方法直接对完整频谱进行处理和评估,不需要进行时窗选取,有利于构建统一的质量评价标准。

为了检验包络评价方法与RMSE评价方法的效果差异,将两种方法应用于受到较强干扰的野外采集数据,该数据波形及频谱如图5所示。

在实际工作中,真实的信号值是未知的,通常采用以下方案计算RMSE:首先,对数据进行短时傅里叶变换;其后,用所有时窗对应幅值的平均值代替真实值;最后,进行RMSE的计算。此方案的计算结果实际上就是全部时窗所对应幅值的标准差。在此次对比中,分别采用50 s、100 s、150 s、300 s、450 s、600 s作为计算短时傅里叶变换的时窗长度,并将结果关于接收幅值进行归一化(即除以对应频率的接收幅值),以与包络评价方法进行对比,为了体现统一性,均以评价指标小于等于5%作为有效频率提取标准,对比结果如图6、7。

为了便于分析对比,将不同时窗长度下归一化后的RMSE值曲线与噪声评价数曲线的对比,按频率截取3段典型评价指标曲线分别绘制(图6)。从图中可见,不同的时窗长度对RMSE方法的最终评价结果具有显著的影响,时窗长度越大,评价结果(RMSE值)越小,即认为数据质量越高。包络评价算法所得的噪声评价数整体上处于300 s时窗RMSE曲线与600 s时窗RMSE曲线之间,较为接近450 s时窗RMSE曲线,这一对比结果表明,包络评价方法并不存在时窗长度较短导致噪声估计值偏高的问题。同时,由于包络评价方法是对完整数据频谱进行包络求取噪声估计值,而非通过参数统计,因而也不存在因统计量过少而失去统计意义的问题。

图7为上述方法所对应的提取后归一化电场值曲线(为方便观察,通过乘以不同系数对曲线进行了平移)。归一化电场值是实际测量电位差与发射端电流强度的比值,单位为mV/A,在使用高阶伪随机信号的广域电磁法实际工作中,数据处理使用的是归一化电场值而非实际测量值。由图中可见,RMSE方法在150 s及更短时窗提取的有效频率个数明显少于其他情况;300 s、450 s时窗对应的RMSE方法与包络评价方法提取结果相近,但提取出的有效频率个数少于包络评价方法;600 s时窗对应RMSE方法虽提取频率数略多于包络评价方法,但多出的频率中存在异常值(6.5 Hz),导致曲线变得不够平滑,并不能满足对优质数据的基本假设;同时,在所有方法的提取结果中,仅有包络评价方法成功提取出了40 Hz的信息,且该频率信息不为异常值,满足了对优质数据的假设。

对比结果说明,RMSE评价方法的结果严重依赖于所选择的时窗长度,且在采用相同的截断阈值时,过长的时窗可能导致本不能满足质量要求的频率被引入提取结果。若人为地给不同时窗情况下的信息提取设定不同阈值,则可能需要对同一批次数据设定多个不同的评价标准,并不利于体现质量评估的可信性。

另一方面,在总时长为3 600 s的接收数据中,以300 s作为时窗长度仅能将数据分为12段,在此情况下计算的RMSE是否仍然具有足够的统计意义,具有一定的争议性。若接收数据的总时长更短,那么为了保证RMSE的统计意义,时窗选取还需变得更短,这就将会导致在相同提取标准下,最终提取频率数的进一步减少。

综上所述,包络评价方法相比RMSE评价方法具有如下优势:不受时窗长度限制,能够准确地提取有效频率信息;不需要为同一批次数据设定多个标准,有利于构建统一的评价流程。

本次以“济南城市轨道交通8号线一期工程专项勘查广域电磁探测”数据为基础,使用39频高阶伪随机信号作为发射信号。研究区位于济南市章丘区,测线布置情况如图8所示。本文所用示例测点均来自GY3测线,该条测线紧邻交通干道南侧,沿EW向分布,共有90个测点,设计点距40 m,总长度约为3 560 m,测量时电极分别布设于相邻两点,取较小点号记录数据名称(如GY3-1000测点数据实际对应点号1000与1040的中点位置)。

图9为受到较强噪声干扰的GY3-1720测点数据的波形及频谱。由图9b可见,研究区背景电磁噪声特征与第一节中所分析的城市环境背景噪声频谱特征相符,组成成分一致,且与低频偏移噪声相比,强于第一节中的分析案例。

在发送高阶伪随机信号时,其频谱中含有许多谐波,这些谐波中有些在设计阶段被确定为有效频率,剩下的大部分谐波常常都被认为不能作为有效频率来进一步处理或利用。但是,如果将其视作潜在的有效频率,并对其进行包络评价,那么就能直观地得知其受噪声干扰的程度。

图10为本次实测发射信号的波形及频谱图。受野外工作硬件及环境的限制,实际发射的波形、频谱与理论波形、频谱存在一定差异,如图10a中的实际发射信号并非仅由正负幅值构成,图10b中的主频能量随频率增大而衰减。由图10b可见,发射信号中共有48个理论有效频率,但明显可以看出仍存在其他能量较大的谐波,这些谐波也具有成为有效频率的潜力。通过对实际发射信号的频谱施加阈值,可以提取出一系列潜在的有效频率。这一阈值可以随意设定,只要它足够大到具有作为有效频率的潜力。在本次工作中,阈值设定为全部理论有效频率中最小的幅值。将所有幅值大于该阈值的频率标记为潜在的有效频率,即可得到如图10c所示的潜在有效频率标记情况,被标记为潜在有效频率的频率共有700个。值得注意的是,这些潜在有效频率仍需要进行信息质量评估,以决定其最终是否能够用于后续反演,并不能直接视作有效频率。

以研究区GY3测线上的测点实测数据为例,验证方法在潜在有效信息提取方面的有效性,对标记的所有潜在有效频率进行评估,应用结果见表1。

GY3-1960、GY3-2320测点噪声评价数曲线见图11a、图12a,其上各点为各个潜在有效频率对应的噪声评价数,以5%作为提取阈值,可以看出并非所有潜在有效频率都能够作为有效频率使用。提取前后的归一化曲线对比见图11b、图12b,曲线上各点为噪声评价数小于5%的频率对应的归一化电场值,可见经过提取流程后,有效频率数大幅提高,提升比例最高可达159%,归一化电场曲线上的细节更为丰富。更为丰富的细节能够为反演结果提供更多约束,根据不同频率电磁波趋肤深度的差异,得到更高的分辨率上限。同时,提取后的归一化曲线上并不存在异常值,曲线整体形态保持平滑,说明提取出的额外信息能够满足质量要求,可以提供给后续反演解释使用。

为了进一步验证包络评价方法对受到较强干扰数据的应用效果,将该方法应用于受到较强干扰的GY3-1720测点接收数据上,该数据波形及频谱如图9所示。将包络评价方法应用于该数据,噪声评价数曲线及提取前后的归一化电场值曲线如图13所示。经过提取步骤,有效频率数提升34%,在丰富细节的同时额外提取出0.5 Hz作为有效频率,大幅增加了有效信息反映的最大深度(根据趋肤深度与频率的关系δ=1μσπf,0.5 Hz的趋肤深度为1 Hz的2倍,式中δ为趋肤深度、σ为电导率、μ为磁导率、f为频率),可见即使在数据受到较强干扰的情况下,该方法仍能提取出较多的潜在有效信息,进一步证明了该方法的有效性。

除上述展示的测点外,也对其他GY3测线上的测点数据应用了潜在有效信息提取,结果与上述展示案例基本一致。对于受干扰相对较弱的数据,有效频率数提升比例基本可以达到100%以上;对于受干扰较强的数据,有效频率数提升比例也基本能达到20%~50%;对剩余的大部分数据,提升比例则位于二者中间。同时,所有提取出的额外有效信息均满足质量要求,不会向提取结果中引入明显异常值。

1)本文提出了一种应用高阶伪随机信号,并结合频谱包络评价的方法。试验结果证明,采用高阶伪随机信号发射技术,能够实现多频率同时发射和高精度的数据采集,在频率域规避工频噪声的影响。通过实现频谱包络评价算法,不但能精确评估原始数据的质量,还能将原本未被视为有效频率的谐波标记出来并提取信息。通过对比包络评价方法与传统RMSE评价方法,验证了包络评价方法的优势。本文提出的方法能够有效提高复杂电磁环境中电磁勘探数据的准确性和可靠性。

2)在济南城市轨道交通8号线一期工程专项勘查广域电磁探测实测数据的有效信息提取中,提高了有效频率数量(最高达159%),并且能够提取靠近50 Hz及谐波干扰的有效信息,获得了较好的规避及抵抗工频噪声影响的效果,为城市环境下CSEM探测提供了一种切实可行的数据处理方法。

3)为进一步提高城市复杂电磁环境下CSEM电磁数据的准确性,未来工作将侧重于开发先进的信号处理策略,旨在更好地规避噪声干扰、优化发射信号选择、提取更多的潜在有效信息,使得在复杂的电磁环境中获得更清晰、更可靠的CSEM数据。