简介:面向水下无人自主航行器(AUV)研究与开发的MATLAB/Simulink仿真程序包,内容详尽,适合船舶海洋工程、机器人控制等方向的研究生、工程师以及相关课程设计学习者作为参考与二次开发基础。程序包以S函数和M文件为核心,覆盖AUV运动建模、控制算法与仿真流程,从Simulink模型搭建到底层C代码均有注释,便于读者理解水下航行器的动态特性并开展仿真实验。压缩包共62个文件,主要包含C源文件、头文件、M脚本、Simulink模型(mdl)、说明文档以及MAT数据文件等,整体仅441KB,结构紧凑、目录清晰,可直接加载模型运行。当前已有1750人学习浏览,因其代码注释详细、配套文件完整,常被作为AUV仿真入门与进阶的实用参考资料。读者可获得可直接运行的仿真模型、S函数示例、底层C源码及中文说明文档,能有效缩短建模与调试周期,是学习水下机器人仿真的不错选择。
1. 为什么AUV仿真一定要在MATLAB/Simulink里先跑通
水下无人自主航行器(AUV)的调试成本不在代码层,而在水池和船时。一套还没定型的控制算法直接灌进嵌入式板卡,出海一次暴露的往往是坐标系定义、数据接口和时序抖动这类基础问题,而这些在桌面仿真阶段用几分钟就能抓出来。MATLAB/Simulink在AUV领域几乎成了默认选项,原因是它把六自由度动力学、控制器设计、传感器建模、故障注入和最后的C代码生成全部放在同一套仿真程序里,共用一个时间基准和数据结构。这篇笔记写给运动控制工程师、水下机器人方向的研究生,以及从陆地移动平台转水下方向的开发者。下面按“动力学建模—闭环控制—参数整定—数据联调—部署验证”的顺序展开,每一段都给出能直接落地的Simulink模块与MATLAB脚本。
2. AUV六自由度动力学模型:从运动方程到Simulink模块
AUV仿真程序的地基是六自由度(6-DOF)动力学方程。水下环境比地面多出一层“流体耦合”:附加质量让被推动的那部分水也参与运动,水动力阻尼随速度非线性增长,重力和浮力的差值又决定了平台的静稳定性。在这层基础没打正之前,控制器做得再复杂,跑出来的曲线也只是在自己骗自己。
2.1 两个坐标系与运动学/动力学方程
在Simulink里搭AUV,第一步不是拖模块,而是把坐标系定清楚。大地坐标系(NED,North-East-Down)描述位置与姿态,船体坐标系(BODY)描述速度、力和力矩。DVL测速、IMU角速度、推力器输出都在BODY系,任务文件里的路径点、超短基线定位在NED系。两套坐标通过姿态角换算:滚转φ、俯仰θ、偏航ψ。运动学方程写作:
η̇ = J(η) ν
其中η = [x, y, z, φ, θ, ψ]ᵀ是NED系下的位置与欧拉角,ν = [u, v, w, p, q, r]ᵀ是BODY系下的线速度与角速度。J(η)是块对角矩阵,位置部分用方向余弦矩阵,姿态部分用欧拉角率变换矩阵。这个方程不含质量信息,只负责坐标系之间的转换。
动力学部分常用的框架是Fossen形式:
M ν̇ + C(ν)ν + D(ν)ν + g(η) = τ + τ_dist
M包含刚体质量/惯量与附加质量,因为物体在水中加速时还要带动周围一坨水一起运动;C(ν)来自科氏力和向心力,低速机动时多数团队会简化;D(ν)是水动力阻尼,同时包含线性项和与速度平方相关的二次项;g(η)是重力和浮力共同产生的恢复力与恢复力矩,鱼雷形AUV的重心一般低于浮心,这个压力差决定了它横滚和纵倾方向的自恢复特性;τ是推进器在BODY系下给出的控制力/力矩,τ_dist是海流与波浪的等效扰动,故障注入也走这个通道。
2.2 用MATLAB Function实现六自由度方程
常见做法是把上述方程写进一个MATLAB Function模块,外部只留一个参数结构体p。这样换船型时只需替换初始化脚本,模型结构完全不动。下面是六自由度方程的一个可运行骨架,放在Simulink的MATLAB Function模块里即可:
function [nu_dot, eta_dot] = auv6dof(tau, eta, nu, p) % AUV六自由度动力学,Fossen形式,适用于中低速小型AUV % 输入: tau 6x1 控制力/力矩; eta 6x1 位姿; nu 6x1 速度 % 输出: nu_dot 船体系加速度; eta_dot 大地系位姿变化率 phi = eta(4); theta = eta(5); psi = eta(6); cP = cos(phi); sP = sin(phi); cT = cos(theta); sT = sin(theta); cY = cos(psi); sY = sin(psi); % 方向余弦矩阵 R: BODY -> NED R = [cY*cT, cY*sT*sP-sY*cP, cY*sT*cP+sY*sP; sY*cT, sY*sT*sP+cY*cP, sY*sT*cP-cY*sP; -sT, cT*sP, cT*cP]; % 欧拉角率变换矩阵 J_ang,theta=±90°时奇异,工程上应限制俯仰角范围 J_ang = [1, sP*tan(theta), cP*tan(theta); 0, cP, -sP; 0, sP/cT, cP/cT]; eta_dot = [R*nu(1:3); J_ang*nu(4:6)]; % 惯性矩阵:刚体 + 附加质量,按对角处理 M = diag([p.m+p.Xud, p.m+p.Yvd, p.m+p.Zwd, ... p.Ixx+p.Kpd, p.Iyy+p.Mqd, p.Izz+p.Nrd]); % 科氏/向心力矩阵,只保留刚体部分 C = zeros(6); C(1,6) = -p.m*nu(2); C(2,6) = p.m*nu(1); C(4,6) = -p.Izz*nu(6); C(5,6) = p.Iyy*nu(5); C(6,6) = 0; % 线性 + 二次阻尼,方向始终阻碍运动 D = diag([p.Xu+p.Xuu*abs(nu(1)), p.Yv+p.Yvv*abs(nu(2)), ... p.Zw+p.Zww*abs(nu(3)), p.Kp+p.Kpp*abs(nu(4)), ... p.Mq+p.Mqq*abs(nu(5)), p.Nr+p.Nrr*abs(nu(6))]); % 恢复力/力矩,假设重心在原点、浮心在(0,0,zB),zB为负 W = p.m * p.g; B = p.rho * p.V; g_eta = [(W-B)*sT; -(W-B)*cT*sP; -(W-B)*cT*cP; -p.zB*B*cT*sP; -p.zB*B*sT; 0]; nu_dot = M \ (tau - C*nu - D*nu - g_eta); end这段代码里,M被简化成对角阵,忽略横荡与转艏之间的附加质量耦合,对定深直航、稳态回转这类场景精度足够;如果做水下对接或急转弯,必须补全附加质量非对角项,否则仿真的横荡响应会明显偏乐观。阻尼矩阵把线性项和二次项放在同一对角阵,每一项都乘当前速度绝对值,保证阻尼力始终与运动反向。恢复力矩基于重心在原点、浮心在竖轴上的假设,不同文献对g(η)符号的定义存在差异,模型联调前第一件事是自检:松开控制后,AUV的纵倾和横滚应自动回到平衡点,行为不对就翻转相应符号。
注意:欧拉角率变换矩阵在俯仰角接近±90°时奇异。做大俯仰机动(比如海底爬坡)时,姿态部分应改用四元数姿态表示,避免仿真中途发散。
2.3 水动力参数表与初值设置
动力学方程里每一个系数都有明确的物理单位,参数表是仿真程序的“户口本”。下面是一组小型鱼雷形AUV的入门参考值,适合用它跑通第一条闭环曲线:
| 参数 | 符号 | 示例值 | 单位 | 物理说明 |
|---|---|---|---|---|
| 质量 | m | 30 | kg | AUV空气中净质量 |
| 附加质量 | X_udot | -10 | kg | 纵荡方向 |
| 附加质量 | Y_vdot / Z_wdot | -30 / -30 | kg | 横荡/垂荡方向 |
| 转动惯量 | Ixx / Iyy / Izz | 0.5 / 1.5 / 1.5 | kg·m² | 船体系惯量 |
| 线性阻尼 | Xu / Yv / Zw | -3 / -12 / -12 | kg/s | 低速线性阻尼 |
| 二次阻尼 | Xuu / Yvv / Zww | -15 / -30 / -30 | kg/m | 高速平方阻尼 |
| 排水体积 | V | 0.03 | m³ | 用于计算浮力 |
| 密度 | rho | 1000 | kg/m³ | 海水约1025 |
符号约定上,附加质量和阻尼系数在Fossen体系里通常取负值,加进M和D以后实际效果是“减少等效质量”和“消耗能量”。这张表是入门起点,不是标准答案。水动力系数最可靠的来源是CFD计算加约束模试验,公开论文的数据作为初值可以,但不能直接拿去做控制增益设计。我一般把整组参数放进初始化脚本init_auv_params.m里写成结构体p,Simulink模型的每个块只保留p作为参数引用,换船只需要替换脚本,模型本身不用改。
3. 搭建AUV的Simulink闭环:控制器、执行器与传感器
动力学模型立住之后,把它包进子系统,就可以开始搭控制闭环。第一次接触AUV仿真的人容易把控制器堆得很复杂,结果跑出来的曲线反而不如一个带抗饱和的离散PID。水下系统的特点是模型参数不准、测量信号脏,控制结构保持简洁,同时把抗饱和、噪声注入和故障接口留好,比追求先进算法重要得多。
3.1 闭环总体结构与MATLAB Function模块
从顶层看,AUV仿真程序至少分三层:制导层负责给期望航路点或深度剖面,控制层根据位姿误差计算六维控制量τ,执行器层把τ映射到各推进器转速并施加饱和。Simulink里往往把制导、控制、动力学、传感器各自封装成独立子系统,信号用向量或总线连接,避免跨层直接连线。
每个子系统内部建议优先用MATLAB Function模块,而不是堆积基础Simulink块。原因是M,C,D,g这类矩阵运算在MATLAB语言里表达最紧凑,改动参数后代码可读性最好。控制层一个很重要的约定是:输入用错误信号或状态信号,输出统一为六维τ(三个力、三个力矩),这样动力学子系统接口不用动,换控制器只替换一个模块。
3.2 离散PID控制器与抗饱和
控制器在嵌入式平台上的实际形态是离散的,仿真里也用固定步长离散PID更接近实装。下面是一个深度+航向双通道PID的MATLAB Function实现,输出直接合并为六维控制量:
function tau = pid_depth_heading(z_err, psi_err, p) % 深度与航向双通道离散PID % z_err: 深度误差,NED系下向下为正 % psi_err: 航向误差,需预先折叠到[-pi, pi] persistent ei_z ei_psi eprev_z eprev_psi if isempty(ei_z) ei_z = 0; ei_psi = 0; eprev_z = 0; eprev_psi = 0; end dt = p.dt; % 积分项累加,带输出饱和冻结 tau_z_raw = p.Kp_z*z_err + p.Ki_z*ei_z + p.Kd_z*(z_err-eprev_z)/dt; tau_yaw_raw = p.Kp_yaw*psi_err + p.Ki_yaw*ei_psi + p.Kd_yaw*(psi_err-eprev_psi)/dt; if abs(tau_z_raw) < p.tau_z_limit ei_z = ei_z + z_err * dt; end if abs(tau_yaw_raw) < p.tau_yaw_limit ei_psi = ei_psi + psi_err * dt; end tau = zeros(6,1); tau(3) = tau_z_raw; tau(6) = tau_yaw_raw; eprev_z = z_err; eprev_psi = psi_err; end这里的抗饱和逻辑是:当控制量逼近执行器限幅时冻结积分累加,防止积分饱和造成深度的长时间过冲。微分项是噪声放大器,真实传感器信号进控制器之前必须先过低通滤波,否则由测量噪声引起的微分尖峰会在仿真里被误判为控制器性能问题。航向误差折叠到[-π, π]这一步不能省,否则目标航向从-179°变到179°时会给出一个绕远路的控制量。
PID初始增益可以从工程经验起步:深度通道Kp按“重力恢复力能够自然压住”的数量级去试,航向通道Kp先从较小值往上加,Ki取Kp的十分之一到五分之一。跑通后再用Simulink的PID Tuner或Control System Toolbox做局部优化。
3.3 深度/航向滑模控制器
PID在有海流常值扰动时能靠积分消除稳态误差,但对参数突变和短时强扰动的响应偏软。希望提升抗扰能力时,我一般会在同一个Simulink模型里再放一个滑模控制器模块,和PID切换对比。深度通道的滑模控制器实现如下:
function tau_z = smc_depth(err, derr, p) % 深度滑模控制:切换面 s = derr + lambda*err % err为深度误差,derr为深度误差变化率 s = derr + p.lambda * err; % 等效项 + 趋近项,tanh代替sign抑制抖振 tau_z = p.m * (p.zddot_ref - p.lambda*derr) - p.eta * tanh(s/p.epsilon); end切换面参数lambda决定误差收敛带宽,eta必须大于扰动项的上界,epsilon控制边界层厚度。epsilon取太小时,控制量会在滑模面附近高频抖振,仿真步长稍大就会出现“锯齿形”推力曲线;取太大则退化成高增益线性控制。实际操作中从epsilon等于0.05开始,先看控制量曲线是否平滑,再逐步减小。滑模控制的优势在于对模型误差有一定宽容度,适合在仿真里验证“参数偏差20%时控制器是否依然稳定”这类问题。
3.4 传感器模型与噪声注入
传感器模型是仿真程序和实机之间最重要的一道桥梁。AUV常见的传感器组合是深度计、DVL(多普勒测速仪)和IMU,它们的噪声特性完全不同。简化传感器模型如下:
function [z_dvl, z_depth] = auv_sensor(eta, nu, p) % 深度计与DVL测量模型 % eta(3)是NED系深度z,向下为正 % DVL输出BODY系对地速度 z_depth = eta(3) + sqrt(p.var_depth)*randn() + p.bias_depth; z_dvl = nu(1:3) + sqrt(p.var_dvl)*randn(3,1); end传感器参数建议按下面这张表设置,它决定了闭环仿真的可信度:
| 传感器 | 更新率 | 噪声方差 | 标定偏置 | 备注 |
|---|---|---|---|---|
| 深度计 | 10 Hz | 0.0004 m² | 0.1 m | 每100个周期用随机游走更新偏置 |
| DVL | 5 Hz | 2.5e-5 (m/s)² | 0 | 用随机丢包模拟海底地形失锁 |
| IMU陀螺 | 100 Hz | 1e-6 (rad/s)² | 0.01 rad/s | 角速度积分前需去偏置 |
丢包可以单独用一个Bernoulli Random Number模块生成触发信号,DVL在丢包周期内保持上一帧输出,也就是零阶保持。控制器微分项对DVL速度噪声格外敏感,仿真里如果PID输出在静止状态抖动明显,先检查噪声方差是否给得过大,再去动增益。
4. 参数整定、数据导入导出与联合仿真
动力学、控制器和传感器都进模型后,仿真程序才真正进入“工程化”阶段。这个阶段的三个高频需求是批量扫参找稳定区间、把实航或水池实验数据导入仿真做对照、以及把Simulink模型和液压、电机或通信设备联起来验接口。
4.1 用脚本批量扫参数并筛选结果
手动改一次Kp跑一次仿真,效率太低。批量扫参的标准做法是把仿真放进MATLAB脚本里循环调用,用sim命令控制模型运行,把结果写入CSV后统一筛选。下面是深度通道Kp与Ki的粗扫描脚本:
% 批量扫PID参数,结果写入scan_result.csv Kp_range = linspace(10, 60, 6); Ki_range = linspace(0.1, 1.5, 5); scan_data = []; for i = 1:numel(Kp_range) for j = 1:numel(Ki_range) p.Kp_depth = Kp_range(i); p.Ki_depth = Ki_range(j); out = sim('auv_sim_6dof.slx', 'StopTime', '120'); z_err = out.logsout.get('z_err').Values.Data; scan_data = [scan_data; Kp_range(i), Ki_range(j), ... rms(z_err), max(abs(z_err))]; end end writematrix(scan_data, 'scan_result.csv');脚本里的评价指标建议固定为四列:RMS误差、最大绝对误差、调节时间、超调量。RMS误差说明整体跟踪质量,最大绝对误差决定是否会发生触底风险。扫描范围先取大步长粗搜,找到稳定区后再缩小范围细扫。数据量大的时候把for循环改成parsim并行仿真,R2020a及以上版本在并行计算工具箱下可以直接复用当前工作区参数。
提示:
out.logsout的信号名必须和模型中信号记录配置完全一致。跑扫描前先在命令行用simOut.getAllSignals检查一遍,避免脚本中含错误信号名导致整个循环白跑。
4.2 实测CSV导入与FFT振荡分析
仿真数据只有和实测数据对得上才有意义。把实航日志或水槽实验的深度数据导出成CSV,然后在MATLAB里读入,用FFT看振荡频点,是排查控制器异常最直接的手段:
% 导入CSV:第一列为时间,第二列为深度 raw = readmatrix('run_20250506.csv'); t = raw(:,1); z = raw(:,2); Fs = round(1 / mean(diff(t))); z = z - mean(z); Y = fft(z); f = (0:length(Y)-1) * Fs / length(Y); figure; plot(f, abs(Y)); xlim([0 3]); grid on; xlabel('Frequency (Hz)'); ylabel('Magnitude');关注频谱图上能量集中的频点。如果0.2 Hz附近出现明显峰值,而模型特征频率分析显示这不是纵摇或深度环的固有模态,那大概率是传感器延迟和控制器增益一起引入的自激振荡。对照方法是在Simulink里给同一工况施加相同初始扰动,导出仿真深度再跑一遍相同FFT,看峰值频率是否一致。频率对不上,先查模型里的传感器更新率设置和PID的离散步长,这两个参数最容易造成仿真与实测频差。
4.3 联合仿真、CAN故障注入与外部模式
AUV仿真程序很少单独存在。推进器液压系统用AMEsim、水下机械臂用Adams、整车级别的任务规划有时还要和Carsim这类平台做联合仿真,Simulink通过S-Function或FMU标准包把AUV本体模型作为被控对象嵌进更大系统,接口上只需约定输入为六维控制量、输出为位姿与速度向量。仿真程序在设计时就应该把这两个端口做在模型最外层,不要埋在子系统深处。
CAN报文层面的故障诊断仿真越来越常见。推进器、舵机、深度计这类节点在水下用CAN总线通信,Simulink里可以用Vehicle Network Toolbox的CAN Transmit/CAN Receive模块收发真实报文,再在总线中间插入故障注入子系统,把某个节点的帧周期拉长或篡改状态位。这样验证的就不是“控制器在动力学上是否稳定”,而是“通信异常后诊断逻辑能否快速隔离故障”。
外部模式(External Mode)是Simulink把仿真程序拖到实时目标机上的第一步。用Simulink工具栏的“Run on Target”将模型编译部署到目标机,在宿主机界面仍然可以实时修改Kp、Ki,曲线立即回传。这个模式下编辑器里的参数修改不会打乱实时任务,适合在实验室水池里做半实物联调。
模型参数固化推荐用Simulink数据字典(.sldd)。它把m、Xu、Kp这类参数从工作区挪进独立文件,避免脚本切换时工作区变量被覆盖。创建方式:
dictObj = Simulink.data.dictionary.create('auv_params.sldd'); addParameter(dictObj, 'm', 30); addParameter(dictObj, 'Xud', -10); saveChanges(dictObj);有人问“Simulink怎么生成sdf文件”,实际要生成的就是这个.sldd数据字典。模型链接数据字典后,工作区同名变量不再参与仿真,所有块统一从数据字典取值,团队协作和版本管理都会清晰很多。数据字典配合Embedded Coder生成C代码时,参数会作为可配置宏或全局变量导出,实装调试时只改头文件,不必动模型。
5. 从Simulink到实装验证的3个技巧
仿真程序跑得再漂亮,最终都要回答一个问题:这套参数放到真机动平台上还能不能站住。这里分享三个从模型往实装过渡时最有用的技巧,都基于前面已经搭好的仿真程序。
5.1 用外部模式做在线调参
外部模式不只用于实时目标机,有时在水池联调阶段也会直接用。做法是先把仿真模型切成固定步长离散求解器,步长和最终嵌入式任务周期保持一致,比如0.01秒。然后将控制器参数Kp、Ki映射到数据字典,在Simulink编辑界面连接目标机后,在线增大Ki观察深度误差的收敛速度,出现等幅振荡就把Ki调回三分之一。这一步能快速排除“算法看起来对,但控制周期内算不完”这一类实装问题,因为外部模式强制按步长实时执行,模型在宿主机上的运行速度不再掩盖计算耗时。
5.2 用数据字典与C代码生成锁定参数
实装阶段最怕有人偷偷改了一个模型参数,导致下水版本的动力学特性和仿真对不上。数据字典在这里的作用是锁定基线。模型编译前使用Embedded Coder生成C代码,把水动力参数和控制器增益作为#define宏导出到头文件,编译脚本检查头文件哈希值是否与仿真验证版本一致。我一般还会把求解器类型强制设为离散,因为连续求解器代码里会引入大量的ODE求解逻辑,实装平台不一定扛得住,固定步长离散求解器生成的代码结构和仿真行为最接近。
5.3 蒙特卡洛鲁棒性验证
单一工况下的“完美曲线”没有意义,AUV的参数本身就存在不确定性。用蒙特卡洛把控制器放在一整套随机扰动里跑,才能看出它是不是真的稳。下面是最小化的鲁棒性验证脚本:
success = 0; trials = 100; for k = 1:trials p.m = 30 * (1 + 0.1*randn()); p.Xu = -3 * (1 + 0.2*randn()); p.Kp_depth = p.Kp_depth * (1 + 0.1*randn()); out = sim('auv_sim_6dof.slx', 'StopTime', '60'); z_end = out.logsout.get('z_err').Values.Data(end-500:end); if max(abs(z_end)) < 0.2 success = success + 1; end end disp(success / trials);脚本里每次迭代同时扰动质量、阻尼和控制器增益,模拟水动力系数辨识偏差与控制增益漂移的叠加效果。成功率低于90%时,先不要急着换控制器结构,回头检查一下参数不确定性范围是否给得过于乐观。把蒙特卡洛结果写成报告附在评审材料里,比单条时域曲线有说服力得多。
本文还有配套的精品资源,点击获取