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

doi:10.11720/wtyht.2022.1466

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

乔玉^,¹, 陈凯^,¹, 阳琴²

1.中国地质大学(北京) 地球物理与信息技术学院,北京 100083

2.中国电子科技集团公司第五十八研究所,江苏无锡 214000

Channel calibration calculation program for OBEMs

QIAO Yu^,¹, CHEN Kai^,¹, YANG Qin²

1. School of Geophysics and Information Technology, China University of Geosciences(Beijing), Beijing 100083,China

2. The 58th Research Institute of China Electronics Technology Group Corporation, Wuxi 214000, China

通讯作者: 陈凯(1984-),男,副教授,主要从事海洋电磁仪器开发及应用工作。Email:ck@cugb.edu.cn

责任编辑: 王萌

收稿日期: 2021-08-24 修回日期: 2022-01-26

基金资助:

国家自然科学基金项目“海底MT的运动海水电磁噪声分离方法研究”(41804071)

Received: 2021-08-24 Revised: 2022-01-26

作者简介 About authors

乔玉(1999-),女,主要从事地球物理仪器开发及应用工作。Email: 2010200023@cugb.edu.cn

摘要

海底电磁接收机(OBEM, ocean bottom electromagnetic receiver)是用于观测海底大地电磁信号的主要仪器。为了验证仪器的功能并提升仪器测量精度,在仪器入水前采用通道标定的方式进行仪器自检。现有的OBME-Ⅲ在进行标定计算时,需将通道标定生成的标定文件从仪器中导出,再借助PC端的Matlab程序进行计算,存在海上作业复杂、计算效率低等不足。针对这一系列问题,笔者开发编写了基于ARM-Linux平台的标定计算程序。通过输入相关命令,便可在OBEM-Ⅲ本地实现通道标定计算,减少了导出数据进行多平台计算引起的麻烦。所采用的混合基快速傅里叶变换计算方法,在保证计算精度的同时,将计算时间从90 s缩短到11 s,大大提升了计算速度,提高了仪器海上作业的效率。在2020年7~8月南海西南开展的大地电磁科研任务中,该程序表现出色,成功获得高质量海底大地电磁测深数据。

关键词： 大地电磁; 海底电磁接收机; 标定计算; 混合基FFT

Abstract

An ocean bottom electromagnetic receiver (OBEM) is primarily designed to measure the submarine electromagnetic signal. In order to verify the function of the instrument and improve the accuracy of result, the instrument is self-checked by a built-in channel calibration before OBEM enters the water. OBME-Ⅲ needs to export the calibration file generated by channel calibration from the instrument, and then use the Matlab program on the PC side to perform calibration calculation, which has shortcomings such as complex offshore operations and low calculation efficiency. In order to solve these problems, the author developed a calibration calculation program based on the ARM-Linux platform. By entering the relevant commands, the channel calibration calculation can be realized locally in OBEM-Ⅲ, reducing the trouble of exporting data for calculation on multiple platforms. The mixed-base fast Fourier transform calculation method adopted reduces the calculation time from 90 s to 11 s while ensuring the calculation accuracy, which greatly improves the calculation speed and improves the efficiency of the instrument's offshore operations. In the magnetotelluric scientific research mission carried out in the southwest of the South China Sea from July to August 2020, the program performed well and successfully obtained high-quality submarine magnetotelluric sounding data.

Keywords： magnetotelluric; OBEM; calibration calculation; mixed radix FFT

PDF (1945KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

乔玉, 陈凯, 阳琴. 海底电磁接收机的通道标定计算程序[J]. 物探与化探, 2022, 46(3): 550-556 doi:10.11720/wtyht.2022.1466

QIAO Yu, CHEN Kai, YANG Qin. Channel calibration calculation program for OBEMs[J]. Geophysical and Geochemical Exploration, 2022, 46(3): 550-556 doi:10.11720/wtyht.2022.1466

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所示,主要由选择开关、采集电路、标定信号产生电路、标定计算电路组成。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 标定原理框

Fig.1 Block diagram of calibration

选择开关为单刀双掷继电器,可通过微控制单元(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卡中的标定文件,通过在频域上与标准信号进行比较计算,输出标定结果文件。

2 计算原理

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

2.1 读取标定文件

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

表1 标定文件的部分参数

Table 1 Partial parameters of calibration file

输入文件名称	采样频率 /Hz	基频 /Hz	数据长度	单周期采样点	FFT计算长度	保留谐波数	最高计算频点/Hz
TSH	2400	24	24000	100	24000	22	1032.00000
TSM	150	2.5	1200	60	1200	14	67.50000
TSL	150	15/160	4800	1600	1600	38	7.03125

新窗口打开| 下载CSV

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

(1)

L + M + H + c h a n n e l_f l a g,

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

(2)

v a l u e = \frac{(H \times 256 \times 256 + M \times 256 + L)}{2^{23}} \times 2.5 (V)

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

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

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

(3)

\begin{matrix} X (k) = \overset{N - 1}{\sum_{n = 0}} x (n) e^{- i \frac{2 π}{N} n k}, 0 \leq k \leq N - 1 \end{matrix}

式中:x(n)为输入序列,X(k)为经过DFT计算后的输出;为了书写方便,常用W_N= $e^{- i \frac{2 π}{N}}$ , ${W_{N}}^{n k}$ 被称为旋转因子。根据欧拉公式,式(3)可表示为:

(4)

\begin{array}{l} X (k) = \overset{N - 1}{\sum_{n = 0}} x (n) \cdot [c o s (\frac{2 π n k}{N}) - i \cdot s i n (\frac{2 π n k}{N})], \\ 0 \leq k \leq N - 1 \end{array}

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

2.2.1 混合基FFT

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

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

Table 2 Comparison of DFT and FFT complex multiplication calculations

	DFT/(N为任意自然数)	基2FFT(N=2^k)	基3FFT(N=3^k)	基5FFT(N=5^k)
复数乘法	N²	(Nlog₂N)/2	2N log₃N	4N log₅N
复数加法	N(N-1)	Nlog₂N	2N log₃N	4N log₅N

新窗口打开| 下载CSV

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

2.2.2 多基多进制排序

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

(5)

\begin{array}{l} N = (n_{r} n_{r - 1} \dots n_{2} n_{1})_{2} \\ = n_{r} \times 2_{r - 1} + n_{r - 1} \times 2_{r - 2} + \dots + n_{2} \times 2^{1} + n_{1} \times 2^{0}, \end{array}

式中:n_r,n_r_-1,…,n₂,n₁均小于2,等号左边为二进制表达方式,括号外的下角标2代表每位都以2为基底;等号右边为其对应的十进制数表示,经过倒序后,其新的位序为(n₁n₂…n_r_-1n_r)=n₁×2^r^-1+n₂×2^r^-2+…+n_r_-1×2¹+n_r×2⁰,新的位序也代表了其新的输入顺序。在多基多进制的情况下,FFT计算长度N可以表示为:

(6)

\begin{array}{l} N = (n_{r} n_{r - 1} \dots n_{2} n_{1})_{x_{r} x_{r - 1} \dots x_{2} x_{1}} \\ = n_{r} (x_{r - 1} \cdot \dots \cdot x_{2} x_{1}) + \\ n_{r - 1} \cdot (x_{r - 2} \cdot \dots \cdot x_{2} x_{1}) + \dots + n_{2} \cdot x_{1} + n_{1}, \end{array}

式中n_r,n_r_-1,…,n₂,n₁分别小于x_r,x_r_-1,…,x₂,x₁,等号左边为多进制表达方式,括号外的下角标x_r,x_r_-1,…,x₂,x₁分别代表每一位的基底;等号右边为其对应的十进制数表示。同理,其倒叙后的新位序N’表示为:

(7)

\begin{array}{l} N' = (n_{1} n_{2} \dots n_{r - 1} n_{r})_{x_{1} x_{2} \dots x_{r - 1} x_{r}} \\ = n_{1} \cdot (x_{2} \cdot \dots \cdot x_{r - 1} \cdot x_{r}) + \\ n_{2} \cdot (x_{3} \cdot \dots \cdot x_{r - 1} \cdot x_{r}) + \dots + n_{r - 1} \cdot x_{r} + n_{r}, \end{array}

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

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 12点混合基FFT流程

Fig.2 12-point flow chart for Mixed-radix

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 混合基FFT程序框

Fig.3 Block diagram for Mixed-radix

2.3 比较计算

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

3 室内测试

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

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 标定计算过程

Fig.4 Progress of calibration

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

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 标定计算结果

Fig.5 Results of calibration

4 海上实测结果

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

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 实测海底大地电磁测深曲线

Fig.6 Measured seafloor magnetotelluric sounding curve

5 结论

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

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Steven

C C

, Arnold

S O

, Michael

, et al.

Marine magnetotellurics for petroleum exploration Part I: A sea-floor equipment system

[J]. Geophysics, 1998, 63(3),816-825.

DOI:10.1190/1.1444393 URL [本文引用: 1]

[2]

何继善

海洋电磁法现状和进展

[J]. 地球物理学进展, 1999, 14(1):7-39.

[本文引用: 2]

J S

The situation and progressof marine electromagnetic method research

[J]. Progress in Geophysics, 1999, 14(1):7-39.

[本文引用: 2]

[3]

赵宁, 王绪本, 秦策, 等.

三维频率域可控源电磁反演研究

[J]. 地球物理学报, 2016, 59(1): 330-341.