简介《多波束水深测量数据处理》面向海洋测绘、水深测量及多波束系统应用人员聚焦面测深数据从原始观测到成果归算中的坐标系统与误差改正问题。内容围绕测量船、多波束探头、姿态传感器、电罗经、水平及测量坐标系展开讨论横摇与纵摇改正、声线处理中的坐标转换、航向归算以及系统偏差和测量误差的分析与处理思路适合需要理解多波束数据处理流程、排查精度问题的技术读者参考。压缩包内含1个doc文档整体约152KB篇幅精炼便于按主题查阅坐标系定义、改正公式与误差影响分析。已有1125人学习下载可为多波束测深内业处理、成果检核和作业方案复盘提供基础参考。1. 多波束面测深的隐形成本一 Ping 上千个波束每个都要单独归算很多人第一次接触多波束注意力都在换能器上——扇面开角多大、波束数多少、最大测深多深。真正决定成果能不能用的是换能器之后那条看不见的归算链路。一个 Ping 打出去上千个波束带着旅行时间和波束单元入射角回来每一个都要先按声速剖面折射成水深和侧向偏距再经横摇、纵摇、航向三次旋转落到水平坐标系最后叠上吃水、潮位、涌浪才成为深度基准面下的高程。这条链路上没有哪一步会主动报错。坐标系定义反了、符号取错、旋转顺序颠倒数据照样出图只是地形安静地偏了几十厘米。面测深提高了效率也把系统性错误的隐蔽性一起放大了。2. 六套坐标系与姿态改正船体一歪水深为什么要重算2.1 一条船为什么要定义六个坐标系多波束系统不是一台仪器是一组传感器捆在一条会摇的船上。除了多波束测深仪和定位设备还有测航向的电罗经、测纵摇横摇的姿态传感器、测船体上下起伏的涌浪滤波器。每个设备都有自己的安装位置、安装姿态和零点输出量的物理含义也各不相同。要让它们协同先得给它们一套共同的语言。测量船坐标系 XbYbZb 是这套语言的根右手系站在船艏方向看向左为 X 轴正向沿船艏艉线为 Y 轴正向垂直向下为 Z 轴正向原点一般选在多波束探头的发射中心。把原点放在探头上而不是船重心是因为所有波束的旅行时间都从探头起算归算起点和测量起点重合能省掉一次平移。围绕这个根派生出五个坐标系。探头坐标系 XcYcZc 与测深扇面垂直Yc 指向船艏Xc 向左Zc 沿中心波束方向原点同在探头发射中心姿态传感器坐标系 XdYdZd 和电罗经坐标系 XeYeZe 的三轴与船体系对应轴同向平行用来描述“传感器装在哪、朝哪”水平坐标系 XaYaZa 与船静止平放时的船体系重合XaYa 平面平行于当地水平面横摇、纵摇、航向倾斜三项改正都以它为基准测量坐标系 XY 一般是高斯平面坐标系水深点经航向归算后最终落在这里。定义这么多套不是为了把公式写复杂。每个传感器只在自己那套坐标系里说话改正在两两之间的变换里做一次做完就进入下一层。只用一个坐标系硬套横摇改正、罗经倾斜改正、安装偏差改正会全部搅在一起符号一错就再难拆开。2.2 横摇纵摇一次旋转和三次旋转是等价的姿态传感器测出横摇 α 和纵摇 β 之后船体系到水平系的变换就是一个旋转。从物理上看船体系绕空间某条直线 l 旋转一个合成角 ρ就能和水平系重合。这是最直观的一次旋转表达ρ 由 α、β 按tan²ρ tan²α tan²β一类关系合成旋转轴的方向则由两个姿态角共同决定。问题在于声线处理时这套公式不好用。声线在探头坐标系里按波束单元入射角展开每个波束要的是自己相对水平面的角度而不是整条船绕某根合成轴转了多少。于是换一种写法先把水平系绕 Za 轴转一个偏角 ε让 Ya 轴与 Yb 轴在 XaYa 平面上的投影重合再分别绕 Y 和 Xa 各转一次得到与船体系重合的姿态。这样大角度被拆成三个小步骤每一步都有明确的几何含义。import numpy as np def rot_x(t): return np.array([[1, 0, 0], [0, np.cos(t), -np.sin(t)], [0, np.sin(t), np.cos(t)]]) def rot_y(t): return np.array([[ np.cos(t), 0, np.sin(t)], [0, 1, 0], [-np.sin(t), 0, np.cos(t)]]) def rot_z(t): return np.array([[np.cos(t), -np.sin(t), 0], [np.sin(t), np.cos(t), 0], [0, 0, 1]]) def hull_to_level(alpha, beta): 船体系 - 水平系分解形式。 alpha: 横摇(rad)beta: 纵摇(rad) 返回 3x3 矩阵 R满足 [xa,ya,za]^T R [xb,yb,zb]^T # 第一步偏转角 eps把 Yb 在水平面上的投影对齐到 Ya eps np.arctan2(np.tan(beta) * np.sin(alpha), np.tan(beta) * np.cos(alpha) 1e-12) # 第二步、第三步两个中间角闭式解见合成旋转公式 beta_p np.arctan(np.tan(beta) * np.cos(eps)) alpha_p np.arctan(np.cos(beta_p) * np.tan(alpha)) return rot_z(-eps) rot_y(-beta_p) rot_x(-alpha_p)这段代码里三个角不是随便取的ε 负责把船体纵轴的水平投影摆正β 是投影摆正之后的纵摇分量α 是残余的横摇分量。三个矩阵从右往左乘顺序对应“先绕 Xa 转、再绕 Y 转、最后绕 Za 转回”。顺序颠倒得到的姿态矩阵就不正交后面所有波束的方向都会带一个缓慢变化的偏斜图上表现为测线边缘的假地形。2.3 为什么声线处理偏爱分解形式因为声线归算需要的是每个波束相对水平面的入射角 θ而换能器给出的只是波束单元入射角 γ——以换能器垂直轴为基准的角度。船一摇探头跟着摇γ 就不再等于 θ。要把 γ 换成 θ必须知道探头相对水平面到底转了多少而分解形式恰好把“转了多少”拆成了 ε、β、α 三个可以单独引用的量。偏转角 ε 还有一个额外好处它可以推迟到航向归算时再统一处理。航向归算本来就要把真航向 A、高斯坐标收敛角 γ 和坐标偏转角放在一起算ε 混进去正好一趟做完不必每个波束都重算一遍反三角函数。这正是几千个波束一个 Ping 时最在意的开销。提示分解形式里的 α、β 与姿态传感器直接读出的 α、β 不是同一个量。代码里如果直接把 α、β 塞进旋转矩阵短测线看不出来长测线会在边缘波束上积累成米级偏差。3. 航向归算光纤罗经和陀螺罗经测的根本不是同一个物理量3.1 光纤罗经三个正交陀螺测的是角速度分量光纤陀螺的原理是 Sagnac 效应光束在环形光路中前进时如果环路本身在转动顺逆两个方向的光程相对静止时都会发生变化光程差正比于环路的转动角速度。三个互相垂直的光纤陀螺装在一起就能分别测出地球自转角速度在三个方向上的分量。船静止水平时地球自转角速度 ωe 可以分解为水平分量 ω1 ωe·cosφ 和垂直分量 ω2 ωe·sinφφ 是当地纬度。由三个陀螺读出的 ωx、ωy、ωz就能解出船舶真方位 A。但船在摇陀螺测到的是倾斜坐标系三轴上的分量不是水平系的分量。必须借助光纤罗经自带的姿态传感器测出的横摇、纵摇先把分量转回水平系再解算真方位。def fiber_heading(wx, wy, wz, roll, pitch): 倾斜系角速度分量 - 水平系 - 真方位(rad) v np.array([wx, wy, wz], dtypefloat) R hull_to_level(roll, pitch) # 船体系到水平系的旋转 v_level R.T v # 把分量转到水平系 return np.arctan2(v_level[0], v_level[1])这里有个容易被忽略的前提姿态传感器和光纤陀螺的坐标系必须一致对应轴平行。一体化设备出厂时通常已经保证分离设备就必须先做静态测定。否则R转的不是同一组轴转完仍然是错的而且错得很有规律——航向误差随横摇角周期性变化看起来像“罗经漂移”实际是坐标系没对齐。3.2 陀螺罗经测的是倾斜面上的航向陀螺罗经走的是另一条物理路线。高速旋转的陀螺转子没有外力矩作用时转轴方向保持不变这叫定轴性有外力矩作用时转轴朝外力矩方向转动这叫进动性。陀螺罗经利用地球自转产生的控制力矩让主轴自动找北再加阻尼力矩让它稳定指北。结构上分三部分灵敏部分是陀螺球负责指北随动部分跟踪陀螺球让主罗经刻度盘的零度线与主轴方向保持一致固定部分除了支承还靠罗经基线来读取航向。关键结论藏在这句话里船舶航向是陀螺主轴和罗经基线在主罗经刻度盘平面上投影线之间的夹角。刻度盘可以绕垂直轴旋转但垂直轴装在固定部分上船一倾斜刻度盘跟着倾斜。所以陀螺罗经读出的是倾斜面上的航向 A不是水平面上的真方位。倾斜改正的做法是在水平系里设陀螺主轴向量为 (xa, ya, 0)在倾斜系里设为 (xb, yb, zb)令 xb sinA、yb cosA用旋转矩阵反算回水平系得到 xa、ya再求真方位 A1最后加上偏转角 ε 得到 A。只考虑横摇、纵摇为零时改正量有一个很整齐的规律横摇角 α最大航向倾斜改正出现航向5°约 0.11°45°、135°、225°、315°10°约 0.44°同上15°约 1.00°同上这几个数字值得记一下。改正量大致随横摇角平方增长而且在四个对角航向上最大正北正东反而为零。这意味着同一条测线来回跑去程和回程的航向误差不一样重叠区会出现随航向缓慢变化的错位不是随机噪声平均不掉。3.3 两种罗经混用会把误差直接灌进平面位置光纤罗经输出的已经是水平系真方位前提是倾斜改正做了陀螺罗经输出的是倾斜面航向除非它内部自带倾斜修正。后处理软件如果默认所有罗经都是同一种把陀螺罗经读数直接当水平系真方位用那 1° 的倾斜改正误差就变成水深点平面位置的系统性旋转。1° 在 200 米斜距的边缘波束上横向偏移约 3.5 米比很多测线的成图比例尺允许的误差大一个量级。分离式设备还要多一道工序先静态测定姿态传感器和罗经坐标系的空间关系。姿态传感器的 Yd 轴到底朝哪罗经坐标系本身是否倾斜这两个问题不搞清楚后面所有航向归算都建在沙子上。常见的做法是把船停在平静水域记录一段时间的姿态和航向读数用均值确定两者之间的固定偏转。这一步做完再去做多波束探头的动态偏差测定顺序不能反。4. 声线弯曲与查表插值斯涅耳定律怎么变成一张可复用的表4.1 声速剖面为什么让声线弯多波束声线计算普遍采用二维模式声速只在垂直深度方向上变化水平方向认为是均匀的。这个假设在开阔水域够用在河口或温跃层剧烈的海域会打折扣但作为归算的基本模型它把三维问题降成了一维积分计算量差了整整一个量级。在这个假设下声线服从 Snell 定律sinθ / v 常数θ 是声线切线与垂线的夹角v 是当地声速。换能器表面声速为 v1、初始入射角为 θ1 时任意深度的入射角满足sinθ (v/v1)·sinθ1。声速随深度变小θ 变大声线向下弯声速变大θ 变小声线向上弯。表层声速计读数如果偏差 1 m/s深水边缘波束的弯曲量就能差出几十厘米。4.2 离散分层积分连续形式是两个积分水深D ∫ v·cosθ dt侧向中心距X ∫ v·sinθ dt积分变量是单程旅行时间。实际声速剖面是一串离散采样点积分变成累加D Σ vi·cosθi·ΔtiX Σ vi·sinθi·Δti每一层的 θi 用 Snell 定律从表面声速和当前层声速递推得到。def trace_ray(svp, theta0, t_total, dt0.002): 二维声线追踪。 svp : [(depth, speed), ...] 垂线声速剖面depth 递增 theta0 : 换能器表面入射角(rad) t_total : 单程旅行时(s) 返回 (D, X): 水深、侧向中心距 D X 0.0 t, theta 0.0, theta0 v_surface svp[0][1] while t t_total: step min(dt, t_total - t) v interp_speed(svp, D) # 用当前深度插值声速 D v * np.cos(theta) * step X v * np.sin(theta) * step t step s v_surface * np.sin(theta0) / v # Snell 不变量 theta np.arcsin(np.clip(s, -1.0, 1.0)) return D, Xdt取 2 毫秒是经验值小到一层之内声速近似不变又不至于让循环次数爆炸。interp_speed按当前累计深度在剖面里插值不能用时间或层号直接索引否则声速跳变的地方会算错一整层。返回的 D 是探头到海底的垂直深度X 是波束脚印相对中心波束的侧向偏距两个值后面都要参与坐标转换。4.3 查表加双线性插值把在线计算搬到离线一个 Ping 少则几十、多则上千个波束每个波束都跑一遍上面的循环实时处理扛不住。通行的做法是预计算以波束入射角为行、旅行时间为列把每个网格点对应的水深和偏距算好存成两张表。实测时按每个波束实际的 θ 和 t 在表里插值。def bilinear(table, i, j, lam, phi): 四邻域双线性插值。 lam: 时间方向权重 (0~1)phi: 角度方向权重 (0~1) d1 table[i][j] lam * (table[i1][j] - table[i][j]) d2 table[i][j 1] lam * (table[i1][j 1] - table[i][j 1]) return d1 phi * (d2 - d1)插值精度取决于表的步长。入射角步长取 2°、旅行时间步长取 2 毫秒时线性内插对偏距和水深的影响不超过 5 厘米。原文给了两组具体估算边缘波束取斜距 200 米、入射角步长一半 1°、入射角 75°偏距误差约 3 厘米、水深误差约 1 厘米近中央波束入射角 1°偏距误差接近零、水深误差约 3 厘米。另外声速在一个时间步长内变化 15 m/s、步长 2 毫秒时声线长度的变化小于 1.5 厘米对垂直波束水深的影响也在同一量级。这几个数字合起来说明一件事查表插值本身引入的误差被压在厘米级比起姿态误差和安装偏差它不是主要矛盾。真正要控制的是表的覆盖范围——如果实测入射角超出表的边界插值会退化成外推误差迅速放大。4.4 波束单元入射角不等于波束入射角换能器各波束单元的入射角 γ 以换能器垂直轴为基准船一摇探头跟着摇γ 就不再是相对水平面的角度。必须先用姿态角把 γ 换算成波束入射角 θ再拿 θ 去查表。换算过程要用到分解形式得到的 α、β、b 三个中间量本质是把波束方向向量从探头坐标系经船体系转到水平系再取它与垂直轴的夹角。这一步的误差会被斜距放大。入射角误差 δθ 对成果的影响近似为δD l·sinθ·δθ、δX l·cosθ·δθl 是声线斜距。取 l 200 米、θ 75°、δθ 0.1°水深误差约 0.33 米。边缘波束对入射角误差的敏感度是中央波束的几十倍这也是为什么姿态传感器的精度指标在深水多波束里比测深精度还关键。5. 常见问题排查安装偏差、时延和粗差多波束最容易翻车的五个地方多波束的坑有个共同特征它们不报错。设备状态正常、数据文件完整、软件流程跑通成果就是不对。下面这几条是从实际作业里反复出现的现象倒推出来的每条按“现象、原因、解决”写顺序也大致对应排查顺序。5.1 往返测线对不上问题多半出在偏差测定顺序现象同一条测线往返测量重叠区同一位置的水深差出一条台阶去程整体偏深或偏浅而且偏差沿测线方向基本恒定。原因纵摇综合偏差没有测定或者测定顺序搞反了。纵摇偏差的影响是让同一水深点在航向上前后移动往返测量时两次移动方向相反重叠区就出现系统性错位。如果先测了艏向偏差再测纵摇艏向偏差会让探头扇面整体偏转一个角度同一船位处的两 Ping 不再平行纵摇的判断依据就被污染了。解决严格按横摇、纵摇、艏向的顺序做。横摇在平坦水域用往返测线测因为平坦水域里航向偏差只影响点位、纵摇偏差对各水深点等幅影响水深水深仍然是平的横摇偏差的影响最突出纵摇在有明显坡度的水域用往返测线测航向垂直于斜坡艏向留到最后用两条平行测线同向测量比较重叠带内同一目标点的位置变化。顺序颠倒一次后面两项都要重测。5.2 平坦水域出现周期性条带横摇偏差没改干净现象平坦水域测出来的水下地形呈波浪状条带条带周期和测线间距接近重叠区尤其明显。原因横摇综合偏差残余。横摇偏差让扇面绕船体纵轴整体倾斜一个固定角度中央波束的水深几乎不受影响越靠边缘的水深偏差越大一个 Ping 内的水深剖面就变成一条斜线。相邻测线拼接时两条斜线在重叠区交叉形成条带。解决回到平坦水域做往返测量用同一位置往返两次的边缘波束水深差反算横摇偏差。这一步要在纵摇和艏向之前做完并先行改正因为横摇偏差会改变扇面的形状后面两项测定的几何前提都依赖扇面是正的。横摇偏差测定对水域平坦度的要求最高选一块几百米范围内水深变化不超过一米的地方比什么算法都管用。5.3 重叠带目标错位艏向综合偏差在作怪现象两条平行测线同向测量重叠带里的沉船、暗礁或其他孤立目标在两条测线上位置不一致沿航向错开一个固定量。原因艏向综合偏差。它同时包含多波束探头和电罗经的安装偏差表现为扇面相对船体中线整体偏转一个角度。同向平行测线测量时两条测线的船位不同扇面偏转导致同一目标在两条测线上的相对位置发生偏移往返测量反而看不出来因为去程和回程的偏移方向相反容易误判成纵摇问题。解决选有明显目标点的水域两条平行测线同向测量比较重叠带内同一目标点的航向位置变化反算艏向偏差。没有天然目标点时可以布设人工反射目标。艏向偏差测定应放在纵摇之后因为它依赖扇面形状已经改正的前提。测完以后探头坐标系到水平坐标系的偏转关系就确定了这个关系同时也是姿态传感器坐标系和罗经坐标系之间的桥梁。5.4 定位点和水深点系统性错开先查定位时延现象水深点的平面位置整体沿航向偏移一个量偏移量随船速增大而增大低速时几乎看不出来。原因定位时延。定位设备解算出的位置对应的是某个时刻声波发射是另一个时刻两者不同步。船在动几十毫秒的时延在 8 节航速下就是几十厘米的位移而且方向固定不会因为多测几遍而抵消。解决在综合偏差测定之前先做定位时延的动态测定。常用做法是沿同一测线往返测量一个孤立目标比较目标在两次测量中的位置差时延造成的偏移在往返方向上相反可以直接解算出来。时延没测就做横摇纵摇测定偏差会被时延污染测出来的“综合偏差”里混进了不属于安装误差的成分。5.5 涌浪改正后水深反而更差杠杆臂放大了姿态系统偏差现象涌浪滤波器读数正常潮位也正常改正后的水深却比不改正时波动更大尤其在船艏或船艉附近作业时明显。原因涌浪滤波器安装位置离多波束探头较远姿态传感器的系统偏差被杠杆臂放大。涌浪改正不是简单地把观测值减掉而是要把滤波器位置的上下起伏换算到探头位置先按安装坐标和姿态系统偏差做一次转换再做倾斜转换。原文给过一组估算滤波器横向偏 10 米、涌浪观测值 0.5 米姿态系统偏差各 1° 时改正值可以差出 0.20 米。滤波器离探头越远这个放大越厉害。解决静态测定姿态传感器和罗经的坐标系偏差先改系统偏差再做倾斜转换。船上不同位置的上下起伏本来就不一样船中后部小、船头大这种差异既引起船体倾斜也让不同位置有不同的涌浪值不能拿一个读数全船通用。一体化设备通常坐标一致性较好分离设备必须逐项测。6. 归算链路的设计取舍误差预算为什么必须卡在厘米级把整条链路的误差摊开看会发现一个不太舒服的事实大部分误差项单独看都不大叠在一起却决定了成果能不能用。误差项性质量级处理手段姿态、航向、声速观测误差偶然设备标称精度误差传播律综合偏差测定误差系统波束角大单次照准精度低大量观测取平均接收时段内姿态变化动态80 毫秒内横摇变化 1°边缘波束水深影响可达 1.66 米模型补偿或缩短 Ping 周期声线查表插值计算厘米级控制步长和表边界定位时延系统随航速线性增长动态测定这张表里最值得琢磨的是第三行。一个 Ping 从接收第一个回波到最后一个回波有时间差30 米水深、单边 3 倍水深本文还有配套的精品资源点击获取
