Matlab在CSAMT数据“去噪”及近场校正中的应用

2019-11-19胡庆辉

物探化探计算技术 2019年5期

胡庆辉

(1.山东省地质矿产勘查开发局第四地质大队，山东潍坊 261021; 2.山东省地矿局海岸带地质环境保护重点实验室，山东潍坊 261021)

0 引言

可控源音频大地电磁测深法(CSAMT),是在大地地磁测深法(MT)和音频大地电磁测深法(AMT)的基础上发展起来的一种人工源频率域测深方法[1]。自上世纪70年代引入国内，在油气勘查、金属矿普查、水文及环境地质等多个领域显示出广阔的应用前景。随着经济社会发展，CSAMT野外数据采集所面临的“电磁环境”越发复杂，勘探信号不可避免地受到“污染”，导致最终反演数据质量不高，加大了资料解译的难度甚至造成失误。此外，受限于CSAMT方法的原理及假设前提，在中、低频率上电磁场进入“近区”或者“过渡区”，卡尼亚视电阻率将发生畸变，无法实现勘查目的[2]。底青云等[3]认为CSAMT资料处理的质量依赖于近场校正的质量。

国内学者针对CSAMT方法的电磁干扰压制及近场校正技术进行了相关研究。伪随机编码技术在供电环节的使用能够有效提高野外观测数据信噪比[4-5]，而雷达等采用信息熵和有理函数滤波相结合的方式对干扰数据进行处理也取得到比较好的效果[6]。通常的近场校正的方法是“三角形校正法”[3，7-9]，而采用“全频域视电阻率法”进行近场校正得到越来越广泛的关注和应用[10-11]。

用户界面(或接口)是指人与机器(或程序)之间交互作用的工具和方法[12]。笔者采用Matlab GUI编程，实现友好的可视化交互式操作界面，剔除“畸变”频点数据，实现CSAMT数据的“去噪”处理。采用数值计算功能，编写了基于全频域视电阻率计算的程序以实现近场校正。程序的应用提高了CSAMT数据的处理效率和效果，保证后续处理的数据质量，在实际生产中得到应用并取得了很好的效果。

图1 CSAMT野外实测干扰数据典型特征Fig.1 Typical feature of interfering survey data of CSAMT(a)卡尼亚视电阻率曲线；(b)阻抗相位曲线；(c)电场分量幅值曲线；(d)磁场分量幅值曲线

1 CSAMT干扰及“近场”数据特征

CSAMT野外数据环节的电磁干扰主要来自天然电磁噪声和人文电磁噪声[13]，其中天然电磁噪声主要来自雷电活动，其引起的电磁场在CSAMT常用工作频段内能量较高。这类噪声一般可采用多次叠加进行压制。人文电磁噪声主要包括工业用电、大功率发射台以及随机的电磁脉冲干扰(如动力电开关转换、铁磁性物体靠近磁探头等)。固定频率的电磁干扰可以通过设置陷波器进行压制，而能量较高、随机的强干扰还不能完全克服，需要野外施工人员正确判断电磁干扰特点，采取必要的规避措施，尽量避开干扰源或者错开干扰时段，以便采集到相对高质量的野外数据。

随机电磁干扰在CSAMT实测数据中主要表现为电场及磁场分量复测离差增大，单频点数据离差加大，数据一致性较差，阻抗相位偏差较大，甚至出现反相，计算得出的卡尼亚电阻率出现“飞点”，曲线形态不能反映真实大地电阻率的变化特征(图1)。

基于上述两个方面的问题，笔者利用Matlab数值计算及GUI功能，设计了针对野外数据预处理和近场校正的相关程序，在CSAMT数据“去噪”预处理及近场校正等方面进行了探索，取得较好的应用效果。

2 CSAMT“去噪”及“近场校正”方法

现阶段，可控源音频大地电磁方法(CSAMT)“去噪”多采用滤波方法，在干扰比较小的情况下，一般能够得到比较好的效果。但随着干扰的加强，滤波方法往往难以取得理想的电磁噪声压制效果。

本次“去噪”处理思路比较直观，即从多次观测数据中挑选“合理”数据，舍弃畸变数据，从而达到“去伪存真”的目的。所谓“合理”，即从电极布设及采集时就保证观测质量，合理避开干扰源或干扰地段，错开干扰时段进行数据采集，使多数观测数据为真实地电反映，其数据分布应具有“集中分布”特征，重复性较好。而少数观测值由于随机干扰产生“畸变”，相关参数如相位、电场分量幅值、磁场分量幅值等偏离合理区间，出现明显的“飞点”特征，通过友好界面对畸变数据予以剔除，从而实现“去噪”处理。

图2 CSAMT数据预处理GUI界面Fig.2 Matlab GUI interface for data pre-processing of CSAMT

数据处理通过鼠标在界面上框选畸变数据实现便捷操作(图2)。用户界面分为三部分：①工具栏；②绘图区；③信息提示区。工具区包括数据加载、重置、显示参数选项、数据操作、数据保存等功能按键或下拉菜单。绘图区则显示相应选项的散点图、拟合曲线图或者拟断面等值线图。界面下方为信息提示区，显示部分文件信息以及当前点号等。

加载CSAMT数据完成后，绘图区默认显示首个测点的卡尼亚电阻率数据，绿色圆点为某个频率重复观测数据点，蓝色曲线为其样条插值拟合曲线。可以看出，在电磁干扰条件下采集的数据存在明显的离散特征，即重复观测数据一致性较差，“飞点”较多。同时，每个频率的观测数据点有重复性较好或相对较为集中的特征，我们认为这样的数据为“合理”数据。用户使用鼠标选择工具栏“删除”按键，在绘图区框选需要剔除的数据点(图2中红色框及红色点)，它们将不参与后续的数据处理。利用下拉菜单选择“阻抗相位”、“电场分量幅值”及“磁场分量幅值”三个选项在相应数据散点上完成类似操作，实现多个参数对同一个测点的“去噪”处理。当该测点数据选择完毕后，通过“保存”按键实现AVG数据保存，为后续处理提供高信噪比数据。

CSAMT数据近场校正是在上述AVG数据的基础上，利用“等效电阻率全频域视电阻率法”，采用迭代拟合方式计算等效电阻率。

通常采用赤道偶极方式进行CSAMT数据采集时，采集参数为电场和磁场的水平分量Ex、Hy，卡尼亚视电阻率可以表示为式(1)。

(1)

(2)

利用ρs是P的单调下降函数，通过实测ρs确定P值，并通过式(3)计算大地电阻率。

(3)

Matlab数值计算的实现过程中，采用迭代方式求取F(P)=0的根。具体思路为：给定取值范围为[a,b] ，给P赋初值c=(a+b)/2，判断F(a).F(c)符号，根据函数的单调下降性质，采用二分法逐步逼近零点，从而求取P值进而得到大地电阻率ρ。

近场校正处理通过Matlab GUI编程实现友好用户界面进行操作。通过主界面“近场及静校正”按键，可以调用NFC(近场校正)模块(图3)。NFC界面包含工具栏及绘图区两部分，工具栏主要包括数据加载、数据列表、数据编辑、近场校正及数据保存按键。在加载数据环节需要同时加载CSAMT实测数据(*.AVG)、观测点坐标文件(*.stn)，用于计算每个测点的收发距及相对角度，计算结果将在相应绘图区显示。绘图区显示测点卡尼亚视电阻率曲线和阻抗相位曲线，用户可以通过进一步的编辑功能对曲线进行处理，调整曲线形态，使其更符合“非波区”曲线特征。计算完成后则显示近场校正前后曲线对比情况，使处理人员能够初步了解校正效果。

图3 CSAMT数据近场校正GUI界面Fig.3 Matlab GUI interface for near-field correction of CSAMT

3 近场数据正演模拟计算

3.1 均匀半空间地电模型

均匀半空间地电模型电阻率为1 000 Ω·m，正演计算收发距r=6 km及r=14 km两种情况下的卡尼亚视电阻率特征(图4)。可以看出，高频段视电阻率能够正确反映模型电阻率特征，当进入到中低频率时，曲线发生畸变，近45°上扬，严重偏离模型电阻率。采用“等效电阻率全频域视电阻率法”校正后，全频域视电阻率曲线已经基本水平，较准确反映了模型电阻率数值，说明利用全频域视电阻率法进行CSAMT近场校正效果非常理想。

图4 均匀半空间CSAMT正演数据近场校正处理曲线对比Fig.4 Near-field correction curve contrast of CSAMT forward data in the homogeneous half-space model

图5 层状介质CSAMT正演数据近场校正处理曲线Fig.5 Near-field correction curve of CSAMT forward data in the layered model(a)H型测深曲线；(b)A型测深曲线；(c)K型测深曲线；(d)Q型测深曲线

图6 某测线CSAMT 1900点去噪及近场校正处理效果Fig.6 Data denoise and near-field correction process results of survey point 1900(a)实测卡尼亚电阻率数据散点及平均曲线；(b)实测阻抗相位数据散点及平均曲线；(c)“去噪”及近场校正后曲线对比；(d)“去噪”后阻抗相位数据点及平均曲线

图7 CSAMT测量L4线视电阻率断面图(预处理前)Fig.7 Apparent resistivity section of CSAMT survey of line 4(before preprocessing)

3.2 层状地电模型

层状地电模型采用三层模型，分别设置了H型、A型、K型及Q型四种测深曲线类型(图5)。在中低频段，各地电模型正演曲线进入“近区”，曲线形态发生畸变。采用“全频域视电阻率法”校正后，测深曲线形态基本反映了真实的地电模型特征。

同时也应该注意到，在均匀半空间模拟中，过渡区出现局部“下凹”，校正值偏离模型电阻率。对于层状模型，尤其是K型测深曲线，中间高阻层虽然趋势较明显，但其值在校正后明显偏低。对于过渡区“校正不足”现象，笔者也做了探索性工作，主要思路是在中低频段数据近场校正时增加修正系数。根据大量正演模型试算数据，对近区和过渡区近场校正修正系数采用多项式拟合。修正系数公式见式(4)及式(5)。

近区修正系数：

K(f)=-0.0071f3+0.0773f2-

0.2017f+1.0634

(4)

过渡区修正系数：

K(f)=-0.0018f3+0.0111f2-

0.0164f+1.6322

(5)

其中，f为校正点频率。在进行近场校正时，合理选择近区及过渡区频率，并在校正后数值基础上进行修正。该方法在正演数据中得到了很好的校正效果，过渡区“下凹”现象得到改善，近区数据也得到适当调整，在野外实测数据中的应用效果有待进一步验证。

图8 CSAMT测量L4线视电阻率断面图(预处理后)Fig.8 Apparent resistivity section of CSAMT survey of line 4(after preprocessing)

4 处理实例及效果分析

某高速公路连接线拟建路线穿越煤炭采空区，为确保拟建项目的安全，在工作区内开展CSAMT测量工作，查清采空区位置、分布范围等特征。观测装置采用赤道偶极标量方式，采集参数为Ex和Hy，发射偶极AB=1 200 m，测量偶极MN=50 m，最大收发距为6.4 km，测量频段为1 Hz～8 192 Hz。由于工作场地为城市近郊，村庄厂房密集且有高速公路和多条高压线穿过，电磁干扰非常严重。在施工过程中采取多种措施保证野外采集数据质量，如加大供电电流、增加重复观测次数、合理避开干扰时段等。

以某测线1900点的数据散点及卡尼亚电阻率和阻抗相位曲线为例(图6)。图6(a)、图6(b)分别为卡尼亚电阻率及阻抗相位实测数据散点及平均值曲线，由图6(a)、图6(b)可以看出，各个频段的数据点均存在重复性差、离散度高等特点，说明该测点在采集过程中受到较强烈电磁干扰，造成数据信噪比下降。中低频段进入“近区”，视电阻率曲线以近45°上扬，且在32 Hz处出现局部极小值，应为过渡区反映。阻抗相位曲线也显示了近区及过渡区信息。

图6(c)、图6(d)为数据处理后数据散点及曲线特征。利用预处理程序，剔除离散度大的点，保留重复性好、相对集中分布的点作为“合理”的实测数据，曲线形态得到适当调整，数据信噪比得到提高。采用处理后的数据进行近场校正(图6(c)红色曲线)，选取1 Hz～8 Hz为近区频段，16 Hz及32 Hz为过渡区频段，增加相应的修正系数，采用“全频域视电阻率法”进行计算，可以看出曲线形态发生较大变化，低频段呈45°上扬段斜率变小，过渡区电阻率抬升，高频段视电阻率没有发生变化，曲线整体形态更加合理，便于反演解释。

选取L4线进行平滑滤波二维反演，对处理前后的数据进行对比，以检验该程序的处理效果。L4线东西向布设，穿过煤炭采区，沿线多为厂区，并与高压输电线和铁路相交，电磁干扰较强烈。采用原始数据进行反演(图7)，视电阻率等值线“下拉”现象明显，剖面后段约-400 m深度近水平展布的低阻等值线出现波状起伏且埋深偏浅，新太古代泰山岩群片麻岩地层起伏形态也发生改变。尤为明显的是石24井位置埋深-200 m处出现团块状高阻异常区，与钻井揭示岩性不符，出现“假异常”。

经过“去噪”预处理及近场校正，断面等值线形态更加合理，基本反映了真实的地质信息，异常形态符合地质规律，地电层划分较明显，其形态与已知地质断面基本一致。浅部为不连续的高阻特征，推测为青山群八亩地组凝灰质砂砾岩电性反映，局部低阻异常推测为断裂构造发育部位。中部存在较连续水平展布的低阻异常带，推测为淄博群坊子组含煤地层，低阻异常为采空区充水所致。深部则以高阻为主要特征，推测为新太古代泰山岩群片麻岩地层，局部低阻凹槽推测为断裂构造发育部位。后经钻探验证，煤炭采空区实际位置及埋深与物探解译结果基本一致，L4线验证孔钻遇采空区顶板深度与解译深度误差小于10%。