海底电磁接收机的通道标定计算程序

乔玉，陈凯，阳琴

(1.中国地质大学(北京) 地球物理与信息技术学院,北京 100083；2.中国电子科技集团公司 第五十八研究所，江苏 无锡 214000)

0 引言

海底大地电磁方法是通过测量海底大地电磁信号，经数据处理及反演得到海底以下地层的电性信息，进而获取海底深部电性结构。20世纪末，Steven等[1]已设计出可在1000 m深水区获得良好的电磁响应的海底采集系统。经过近些年的发展，在设备研制、资料处理方法上取得了长足进展，但国外在设备研制方面一直对我国进行技术封锁。为了打破国外在OBEM设备方面的垄断，国内同行从1998年开始进行海洋大地电磁相关的方法理论、仪器设备、资料处理、反演解释等研究工作，近年来获得了显著的成功[2-7]。中国地质大学(北京)在多期“863计划”课题及国家专项的支持下，建立了基于ARM 的海底大地电磁数据采集系统——OBEM-Ⅲ，初步实现了低功耗、小型化以及智能化[8]，在油气勘探、水合物调查、深部构造研究等领域取得了一定的成果[9-10]，仪器的主要技术指标已经与以美国SIO为代表的国际先进同行齐平。在国内外同行的推进下，海底大地电磁方法技术在海域油气勘探[11]、地下水调查[12]等资源环境领域和洋脊扩张、板块俯冲[13]、海底火山[14]等构造地质研究领域取得了显著的成果。

在野外测量之前[2]，为确认仪器性能并校正测量数据中包含的仪器通道响应，通常需要对仪器开展通道标定。通道标定的过程包括产生标准方波并采集存储和标定计算两个过程，涉及硬件电路控制以及软件数据处理：① 在硬件电路方面，通过相关命令控制集成在采集舱内部的标准信号发生模块产生3组不同频率的方波，并送至采集通道的最前端进行采集存储，生成通道标定文件(简称标定文件)；② 在数据处理程序中，读取标定文件中的数据，将该数据与原标准方波(标定信号输入)进行比较计算，得到通道的标定结果，比较计算的过程，称为标定计算。在无需外接任何信号发生器的情况下，由标定结果可判断仪器工作状态，并为后期资料处理提供通道响应改正数据。

OBEM-Ⅲ是基于ARM-Linux平台开发的数据采集系统，但现有的标定计算需借助Windows平台开展，每次标定计算前需将标定文件下载到本地计算机，再通过Matlab程序进行计算，存在用户操作繁琐、作业效率低等缺点。为提高海上作业效率，亟需开发一个标定计算程序，并将其封装为可在ARM-Linux平台上运行的APP，可在海底电磁接收机本地端自主完成标定计算。

1 硬件电路原理

标定过程的硬件电路原理框如图1所示，主要由选择开关、采集电路、标定信号产生电路、标定计算电路组成。

选择开关为单刀双掷继电器，可通过微控制单元(MCU，microcontroller unit)控制采集电路接入信号，通道标定时切换到标准信号产生电路，实际信号测量时切换到外部传感器信号输入。

标定信号产生电路由复杂可编程逻辑器件(CPLD, complex programmable logic device)、晶振、斩波电路以及衰减电路组成。标准信号的产生过程如下：CPLD 对晶振产生的信号进行分频，输出峰峰值为3.3 V的单极性方波，再经过斩波电路和衰减电路，得到峰峰值为100 mVpp 的双极性方波，即标定信号。斩波电路运用比较器的原理，将CPLD输出的单极性方波转换为与其频率一致的双极性方波，通过改变CPLD对晶振信号的分频系数，即可输出不同频率的标定信号，晶振的频率准确性与斩波电路的性能，直接影响产生标定信号的准确性。

采集电路包括运放(amp, amplifier)、模数转换器(ADC，analogy-to-digital converter)、CPLD、MCU、GPS、SD 卡等，对输入信号进行采集、存储。电路进行通道之前，由MCU获取当前GPS时间并记录，同时仪器本地时间以此为准继续运行。通道标定开始后，由标定信号产生电路生成的信号通过运放进行滤波、放大，经ADC进行模数转换、 CPLD 完成数据整合，最后由 MCU将采集到的数据写入SD卡内，生成对应的标定文件。后期用Matlab 对采集到的标定文件进行复原，其在时域上表现为方波且无明显畸变，幅值准确性达±1%；频域上基本满足标准方波的频谱分布，相位精度小于±0.37152°。

标定计算电路包括MCU 和SD 卡，MCU 读取SD卡中的标定文件，通过在频域上与标准信号进行比较计算，输出标定结果文件。

图1 标定原理框Fig.1 Block diagram of calibration

2 计算原理

标定计算程序以二进制原始时间序列为输入，以十进制标定计算结果为输出，程序设计包括3部分：读取标定文件、对标定数据进行频谱分析、比较计算得到输出结果文件。

2.1 读取标定文件

标定文件包含TSH、TSM、TSL这3组数据，分别对应不同的采样率以及标定信号频率，详细的参数见表1，其中基频代表输入标定信号的频率，保留谐波数是频谱分析时保留的最高次谐波。在进行数据读取之前，需判断输入标定文件的类型并设置对应的参数。

表1 标定文件的部分参数

每个标定文件包括16 B的块头以及1MB的数据体。块头中包含数据块 ID 和相对时间δt，ID 表示当前数据块的编号，δt代表与对钟时间相间隔的时间。数据体中包含8CH×32768个采样点的数据，每个采样点数据占4 B，以小端模式存储，组织形式为：

L+M+H+channel_flag,

(1)

式中：L,M,H分别对应低、中、高位数据，channel_flag为通道标志位。单通道数据的实际值可表示为式(2)，单位：V：

(2)

2.2 对标定数据进行频谱分析

从标定文件中读取出来的数据，均被放置在有限长的数组中，等待进行下一步的数据处理。标定计算是主要完成对采集信号的频谱分析。对于有限长时间序列，通常采用离散傅里叶变换(DFT , discrete fourier transform)进行频谱分析。

对于N点序列x(n)，其DFT变换的定义为：

(3)

0≤k≤N-1(4)

从式(4)可以看出，N点DFT计算需要N2次复数乘法运算、N(N-1)次复数加法运算。因为实现一次复数相乘，需要四次实数相乘和两次实数相加，所以复数乘法的计算次数直接影响程序的运行时间，特别是当N很大时，计算量呈指数增加。由于需要处理的数据较多，且程序的目标运行平台计算速度不高，直接套用式(3)，运算效率低，耗费时间多，不利于实时对数据进行处理。为了满足现场作业的需求，选用快速傅里叶变换(FFT, fast fourier transform)对数据进行处理。

2.2.1 混合基FFT

FFT计算(又称蝶形计算)是式(3)经过变形之后的图形表达，其结构可清楚的表现数据之间的依赖关系。每完成一组蝶形运算，即可释放输入值的存储空间给输出值，实现同址计算，减少存储空间消耗。常见的运算对有基2、基3等，表2为长度为N时FFT 计算与DFT 计算复数乘法的运算量比较情况，由表可以看出，随着N的不断增加，与直接DFT计算相比，固定基算法计算量的优势越来越明显，但常见的固定基FFT运算只能用于处理长度为Qk(Q为基底数，k为整数)的数据，由表1可知，3 种文件的FFT计算长度并不能满足这一要求。在Windows系统上，可通过末位补零法，增加数据长度，使其达到固定基FFT对输入序列的长度要求但相同的方法在ARM-Linux 平台上会大大增加运算时间。

混合基FFT运算是固定基FFT运算推广后的一般情况，受到计算点数的限制更少，可以对点数为非质数的DFT 进行快速计算，是标定计算的最佳选择。本程序中按照时间抽取的方法，先将X(k)逐层展开到x(n)，将x(n)进行相应的排序后，再从x(n)逐层计算回X(k)。

表2 DFT与FFT复数乘法运算量比较

2.2.2 多基多进制排序

由于对x(n)需进行逐次抽取，必须对输入序列倒位序，得到的次序相当于是原始序列次序的码位倒置，这一过程称为多基多进制的排序。以二进制数为例，任意整数N都可表示为：

(5)

式中:nr,nr-1,…,n2,n1均小于2，等号左边为二进制表达方式，括号外的下角标2代表每位都以2为基底；等号右边为其对应的十进制数表示，经过倒序后，其新的位序为(n1n2…nr-1nr)=n1×2r-1+n2×2r-2+…+nr-1×21+nr×20，新的位序也代表了其新的输入顺序。在多基多进制的情况下，FFT计算长度N可以表示为：

(6)

式中nr,nr-1,…,n2,n1分别小于xr,xr-1,…,x2,x1，等号左边为多进制表达方式，括号外的下角标xr,xr-1,…,x2,x1分别代表每一位的基底；等号右边为其对应的十进制数表示。同理，其倒叙后的新位序N’表示为:

(7)

倒序之后的序列按照对应的十进制大小从小到大进行排序。经过排序后的新序列逐层抽取计算，便可得到最终的FFT计算结果。由于旋转因子WNnk具有周期性，计算量可进一步缩小。图2为N=12=3×2×2时的逐层抽取的过程，其中x’(n)为经过倒序后新序列，X11(k)为x’(n)经过第一次基3抽取的结果，X1(k)序列为X11(k)经过基二抽取的结果，X(k)为X1(k)经过基二抽取得到的最终结果。图3为整个混合基FFT的流程框，其中index2, index3, index5分别为2,3,5的幂，即混合基计算的长度LEN=2index2×3index3×5index5，基2FFT、基3FFT与基5FFT计算采用式(3)进行。

2.3 比较计算

对标定文件完成FFT计算之后，需与输入的标定信号进行比较分析得到通道输入响应。在程序中定义幅值分别为1250以及1075000的理想方波来代替输入的标定信号，其长度与FFT计算长度相同，也进行同样的频谱分析。进行完FFT计算后的理想方波与标定文件在相同频点处进行幅值相除、相位相减，再存储到SD卡中，便生成最后的结果文件。

图2 12点混合基FFT流程Fig.2 12-point flow chart for Mixed-radix

图3 混合基FFT程序框Fig.3 Block diagram for Mixed-radix

3 室内测试

编写好的程序通过交叉编译后上传至SD卡，输入相关命令，即可在OBEM-Ⅲ进行本地标定计算。图4为本地标定计算的实际过程，从计算开始到结束只需11 s，计算结果存储到SD卡的新建文件中,即图4中的PM21L.asc文件。

图5为标定计算结果的幅值相位图，由图可以看出，电道与磁道在通道增益以及相位方面的一致性较好，标定计算的结果有效可靠。

图4 标定计算过程Fig.4 Progress of calibration

4 海上实测结果

2020年7～8月，在南海西南次海盆海底大地电磁探测科学任务中，投入了23台海底电磁接收机，工作水深2 508～4 443 m，完成了46站位海底大地电磁采集任务。海上任务开展前，均完成了通道标定工作，图6给出了现场经预处理得到的测深曲线(S33站位)，有效观测频段为30～30 000 s，获得了高质量海底大地电磁测深数据。

5 结论

针对海底电磁接收机标定计算程序依赖于 PC 端的 Matlab 平台的不足，开展了基于 ARM-Linux 平台的计算程序开发，设计最优计算参数，获得了通道频率响应，得到了标定文件。经室内测试验证了实测结果的正确性，计算耗时约为11 s。该计算程序应用于海上数据采集任务中，获得了高质量海底大地电磁测深数据，有效提升了海上现场作业效率。