目前针对海底地形进行常规测绘的主要手段包括电子海图、单波束回声测深系统、多波束回声测深系统、侧扫声呐等。电子海图中包含等深线、海水深度等基本二维海底地形信息, 但不足以直观显示海底地形, 在对海底地形进行分析时难度较大^[1]。单波束回声测深系统较多波束系统而言, 设备体积小、操作简单, 但获取的数据样本稀疏, 无法对研究对象进行精细刻画, 黄谟涛等^[2]将单波束数据作为约束, 进而提高多波束测深数据的测量精度。多波束回声测深系统具有高效、精度高、分辨率高的特点, 在高性能计算机与多样化传感器等硬件技术的加持下, 多波束测深技术发展迅速^[3], 但存在单次测量获取的数据体量大、数据采集过程中易受外部环境干扰等情况。张目宁等^[4]提出将多波束水体数据根据入射角范围进行分区, 利用适当函数对分区内的反向散射强度波形进行拟合, 提高了拟合精度, 有利于后续进行对水体波形的提取工作, 丰富了多波束数据的应用方向。侧扫声呐系统通过记录海底反向散射回波强度, 生成灰度不一的声呐图像, 具有分辨率高、施工成本低、可以获得连续海底图像的优点, 但不能给出海底确切深度^[5], 赵建虎等^[6]针对多波束与侧扫声呐图像融合进行研究, 使得融合后的图像在分辨率和空间位置精度上优于原始图像, 但此方法受数据采集时船体运动姿态、海况等因素影响较大, 具有一定局限性。

科研鱼探仪过去主要应用于渔业领域, 是渔业领域的主要统计手段^[7], 在渔业资源调查与评估中具有重要作用。科研鱼探仪发射超声波信号, 遇到鱼体发生反射并由接收换能器接收鱼群回波信号^[8], 通过对接收数据进行滤波等快速处理, 可以大致判断鱼群的位置、规模和数量。由于渔业捕捞的即时性, 科研鱼探仪在使用过程中通常对采集到的数据进行实时处理分析, 并将处理结果在显示屏上实时显示, 指导捕捞作业, 数据具有体量大、来源丰富等特点。超声波回波信号除包含鱼群信息外, 也携带海底地形与海水信息, 因此利用科研鱼探仪数据进行地球物理分析成为可能。万树杰等^[9]对科研鱼探仪原始数据进行位置、高程校正等预处理工作后, 利用Surfer软件将获得的测深数据进行插值并网格化, 最终实现渔场海域海底地形的三维可视化。陈勃等^[10]利用科研鱼探仪电路信号监测数据提取脉冲幅值与发射间隔, 并对回波图像进行重构, 随后对重构后图像进行增强与去噪, 利用自适应跟踪手段提取海底地形信息, 实现了对海底地形的反演。目前在科研鱼探仪数据后续处理中, 主要使用Echoview软件进行数据提取。该软件为水声数据处理软件, 可以实现声呐设备数据的读取, 并实时显示回波图像。该软件读取的回波数据虽然包含水体深度与GPS信息, 但并不支持输出任意参数, 使研究人员无法对数据进行噪声压制等处理, 进而限制了数据在海洋地球物理等其他研究中的应用。

当采集与处理软件作为一个整体时, 难以获得原始数据, 或难以提取特定参数。在这种情况下, 需要对数据进行解编处理, 即通过使用头文件和其他相关信息来还原数据的格式。然而, 由于采集数据时所使用的执行标准与通用数据格式标准不一致, 这可能导致解编转储工作变得困难, 这也是地震探测早期经常遇到的问题。而为解决该问题, 彭波^[11]遵循SEG协会标准, 针对哈萨克斯坦二维非标准原始地震数据的解编转储工作进行研究, 使得非标准原始地震数据可以适用于Promax等主流地震处理软件; 洪天聪等^[12]针对DFS-V地震仪磁带记录SEGB格式进行解编, 并升级为微机记录通用SEGY格式, 便于数据备份与传输。科研鱼探仪数据在过去通常用于实时观测指导捕捞活动, 虽然Echoview软件可以读取科研鱼探仪数据, 但并不支持输出海洋地球物理研究所关注的船速、船向、采样间隔、采样点个数等特定参数, 无法进行噪声压制等处理工作。为了将科研鱼探仪数据应用至海洋地球物理研究中, 拓宽研究海底地形地貌与海水温盐结构的数据来源渠道, 针对科研鱼探仪数据的解编处理工作亟待进行。

本文将基于MATLAB程序对科研鱼探仪采集数据进行解编处理, 同时对包含不同信息的数据文件重新分类, 有针对性地提取与海底地形地貌、海水温盐结构等海洋地球物理研究相关的数据信息, 如水体声速、采样数、GPS信息等, 并进行研究分析。以期将时测时用的科研鱼探仪数据应用到海洋地球物理研究中, 降低数据采集成本、避免资源浪费, 为后续噪声压制等处理工作提供数据来源。

1 研究方法

图 1所示为科研鱼探仪数据采集原理示意图。科研鱼探仪接收传感器应固定在刚性基础上, 尽可能远离机舱、螺旋桨以减少采集过程中的干扰噪音。在进行采集工作时, 科研鱼探仪收发器发射超声波信号, 信号传播过程中遇到鱼群、海水结构或海底等声阻抗界面时, 产生回波信号并由接收器接收, 实现寻找鱼群的目的, 这也使得海洋地球物理工作者可以从科研鱼探仪数据中获取海底地形特征、海水结构等信息。本文选用的EK80科研鱼探仪工作频率分别为18、38、70、120、200和333 kHz, 工作频率降低时, 观察范围变大, 与此同时分辨率降低, 而探测深度随功率增加而增加。假设海水盐度为3.5%、温度为+10 ℃、探测条件为表面后向散射强度–20 dB的平均底部, 噪音为典型船舶噪声水平, 在这种情况下, 工作频率38 kHz时, EK80科研鱼探仪的探测深度在2 100~3 000 m; 工作频率200 kHz时, 探测深度则在450 m左右。EK80原始数据文件中的数据报序列及部分参数值取决于采集过程中的频道设置。在本文使用的原始数据中, 工作频率为18 kHz时采样间隔(SI)为2.8×10^–6 s, 38和70 kHz时采样间隔为4.8×10^–6 s, 120 kHz时采样间隔为4.0× 10^–6 s, 200 kHz时采样间隔为3.2×10^–6 s, 333 kHz时采样间隔为2.4×10^–6 s。

解编工作前期主要借助SIMRAD EK80科研鱼探仪参考手册, 研究二进制数据文件的封装方案, 对科研鱼探仪原始数据的数据格式、数据结构和存储类型进行分析, 并明确手动设置参数对EK80原始数据存储的影响。在此基础上依托Matlab平台进行编程作业, 选用恰当的数据类型对原始数据中的参数进行解编转储, 最终实现科研鱼探仪原始数据解编、特定参数的提取与保存。

1.1 数据结构

EK80数据文件根据数据采集时的频道设置, 以表 1所示的基本文件格式, 存储科研鱼探仪工作时相关数据信息, 具体分为数据类、环境类、船体类。数据类包含Configuration数据报、Filter数据报、Parameter数据报、Sample数据报和Annotation数据报。Configuration数据报为原始数据文件中首个数据报, 记录收发器、通道和接收传感器的相关信息, 相关参数在原始数据文件中不再发生更改。收发器与接收传感器之间的对应关系称之为通道, 通道名称用于标识不同数据报中同一频率下采集的数据。每一通道包含两个Filter数据报, 首个数据报存储来自收发器的滤波参数, 其次存储来自处理单元软件的滤波参数。样本数据由Parameter和Sample数据报组成, Parameter数据报包含有数据传输过程中每道所需的相关参数; Sample数据报存储样本数据, 数据类型取决于DataType参数。Annotation数据报存储用户备注信息。环境类包含Environment数据报, 存储环境参数如水深、酸度、盐度、声速、温度等, 由用户输入用于获取更为准确的深度和回波信息。船体类包含NMEA数据报、MRU数据报。NMEA数据报存储GPS导航设备相关信息, MRU数据报存储特定时间内运动传感器信息。

表 1中最后4项数据报与样本数据异步时, 可能会交错存储, 给解编工作造成了困难。

图 2所示标准封装方案应用于所有二进制数据文件, 二进制数据报开头与末尾均设置相同类型的长度标记, 用于判断计算机读取数据文件的顺序与完整性。所有二进制数据报使用相同的报头, DatagramType参数前3个字节标识数据类型, 最后1个字节标识数据报版本。DateTime表示自1601年1月1日以来100 ns间隔数量, 是Windows NT操作系统使用的内部文件时间。XML数据报较二进制数据报而言更加灵活, 且不会被压缩。XML数据报类型由首个标签名称决定, 长度由封装时的数据参数确定。

以典型数据信息NMEA(National Marine Electronics Association, NMEA)为例, NMEA-0183为GPS导航设备统一RTCM标准协议, 原始数据中NMEA结构体定义如图 3a所示, 其中包含数据报报头及GPS相关信息。＄GPGGA、＄GPGST、＄GPHDT、＄GPRMC、＄GPVTG、＄GPZDA均为NMEA-0183协议定义的语句, 对应为定位信息、伪距噪声统计、航向信息、推荐定位信息、矢量跟踪与对地数据、时间信息。

解编字段内容如图 3b所示, 以＄GPGGA为例: 原始数据中字段之间以逗号间隔, ＄GPGGA表明为GGA格式(Global Positioning System Fix Data)的GPS定位信息; 134055.89为UTC时间, 时分秒格式hhmmss.ss; 2127.460836为纬度, 度分格式ddmm.mmmm; N代表北纬; 11558.243914为经度, 度分格式dddmm.mmmm; E代表东经; “2”为GPS状态, 代表差分定位(DGPS), 另外有0为不可用(FIX NOT VALID)、1为单点定位(GPS FIX)、3为无效PPS、4为实时差分定位固定解(RTK FIX)、5为实时差分定位浮点解(RTK FLOAT)、6为正在估算, 共计7种GPS状态; “12”为正在使用的卫星数量, 范围为00~12; “0.7”为HDOP水平精度因子, 范围为0.5~ 99.9; “1.89”为海拔高度, 范围为–9999.9~99999.9, M代表海拔高度的单位为米; “5.89”为地球椭球面相对大地水准面的高度WGS84水准面划分, M代表WGS84水准面划分的单位为米; “33.4”为差分时间, 即从接收到差分信号开始的秒数; “0023”为差分站ID号, 取值范围在0000~1023。当GPS状态不为差分定位时, 与差分有关的参数将定义为空。其中GPS时间、经纬度、海拔高度等信息在后续的资料处理中可以用于标明特殊现象所在位置, 进而通过换算得到对应现象的尺度。

1.2 数据提取

根据采集时的系统设定, 在科研鱼探仪原始数据文件中包含有仪器版本、IP地址、脉冲时长、采集频率、环境参数、运动姿态、GPS等众多信息。海洋地球物理的相关研究主要关心与海底地形地貌、海水温盐结构、数据采集位置等有关的数据信息。所以在数据解编过程中需要对数据进行筛选, 提取出与所研究问题密切相关的参数数据。在科研鱼探仪数据采集过程中, 信号通过收发器以特定频率激发, 接收传感器以一定采样间隔接收超声波回声信号。本文选用的科研鱼探仪原始数据包含6组激发频率, 对应6通道, 工作频率分别为18、38、70、120、200和333 kHz。

基于Matlab软件对原始数据文件进行解编。首先定义方法, 输入为原始数据文件, 输出为原始解编文件与经提取后的解编文件。根据数据结构中定义的参数数据类型, 在Matlab中计算相应参数占用的字节数并利用fread函数进行数据读取。XML类型数据报和字符串参数则利用fcanf函数进行读取。原始数据解编思路如图 4a所示, 依次从原始数据中读取XML类型的Configuration和Parameter数据报。因每一通道包含分别存储收发器和处理单元软件滤波参数的两个Filter数据报, 故设置循环参数为12, 依序读取Filter数据报。SDVLW为“NME0”类型结构体, 存储GPS导航设备相关信息。依次读取存储有温度、水深、盐度等环境信息的Environment数据报和存储有传感器类型与编码方式等信息的MRU数据报。每个通道数据均含有1个NMEA数据报, 存储当前的位置信息, 故设置循环参数为6。最后从Motion数据报中提取数据采集时的运动姿态信息。至此将原始数据文件中与采集环境、位置、船体姿态等相关的信息解编至Datagram结构体, 原始解编数据如图 4b所示, 随后可以提取采样间隔、样点能量值等感兴趣的参数。

对Parameter数据报中所需参数进行提取的流程如图 5所示。原始数据中Parameter数据报为XML类型数据报, 经解编后以char类型字符向量存储在InitialParameter变量中。具体流程为: 首先利用strfind函数对InitialParameter变量中标签值"ES"和结束符号“/> ”进行检索, 返回值为对应字符向量的索引值, 分别存储在idxStart和idxEnd数组中。idxStart数组长度应与原始数据采集时定义的通道数相同。所以首先需要判断idxStart是否不为空, 当idxStart不为空时, 以idxStart数组长度作为循环次数, 以通道名称为判断标准, 对Datagram中InitialParameter变量进行拆分提取, 分别将6个通道的通道频率、采样间隔与海水声速信息提取至对应解编文件。

对样本数据进行提取的流程如图 6所示。首先读取原始数据文件中的PingSequence数据报。以接收传感器对应道数为依据进行数据归类, 将每道对应的6组Ping ChannelID存储在Datagram结构体的PingSequence元胞中。之后进入循环提取每道样点数据, 判断文件类型并读取每道对应的Parameter数据报, 存储在Datagram的Parameter元胞中。判断并读取“RAW3”类型数据报, 提取通道名称判断并赋值给SeqNo变量, 之后根据道数和通道名称将读取的数据类型等信息以结构体的形式存储至Sample元胞对应位置。每道样本个数取决于数据类型与记录点数count, 本文选用的数据中每个记录点包含4个来自四分量接收传感器的复数, 故计算每道应读取多少字节时应考虑实部与虚部, 采用4*2*count作为读取字节数。定义元胞Ping用于存储经提取的所需参数。通过循环赋值可以确定各通道的数据道数并存储至元胞Ping相应位置, 将循环中计算的该道样本个数存储至元胞Ping中指定数组。根据通道名称的不同, 分别判断并读取每道数据采集时对应的GPS信息与船体姿态信息。进入下次循环的标志为数据类型为“MRU0”的数据报读取完毕, 并且在下次循环开始时首先判断文件读取是否结束, 如果文件读取完成则循环结束, 若文件未读取完成, 则继续读取文件信息。

图 7a所示创建样本能量值矩阵的流程图中, 循环次数为原始数据通道数, 首先判断对应通道的样本个数是否为零, 若样本个数为零, 则将对应Ping元胞中的SampleNumber与SampleData变量设为空, 否则将样本个数数组中的最大值赋给SampleNumber, 并创建矩阵维度为样本个数*道数的空矩阵SampleData。

图 7b所示对样本能量值及位置信息进行提取的流程图中, 分别判断各个通道SampleNumber值是否为零, 当SampleNumber不为零时, 判断Ping元胞样本个数数组中的各个值是否均与SampleNumber相等, 若均相等, 则对每道数据进行赋值, 数值为对应位置复数的模量。同时利用NMEA中存贮的GPS信息读取各道采集时所在的经纬度。若不相等, 则根据各道对应样本数进行样本数据赋值与位置信息读取。

图 7c所示流程图含义为将各个通道SampleData矩阵中的样点能量值数据转换为分贝值。至此完成科研鱼探仪原始数据解编与特征参数解编提取工作, SampleData矩阵中存储对应采集位置处的样点能量值信息。

2 地球物理应用

经过上述解编与提取, 不同通道所包含的位置信息、样点能量值等数据保存至Ping元胞中, 如图 8所示。其中通道名称ChannelID为18 kHz, 即采集时的工作频率; ChannelMode为0, 代表该频道的接收模式; SI为数据采样间隔; V为水体声速; ntr为该通道的数据道数; ns存储每道对应的样本点数; SampleNumber为各道中最大样本数; SampleData为样点能量值矩阵; Position记录经纬度信息。

本文选用的科研鱼探仪数据由科考船“东方红3”上搭载的Simard EK80型科研鱼探仪采集得到, 6个换能器的工作频率分别为18、38、70、120、200和333 kHz。调查时间为2022年5月11日—2022年5月18日, 数据采集期间船速范围为4~6 kn, 海况良好, 采集时多波束回声剖面仪、浅地层剖面仪同时工作, 可能会对科研鱼探仪数据产生干扰, 具体见后续分析。调查地点为F站位冷泉, 采集位置如图 9所示。

海底冷泉是指来自海底沉积层或以下的低温气体或流体, 以喷涌或渗漏的方式逸出海底, 与周围海水相互作用, 形成的一种海洋地质现象^[13-14]。目前可使用海底可视化技术、地震探测等方法进行海底冷泉观测, 其中海底可视化技术能够对已知位置的冷泉系统进行直观实时的观测, 但不适用于大规模普查; 地震方法中浅地层剖面系统能够捕捉到冷泉羽状流产生的特征响应, 但精度较低^[15]。科研鱼探仪发射的超声波信号遇到海底、冷泉羽状流、鱼群等声阻抗界面时会产生回声, 可以作为常规冷泉探测方法的补充, 丰富冷泉探测数据来源。

对科研鱼探仪解编数据中样点能量值矩阵SampleData利用Matlab绘图, 结果如图 10所示, 根据采样点数、采样间隔、水体声速可以计算出水体深度信息, 区域水深在1 000~1 500 m。在横向上, 由于船体具有一定速度, 水体结构如海洋内波、冷泉羽状流等也具有流动速度, 二者速度之间可能存在夹角, 后续在横向距离换算上需要考虑多普勒效应, 结合船速、水速、方位角等因素进行计算。

图 11中可以看到形似羽流的声学响应, 结合数据采集区域为自生碳酸盐岩区, 存在冷泉气体逸出现象, 满足冷泉存在的条件, 故认为该处羽状扰动为冷泉羽状流在科研鱼探仪采集数据中造成的声学响应。在图 10、图 11中均可以看到背景存在线性干扰, 结合数据采集时的情况, 分析可能为多波束回声剖面仪、浅地层剖面仪等仪器造成的信号干扰, 或海洋涌浪造成的环境噪声干扰, 后续将尝试通过线性干扰压制等方式进一步处理分析。

注意到鱼探仪所记录的温度、盐度、密度等信息均为浅层信息, 无法直接拓展到整个剖面。但该部分信息可以用于我们后续反演的约束, 即基于鱼探仪提供的振幅剖面, 计算其反射系数, 进而实现对水体的温盐密的反演。此时, 浅层测量的数据可以作为反演的约束, 保证反演的准确性^[16]。上述内容为正在进行的研究。

3 结论

本文针对SIMRAD EK80科研鱼探仪原始数据进行分析, 在详细论述原始数据存储格式的基础上, 对解编思路进行阐述。在解编过程中需要明确原始数据中对二进制数据采取的封装方案, 按照其中规定数据类型占用的字节数进行读取。对于占用存储空间不固定的参数需要根据参考手册中的说明, 计算相应参数占用的字节数。针对XML类型与字符串类型的参数则使用fscanf函数进行读取。最终实现原始数据解编和特征参数读取。通过对过去时测时用的科研鱼探仪数据进行解编转储, 并对解编数据中的样点数据进行绘图, 同时结合采样点数、水体声速等信息计算深度, 可以得到垂向上准确的高分辨率海水及海底剖面。但由于船速、水速等因素导致横向上存在多普勒效应, 需要在后续研究中结合方位角、速度等信息进行距离校正。对于图 10、图 11中存在的线性干扰, 推测为测量过程中多波束回声剖面仪等仪器造成的信号干扰, 针对科研鱼探仪数据进行去噪处理也将是后续的研究方向。针对科研鱼探仪数据的解编转储与特定参数提取工作, 可以充分利用数据资源, 节约数据采集成本。而且通过对解编数据进行噪声压制和处理, 可以得到高分辨率的海水、海底剖面, 为海洋地球物理相关研究提供数据支持, 丰富海洋地球物理研究资料的数据来源。