《1 前言》
1 前言
深海钴结壳富含钴、铂等稀贵金属,是一种重要的战略矿产资源。它以薄薄的一层覆盖在水深 500 ~3 000 m 的平顶海山、海台顶部或坡上,厚度一般为 2 ~6 cm,并对其下层岩石有很强的粘附性,很难轻易被剥离[1 ~3]。因此在开采过程中,获取地表高程信息来计算最佳采集深度是一项关键技术。根据 John E.Halkyard 的钴结壳开采方案,探测装置被安装在海底采矿车的底部,并具备探测精度高、实时性强的特点[4]。笔者结合钴结壳的特殊赋存环境和开采方法,对比研究了机械式微地形探测仪、单波束回声测深仪、多波束条带测深仪、侧扫声纳、测深侧扫声纳、水下激光测深系统、激光微地貌扫描仪、水下立体摄影和水下电磁波探测系统等一系列基于声、光、电的水下地形探测方法的优缺点,首次提出了摆动式超声单波束探测方法[5]。为了验证该方法的可行性,主要从探测装置研制和数字高程数据处理两方面展开实验研究。
《2 摆动式超声单波束探测实验台》
2 摆动式超声单波束探测实验台
《2.1 工作原理》
2.1 工作原理
由图 1 所示,收发一体水声换能器发射一定频率的超声波照射到水底面上,其横向波束角为 β(纵向等同),当碰到底面时产生反射或散射回波。回波信号接收电路在超声波发射同时就进入工作状态,一并接收发射信号和回波信号。当每次接收完一个探测点的回波信号后,由步进电机驱动水声换能器绕摆动中心转动一定角度进行下一次探测。由此半个摆动周期探测一条横向测线。在水声换能器工作的同时,行走机构一直沿着导轨前进。数据采集与处理软件将采集到的发射信号作为参考信号,与回波信号进行振幅相关性分析计算出超声波的渡越时间以及探测距离,再以摆动中心为坐标原点,并结合摆角和行走速度两参数计算出每个测点的三维坐标值和高程值。
《图1》
图1 摆动式超声单波束探测实验台示意图
Fig.1 Sketch map of single -beam swing ultrasonic laboratory bench
《2.2 结构组成》
2.2 结构组成
实验台按物理结构分类,可分为超声探测装置、机械传动装置、数据采集装置、行走机构和微地形模拟料五大部分,如图 2 所示。
《图2》
图2 摆动式超声单波束探测实验台实物图
Fig.2 Physical map of single -beam swing ultrasonic laboratory bench
机械传动装置采用圆柱凸轮机构设计,其优点是结构紧凑,占据空间小,无急回特性,安装便利。笔者采用滑块行程为 80 mm,摆杆摆角为 2°,要求摆杆能够每摆 0.5° 停顿一次进行设计。结果计算出摆杆的摆动角度最大偏差只有 0.01°,可认为摆杆做等角度(0.5°)摆动。
超声探测装置主要包括水声换能器和水下超声测量仪两大部分。水声换能器采用压电陶瓷材料,并按工作频率为 500 kHz,波束角为 3° 进行设计。水下超声测量仪包括超声波发射系统、回波接收系统两大部分。超声波发射系统主要包括脉冲信号产生电路、功率放大电路、换能器阻抗匹配网络;回波接收系统主要包括前置放大电路、带通滤波电路、时间增益控制电路和数据通信接口。其工作原理如图3 所示。
《图3》
图3 水下超声测量仪工作原理图
Fig.3 Principle diagram of underwater ultrasonic measuring instrument
数据采集装置采用 PCI -1714 高速采集卡,在工控机工作的数据采集与处理软件是由 VC ++ 与 MATLAB 语言混合编程而成,并以动态链接库的形式将数据采集和数据处理程序集成到一起。它主要包括数据通信和存储、数字滤波以及振幅相关性运算等功能。行走机构采用步进电机驱动,并由 AT89C52 单片机控制,实现行走速度可调。构造水底微地形的材料有普通河沙、模拟海泥和硬质水泥三种。
《2.3 垂直探测精度检验》
2.3 垂直探测精度检验
当实验台建成后,利用普通河沙、模拟海泥和硬质水泥三种材质的微地形进行垂直探测精度检验。探测高度从 400 mm 到 1 400 mm,每隔 100 mm 探测一次,每一种底质测量了 11 个不同高度值,测量结果由表 1 所示。
《表1》
表1 三种材质的精度测量结果比较
Table 1 Three kinds of materials compared the precision of measurement results
mm
将三种材质的测量误差合在一起进行统计分析,样本大致服从正态分布。利用极大似然估计方法对该正态分布的未知参数进行估计,得出其置信度在 95 % 的各参数值:均值 μ 为 1.080 3,均方差 σ 为 4.272 5, μ的置信区间为 [ -0.434 7,2.595 3 ],σ 的置信区间为 [ 3.435 9,5.651 2 ]。根据 3 σ 规则,对于正态随机变量来说,样本以 95.45 % 的概率落入区间 [ μ - 2σ,μ + 2σ ] ,即区间 [ -7.464 7 mm,9.625 3 mm ] 内。由此发现该区间的测量误差值全部处在区间 [ -10 mm,10 mm ] 范围内。结果表明实验台的垂直探测精度可达到厘米级。
《3 数字高程数据处理技术》
3 数字高程数据处理技术
《3.1 基于改进型数据加窗法的高程异常值剔除方法》
3.1 基于改进型数据加窗法的高程异常值剔除方法
超声信号因掠射角、底质特性、外界噪声等干扰因素的影响,其反射回波会远离水声换能器的接收位置,只能接收到探测点的散射回波。当探测点的散射回波隐藏在干扰因素造成的回波中时,利用振幅相关法检测到的最大相关值有可能不是探测点真实回波产生的,这样不可避免地导致高程异常值出现。吴英姿博士采用数据加窗法处理 H/HCS -017 型条带测深仪的水深数据时,依据地形变化的坡度小于 45° 的经验准则[6]。摆动式单波束探测与多波束探测的工作原理不同,但它通过水声换能器的摆动使之达到横向大扇区、全覆盖的窄带探测,并且水声换能器自身的高指向性等价于多波束的预成波束,因此该经验准则也适合于摆动式单波束探测。笔者依据该准则,并根据摆动式单波束探测原理和振幅相关检测法工作特点,提出了一种改进型数据加窗法。其具体算法如下:
3.1.1 确定高程参考值 z′0
1) 由图 1 所示,设摆动中心为坐标原点,则 都为负值。 i 表示 X 方向(探头摆动方向)的测点序列, j 表示 Y 方向(行走方向)的测点序列。
2)计算 Z 值的均值 μz 与均方差 σz :
理论上应采用 Min() 作为 z′0 ,但 Min()有可能是异常值。根据 3 σ 准则,当 Min() ≤ μz -2σz ,则 z′0 取 Min() ;否则, z′0 取 μz -2σz 。
3.1.2 确定 X 方向数据窗尺寸
设同一测线上两相邻测点在 Z 方向的坐标值为 ,其对应的摆角为 θi , θi +1 ,依据地形变化的坡度小于 45° 的经验准则, 应满足:
式(3)中,m 为冗余量系数,则
令 ,则设定 X 方向数据窗宽为 X CH,窗高为 2X CH。
3.1.3 确定 Y 方向数据窗尺寸
假设行进速度为 v ,水声换能器摆动周期为 T ,沿 Y 方向的相邻探测点的坐标值 应满足:
式(5)中,m 为冗余量系数。 令,则设定 Y 方向数据窗宽为 Y CH。
3.1.4 确定回波能量阈值
由于导致高程异常值产生以掠射角、底质特性和外界噪声影响为主,一条测线所对应的地形范围极小,其底质特性可认为是均质;而掠射角受微地形随机起伏影响,无规律可循,即使水声换能器摆角为 0° 时探测的高程值也会出现异常,由此在数据加窗判别高程异常值中恒定选取正常的高程值作为起始加窗点是一个难题。
一般来说,如果水声换能器接收到探测点的回波是反射回波,则认为其高程值是正常的;而散射回波容易受外界因素干扰,需要借助其他方法判断其高程值的正确性。故利用回波能量这一特点参与到数据加窗中实现起始加窗点的选取工作。但是反射回波与散射回波没有明确的界限。为提高数据加窗判别的准确性,计算每条测线上的回波能量的平均值 μW :
式(6)中,n 为测线上的测点数;W 为测点的回波能量值;将 μW 设为阈值,认为回波能量低于该阈值时,其对应测点的高程值可能会出现异常。
3.1.5 高程异常值判别步骤
1)将坐标值 排序,以接近均值 μz 为原则,由此将最接近均值 μz 的坐标值作为起始加窗点。
2) 引入回波能量,若起始加窗点的回波能量不低于平均值 μW ,则以该点为中心向其测线两边的坐标值依次加窗判别;反之,则将仅次于它接近均值μz 的坐标值作为起始加窗点,重新引入回波能量进行判别,依次类推。
3) 以加窗点 为窗口中心,当 与 都不在窗口内,则认为 为异常值,并剔除;当 或 在窗口内,则借助 Y 方向数据窗口进行判断,若不满足式(5)则认为 为异常值,并剔除。其他情况可认为 是正常值。
《3.2 基于支持向量机的高程值预测方法》
3.2 基于支持向量机的高程值预测方法
当高程异常值被剔除后,会留下测点空缺;此外,摆动式单波束探测存在一个不可避免的问题,就是背向超声波的地形坡面无法被探测,只能通过一些理论方法去预测。支持向量机是在小样本情况下发展起来并建立在 VC 维和结构风险最小原理基础上的统计机器学习理论,具有严格的数学基础[7]。但在运用支持向量回归算法预测时需要解决以下两个主要问题。
3.2.1 线性回归与非线性回归的选择
支持向量回归分为线性与非线性两种类型。当被测微地形是平地面或斜坡面时,支持向量线性回归的预测精度高于非线性回归,尤其是在高程值存在随机误差或高程异常值未被完全剔除情况下。笔者对此采用线性相关系数法,该方法计算简单,运算量小。
设给定一条测线的高程样本集:
线性相关系数 r 计算表达式如下:
对于离散的高程数据样本,线性相关系数 r 的计算表达式如下:
式(8)中
根据经验,当>0.8 ,支持向量机采用线性回归比较合适。
3.2.2 核函数及其参数的确定
事实上,水下微地形较陆地免受风化、雨水侵蚀以及人为活动等诸多因素影响,并加上沉积作用,凸峰地形具有较为明显的对称性。高斯型函数能很好地描述它,其基本表达式如下:
式(9)中, A , x0 和 σ 分别代表峰高、峰位和峰宽[8]。故采用高斯径向基作为支持向量非线性回归的核函数,其参数 σ 是确定对象。因参数 σ 可直接反映凸峰地形峰宽属性。笔者采用 2 σ 作为凸峰地形的坡长,这一点可由正态分布高斯函数证明:
高斯核函数中的参数 σ 与正态分布高斯函数中均方差参数 σ0 存在如下关系:
当正态分布服从(0,1)分布时,即 σ0 =1 ,则 2σ≈ 2.828 ,查表可得 2 σ 的概率约为 99.767 %。由图 4 所示(图 4 中 p1 表示 σ),当 2σ =4 ,函数 f (x) =4· exp( -x2/4) 的坡长近似为 4,并且在其他参数不变情况下, σ= 的回归曲线更接近函数 f (x) 。 由此知道了测线上的坡度、坡长等信息就可以确定参数 σ 的值。另外,误差惩罚参数 C′ 是特征子空间中调节经验风险与置信范围之间的比例[9]。通过它可以控制训练精度,在训练错误和推广能力之间取得折中。 C′ 越大表示对超出回归允许最大误差 ε 的样本的惩罚越大,而 ε 由探测精度决定,因此参数 C′ 和 ε 的选取需要先验知识确定。
《图4》
图4 基于高斯核函数的非线性回归示意图
Fig.4 Nonlinear regression diagram based onGaussian kernel function
《4 实例》
4 实例
图 5 是一块长 900 mm,宽 700 mm 的真实微地形,经过探测实验后得到图 6 所示的网格地形图。
《图5》
图5 真实的随机微地形
Fig.5 True random micro -terrain
可以看出图 6 中存在 4 个明显的高程异常值。先用改进型数据加窗算法(冗余量系数 m =1.2)剔除异常值后,再由支持向量机回归预测(采用 ε -insensitive 损失函数,核函数选用高斯径向基函数,参数 C′, ε 和 σ 分别取 1 000,2 和 30,线性回归不考虑参数 σ),在每条测线上得到间隔为 25 mm 的高程值,由图 7 所示。鉴于网格地形图可读性差,利用 AVS/Express 软件将图 7 渲染后,如图 8 所示。通过图 5 与图 8 外观比较可知,数字微地形能如实反映真实微地形的起伏特征。
《图6》
图6 原始的数字微地形
Fig.6 Original digital micro -terrain
《图7》
图7 高程数据处理后的数字微地形
Fig.7 Digital micro -terrain of elevation data post -processing
《图8》
图8 渲染后的数字微地形
Fig.8 Digital micro -terrain of post -rendering
《5 结语》
5 结语
以确定钴结壳微地形探测方法、实施数字高程数据处理为研究内容,取得以下成果:
1)针对摆动式超声单波束探测方法搭建实验台,其垂直探测精度可达到厘米级。
2)根据探测装置的工作原理,提出了一种改进型数据加窗算法,进一步提高了高程异常值被剔除的准确性。其改进的主要措施有:a. 将回波能量作为探测点的内在属性参与确定起始加窗点;b. 加窗点采用实际探测点,并要求起始加窗点尽可能靠近一条测线上实际探测点的高程平均值。
3)在基于支持向量回归算法的高程值预测中:a. 提出采用线性相关系数法解决线性与非线性回归的选择问题;b. 确定高斯径向基作为核函数,提出并证明了参数 σ 可根据测线上的坡度、坡长信息来确定。