Aether大地电磁测量系统是由美国 CG 公司(Crystal Globe Geophysical Research & Service, LLC)研发的新一代高频低频同时采集的电磁系统,其特点是全频段采集记录、多通道电磁兼容、高低频同时连续采集等,与MTU系列仪器相比提供了更多的信息、更高的效率,同时也带来了数据量大、处理费时、占据内存等问题,目前只能应用CG公司开发的PrMT软件进行数据处理。PrMT软件主要包括原始时间序列的傅立叶变换、Robust估计以及大地电磁的电场、磁场阻抗谱或倾子信息的提取和转换等功能。想要进行更深入的研究,比如长时间雷暴干扰的去除、功率谱估计算法、强电磁干扰环境下有效信号的提取、大地电磁数据特征等^[1-5],就必须首先进行数据格式的研究。而研究数据格式,就必须要用到高级编程语言与好的编译器,C#编程语言和VS编译器是目前普遍应用的面向对象编程语言和综合编译器^[6-8]。

前人对MTU系列仪器文件的研究比较多,柳建新等^[9]做了V5-2000大地电磁测深仪文件头数据格式研究,仇根根等^[10]基于MATLAB平台对V5-2000大地电磁测深仪TBL文件进行了分析和研究^[10],郭泽秋等^[11]基于MATLAB语言做了SSMT-2000大地电磁资料批量处理技术的研究,雷清等^[12]基于python语言对MTU和V8系列大地电磁测深仪TBL文件进行了研究,贺景龙等^[13]基于C#编程语言做了MTU系列大地电磁测深仪数据文件的研究与应用,杨凯^[14]基于Python编程语言研究了MTU系列大地电磁测深仪时间序列。除此之外,谢维等^[15]基于VC++编程语言研究了GPD大地电磁测量系统的AVG文件数据格式。前人的研究大多都是针对MTU系列大地电磁测深仪器数据文件,偶有针对GDP系列大地电磁测量仪器数据文件的研究,且国内几乎没有对Aether大地电磁测量仪器数据文件的研究,使得研究Aether大地电磁测量系统十分必要且非常紧迫。

Aether大地电磁系统是新一代高频低频同时采集的电磁系统,并且高频低频是全时段采集。目前多使用美国Crystal Globe公司开发的prMT处理软件进行数据的预处理,主要包括原始时间序列的傅立叶变换、阻抗估计以及大地电磁的电场、磁场阻抗谱或倾子信息的提取和转换。

Aether大地电磁测量系统数据文件主要包含原始时间序列A24文件、探头校准apial文件以及标准时间序列atts文件。

A24文件是野外采集的原始数据文件,对于5道(E_x,E_y,H_x,H_y,H_z)采集MT数据而言,除了采集结束的最后一个A24文件大小不确定外,A24文件几乎占据了32 016 kB大小的存储空间。每个A24文件都包含了若干个数据块,每个数据块由数据头和数据体组成,每个数据头占2 048个字节,每个数据体占92 160个字节。根据正常A24文件的大小,可以计算出正常的A24文件包含了348个数据块。另外,最后采集的A24文件也可以根据其大小计算出数据块的数量。

数据头的2 048个字节表示的参数以及详细内容如表1所示。

A24文件的数据体是24位有符号整数记录的时间序列,每个数据占3个字节。对于5道采集MT数据来说,每个数据体包含的数据量为:92 160÷3=30 720,每道数据量为:30 720÷5=6 144。

atts文件是Aether大地电磁仪GeoTSHub程序合成的一种标准时间序列文件。单个测点的若干个A24文件可以合成一个标准时间序列atts文件。atts文件由数据头和数据体组成,数据头占据2 048个字节,剩余字节均为数据体。数据头的2 048个字节表示的参数及详细内容如表2所示。

对于atts文件数据头的通道类型占据54个字节,其中字节1~2是空格符,字节3~6是整数类型,代表通道编号,字节7~10是整数类型,代表极距,11~22位是字符串类型,代表通道名称及探头编号,字节23~30是空格,31~34位是整数类型,代表通道增益,字节43~46位是浮点类型,代表通道方位角,字节47~50位是浮点类型,代表通道接地电阻。

apial文件是Aether大地电磁仪的磁探头校准文件,Aether系列电磁系统配备4种频段不同的磁探头,分别为CMT-01系探头、CMT-03系探头、CMT-04系探头以及CMTS-10系探头。CMT-01系探头频段为0.000 5~10 000 Hz, CMT-03系探头频段为0.001~1 000 Hz,CMT-04系探头频段为1~100 000 Hz, CMTS-10系探头频段为0.000 1~10 000 Hz。

apial文件由数据头和数据体组成,数据头占64个字节,其余为校准数据。校准数据每16个字节为一个频点的校准数据,数据类型均为浮点型。每个频点的校准数据1~4字节代表频率,5~8字节为振幅数据,9~12字节为相位数据,13~16字节为保留占位符。数据头的64个字节代表的参数及详细内容见表3。数据头系列编号包括两类CMT和MTS,类型编号包括低频、高频、E类型、R类型、C类型等。

利用C#编程语言实现Aether电磁系统时间序列数据的读取需要考虑两个重点问题:一个是多个文件的合并问题,另一个是大量数据的读写问题。

由于时间序列数据必须按照时间先后顺序进行排列,所以文件也需要按照先后顺序进行排列,合并也是按照顺序进行合并。考虑到时间序列文件可能会经过如复制、拷贝等改写过程,其改写时间信息可能会与原本写入时间有出入,因此本文所用的方法是将文件名称以冒泡排序的方法进行文件排序^[16]。原始时间序列名称是由仪器号+日期+点线号+顺序号组成,考虑到存储的更新,日期更改,顺序重新计数,因此冒泡排序需要外循环以日期进行排序,内循环以顺序号进行排序,按照排序好的文件顺序,进行顺序读取即可。

大量数据(超过2 GB)读取会增加计算机内存负担,超过计算机内存可能会造成程序卡死现象^[17]。由于Aether大地电磁仪器是连续采样记录的,在长时间或大采样率进行采样时数据量就会特别大。要想快速读取大量数据,就得利用内存映射的办法分块读写数据。本文采用内存映射的办法,固定映射内存大小为20 MB,这样不会造成内存卡顿现象,按照每个原始时间序列文件大小(32 MB)为分块大小进行数据读写,可以有效地提高数据的读写效率。

由于时间序列数据文件每个数据块都有1个数据头,故时间序列数据文件有大量的数据头,根据数据头各个字节的代表意义,建立一系列的数据列表,如下:

List<string> SystemVer = new List<string>();//系统版本,4字节;

List <int> HeaderLength = new List<int>();//数据头长度(16位整数),2个字节;

List<string> BoxNumbers = new List<string>();//仪器编号(字符串),16个字节;

List <int> NetworkAddress = new List<int>();//网络地址,2个字节 24-25;

List <int> NetworkWifi = new List<int>(); // 网络Wifi,4个字节;

List <int> SendPacketID = new List<int>();//发送包ID,2个字节;

List <int> SendHeaderPacketID = new List<int>();//发送头包ID,4个字节;

List <int> SavePacketID = new List<int>();//接收包ID,4个字节;

List <int> FirmwareVersion = new List<int>();//固件版本,2个字节;

List <int> FPGAVersion = new List<int>();//FPGA版本,2个字节;

List<string> GPSPosInfo = new List<string>();//GPS位置信息,32个字节;

List<string> GPSTimeInfo = new List<string>();//GPS时间信息,24个字节;

List<string> GPSsynclk = new List<string>();//GPS同步时钟,8个字节;

List <int> SamplingRateH1 = new List<int>();//采样率上限,4个字节;

List <int> SamplingRateH2 = new List<int>();//单位采样率,4个字节;

List<string> ChannelNumber = new List<string>();//通道号,8通道64字节;

List<string> SensorNumber = new List<string>();//探头号,8通道64字节;

List<string> ChannelDirection = new List<string>();//方向,8通道16字节;

List <int> AzimuthAngle = new List<int>(); //方位角,8通道32字节;

List <int> PoleLength = new List<int>();// 极距,8通道16字节;

List <int> ChannelGain = new List<int>();//增益,8通道16字节;

List <int> SysStatus = new List<int>();//系统状态,8字节;

List <int> SDTotalSpace = new List<int>();//SD卡大小,4字节;

List<string> SiteName = new List<string>();//测点名称,32字节;

List<string> FileName = new List<string>();//文件名称,64字节;

List <int> DataLength = new List<int>();// 数据长度,4字节;

List <int> SPIBuffInt = new List<int>();// ADC位数,4字节;

List <int> SDFreeSpace = new List<int>();//SD卡剩余大小,4字节;

List<string> GPSAltitude = new List<string>();//GPS高程信息,12字节;

List <int> LsbValue = new List<int>();//ADC精度指标,4字节;

List<string> GPSTime = new List<string>();//GPS时间,32字节。

图1为基于C#编程语言编写的A24时间序列数据文件读写功能模块的WPF界面。该模块主要完成的功能有A24多文件的读写、拼接、合并以及数据头的显示,最终将多个原始时间序列文件合成转化为一个atts标准文件。

数据读取模块的核心代码如下:

foreach (var fileName in fileNames)

{ step++;

int speed = (int)Math.Floor(Lengthstep * step);//进度条进度

worker.ReportProgress(speed, fileName);//异步任务

lock (locker)

{//使用 lock 加锁,防止发生源数组长度不足问题

var memoryMap = new MemoryMap();//内存映射

bytes = memoryMap.GetFileContent(fileName);//读取数据到映射内存

int datanum = bytes.Count/datalen;//数据块数量

int templen = datanum * datalen;//完整数据块数量

int lastdatalen = bytes.Count - templen;//最后一个数据块的长度

for (int i = 0; i < datanum; i++)

{ headersarray.AddRange(bytes.GetRange(i * datalen, hd));//数据头

datasarray.AddRange(bytes.GetRange(i * datalen + hd, len));//数据体

if (i == 0 && step == 1)

{ //写入文件

FileHelper.ByteToFile(headersarray.ToArray(), writepath);

}

} if (lastdatalen != 0)

{ //结尾的数据块

headersarray.AddRange(bytes.GetRange(datanum * datalen, hd));

datasarray.AddRange(bytes.GetRange(datanum * datalen + hd, lastdatalen-hd));

}

}

}

时间序列的编辑模块主要是完成原始时间序列的显示、编辑、整块切除、限频滤波等功能(图2)。读取格式为Aether大地电磁系统标准时间序列atts文件,atts文件的数据是已经过格值校准后的数据,格值校准是将单片机采集的原始整数值通过单片机采样位数和采样精度换算为采样电压值的过程^[18]。标准的24位单片机的电压值换算为:

式中:ADC__Value为采样的整数值;V_REF为参考电压;m为参考电压的范围;Gain为增益。

时间序列编辑功能模块核心代码如下:

while (endLength > 0)

{ fileStream.Position = start;//开始读取位置

if (endLength < length)//分块读取数据

{ bytesRead = bufferedStream.Read(array,0,Convert.ToInt32(endLength));

for (int i = 0; i < bytesRead; i += 15)

{ step++;

//各分量格值数据的读取

double dataEx = Int24ToInt32_2(array[i], array[i + 1], array[i + 2]);

double dataEy = Int24ToInt32_2(array[i + 3], array[i + 4], array[i + 5]);

double dataHx = Int24ToInt32_2(array[i + 6], array[i + 7], array[i + 8]);

double dataHy = Int24ToInt32_2(array[i + 9], array[i + 10], array[i + 11]);

double dataHz = Int24ToInt32_2(array[i + 12], array[i + 13], array[i + 14]);

//各分量电压数据的转换

datasEy.Add(dataEy / lsb);

datasHx.Add(dataHx / lsb);

datasHy.Add(dataHy / lsb);

datasHz.Add(dataHz / lsb);

tsTimes.Add(step * timeternal);

}

} else

{ bytesRead = bufferedStream.Read(array, 0, length);

for (int i = 0; i < length; i += 15)

{ step++;

double dataEx = Int24ToInt32_2(array[i], array[i + 1], array[i + 2]);

double dataEy = Int24ToInt32_2(array[i + 3], array[i + 4], array[i + 5]);

double dataHx = Int24ToInt32_2(array[i + 6], array[i + 7], array[i + 8]);

double dataHy = Int24ToInt32_2(array[i + 9], array[i + 10], array[i + 11]);

double dataHz = Int24ToInt32_2(array[i + 12], array[i + 13], array[i + 14]);

datasEx.Add(dataEx / lsb);

datasEy.Add(dataEy / lsb);

datasHx.Add(dataHx / lsb);

datasHy.Add(dataHy / lsb);

datasHz.Add(dataHz / lsb);

tsTimes.Add(step * timeternal);

}

} start += bytesRead;

endLength -= bytesRead;

}

数据的读取采取循环分块的方式,数据的显示采用开源的scottplot控件,由于计算机配置以及scottplot控件的功能限制,量级大于百万级的数据显示操作会延时卡顿,目前国内没有很好的开源项目解决此问题,国外有比较先进的商业技术,比如借助GPU加速的LightningChart图表控件,能够解决亿级数量的显示和操作。

Aether系列大地电磁仪器的磁探头是按照不同的探测频段分类的,文件扩展名均为.apial, 探头的校准文件反映了探头的频点个数以及各个频段数据振幅的相位数据。校准需要归一化处理,对于校准文件里面没有的频点需要进行插值处理。插值方法包含了线性插值和非线性插值等。

图3是基于C#编程语言编写的磁探头校准文件分析模块界面,探头校准文件显示分析模块的核心代码如下:

for (int i = 0; i < (buffer.Length-64)/16; i++)

{ //频率

float fre = BitConverter.ToSingle(buffer, 64 + i * 16);

fre_Cal.Add(fre);

//振幅

float amp = BitConverter.ToSingle(buffer, 68 + i * 16);

amp_Cal.Add(amp);

//相位

float pha = BitConverter.ToSingle(buffer, 72 + i * 16);

pha_Cal.Add(pha * millRadianToDegree);//毫弧度到度的转化

}

本次选取2023年7~8月在新疆和青海实测的5分量大地电磁数据对编写的各个功能模块进行实测数据的应用测试,并与prMT软件进行对比验证,检验程序的准确性。L01P011是在新疆和田地区采集的大地电磁数据,L02P011是在青海果洛州采集的大地电磁数据。数据的观测时间均大于24 h,采样率为2 000 Hz。

图4a、图4b分别是用prMT和作者编写的DataWindow模块对L02P011数据文件读取的第一个数据头信息的对比。可以发现信息基本正确,仅在处理一些版本信息和需要转换的电压、状态等信息格式时存在不同,作者编写的DataWindow模块可以同时读取多个文件的数据头信息,也可以通过切换来查看多个数据头信息。

图5a、图5b分别是用prMT和作者编写的attsWindow模块对L01P011数据文件时间序列傅立叶变换的对比。图5a是PrMT傅立叶变换之后各分量的自功率谱,图5b是attsWindow模块傅立叶比变换之后各分量的自功率谱。为了对比程序傅立叶变换的可靠性,本次对傅立叶变换的参数做了统一,其中窗口大小为81,数据对均为100个,PrMT是否进行过滤波算法未知。从计算的结果和各分量自功率谱的情况来看,prMT程序对100个数据对的处理不是简单的平均,也不是求取中位数,更不是应用了简单的挑谱算法(峰值挑谱、能量占比等),其应该是采用了高通滤波和低通滤波的方式进行了处理,这样的曲线看起来比较平滑。attsWindow模块对100个数据对的处理采用的是高斯滤波法进行谱挑选。从图5的自功率谱曲线对比情况来看,曲线形态基本一致,attsWindow模块的自功率谱曲线细节更加丰富。

PrMT软件对于时间序列的处理可能在滤波算法方面更加优秀,其傅立叶变换和功率谱估计均采用了级联抽样和级联萃取的方法,从软件的运行情况来看,高频数据的谱数量上限为99 999,这是为了节省计算时间,同时节省计算机内存占用。但高频数据采集却是连续的,必然会存在数据冗余。

图6a、6b分别是用prMT和作者编写的attsWindow模块时间序列显示的对比。从结果来看,两个软件显示的时间序列数据形态基本一致。个别点数据大小不一样主要是因为抽稀显示的原因。prMT采用的是百万级数据的显示,即每一道显示1百万的数据,作者采用的是亿级数据显示,即每一道显示1亿的数据点。两个数量级相差100倍,虽然在数据显示上没有太大的区别,但是后期傅立叶变换和谱估计就会有很大区别,因为参与估计的数据量大,就会有更加丰富的频谱信息,特别是低频段谱估计时,更多的估计数据能够提供更加准确的频谱信息。

图7是时间序列显示1 s内的数据量对比,图7a是PrMT软件1 s内显示的数据量,由图中大致可以看出1 s内大概有10个数据,图7b是作者编写的attsWindow软件模块1 s内显示的数据量,由图中大致可以看出1 s内的数据量很多。由于attsWindow模块的显示级别为亿级,抽稀因子为3,数据原始采样率2 000 Hz,经过计算可以得出attsWindow模块1 s内的数据量大致为666个。

数据量的增加会占用大量存储空间,数据量显示的增多也会占用更多内存,本次测试所用数据的时间序列采样率为2 000 Hz,时间长度为89 min。attsWindow使用的scottplot插件不用对数据进行抽样就可以完整显示,内存占用量在1.4 GB左右。降采样显示的确可以降低内存占用空间,但同时会存在混叠效应,对于数据的处理以及干扰的判别会有一定的影响。

本文通过研究Aether大地电磁系统的数据文件,基于C#编程语言编写Aether大地电磁系统数据读写模块、时间序列显示模块、探头校准文件分析模块,通过实测数据的应用,检验程序运行的可操作性,与prMT软件运行结果进行比较,基本验证了几大功能模块运行结果的可靠性。

通过对比prMT软件与作者编写的功能模块的运行结果,发现prMT软件在时间序列的显示上速度较快,但信息量较少,attsWindow模块显示速度较慢,但信息量较大。因此后期的傅立叶变换和频谱分析,两者也不一样。傅立叶变换结果显示,attsWindow模块也有更加丰富的频谱信息。以上对比证实,attsWindow模块能够显示更加丰富的频谱信息,对求取的功率谱未进行圆滑处理,能够对实际的电磁信号有更加全面真实的反映。