课题学习(八)----卡尔曼滤波动态求解倾角、方位角

文章详细介绍了卡尔曼滤波在钻柱动力学中的应用,包括KF-1和KF-2两种滤波器的设计,用于处理加速度计和磁强计数据,计算重力加速度并平滑井眼轨迹。文中涉及状态空间模型、测量噪声处理和动态角度计算方法。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

一、 卡尔曼滤波

   卡尔曼滤波的应用要求系统和底层过程的测量模型都是线性的。离散时间线性状态空间系统的描述为: x k = Φ k , k − 1 x k − 1 + G k − 1 w k − 1 x_k=\Phi_{k,k-1}x_{k-1}+G_{k-1}w_{k-1} xk=Φk,k1xk1+Gk1wk1
   式中 Φ k , k − 1 \Phi_{k,k-1} Φk,k1为状态转移矩阵, x k x_k xk为状态向量, G k − 1 G_{k-1} Gk1为噪声分布矩阵, w k − 1 w_{k-1} wk1为过程噪声向量,k为测量历元。
   系统的测量方程由下式给出 z k = H k x k + η k z_k=H_kx_k+\eta_k zk=Hkxk+ηk
   式中 z k z_k zk为系统输出的测量向量, H k H_k Hk为观测或设计, η k \eta_k ηk为测量噪声。系统噪声w和测量噪声 η k \eta_k ηk是具有确定的自协方差函数的非相关零均值白噪声过程。
   通过对钻柱动力学的分析,提出了一种新的状态空间模型算法。在计算方法中,我们定义KF-1的输入矢量为 X = [ a x a y a z m x m y m z ] X=\begin{bmatrix}a_x&a_y&a_z\\ m_x&m_y&m_z\end{bmatrix} X=[axmxaymyazmz],由三轴磁强计和三轴加速度计测量。KF-2的输入向量为井眼倾角和方位角,定义为 X = [ I A ] X=\begin{bmatrix}I\\ A\end{bmatrix} X=[IA],其中 I I I为井眼倾角, A A A为井眼方位角。
   求解过程如下图所示:
在这里插入图片描述
   在此过程中,设计两个卡尔曼滤波器。在KF-1之后,我们可以得到更精确的重力加速度信号gx, gy和gz,定义为KF-1的输出。利用重力加速度,通过建立钻柱旋转时的方程,可以得到钻柱的倾角和方位角。然后使用KF-2进一步平滑钻井轨迹。KF-2的输出定义为 M ′ = [ I ′ A ′ ] T M'=\begin{bmatrix}I'&A'\end{bmatrix}^T M=[IA]T,这是更精确的。

1.1 KF-1的状态空间模型

   传感器安装在钻柱的中心,在旋转过程中,x轴和y轴的测量信号呈现正弦波。理论上,加速度计和磁力计的信号有相同的规律。在实际钻井过程中,钻柱的振动对磁强计信号的影响较小。也就是说,磁力计信号是用来校准加速度计信号的。通过实验室测试,我们可以得出磁通门信号的变化与重力加速度信号的变化是一致的。
   假设角速度为 w x , y , z w_{x,y,z} wx,y,z,采样间隔为Δt,则
在这里插入图片描述
   在KF-1中,我们将状态向量定义为 X = [ g x ( k ) g y ( k ) g z ( k ) ] X=\begin{bmatrix}g_x{(k)}\\g_y{(k)}\\g_z{(k)}\end{bmatrix} X= gx(k)gy(k)gz(k) ,从 a x , a y , a z a_x,a_y,a_z ax,ay,az的振动信号可以得到重力加速度信号 x g ( k ) x_g{(k)} xg(k)。因此,系统输出的测量矢量为,当钻柱旋转时,除了z轴信号外,x轴和y轴的测量信号将呈现正弦波。因此,变换矩阵定义为 z ( k ) = [ a x a y a z ] z(k)=\begin{bmatrix}a_x\\a_y\\a_z\end{bmatrix} z(k)= axayaz , H k = [ m x ( k − 1 ) m x ( k ) 0 0 0 m y ( k − 1 ) m y ( k ) 0 0 0 1 ] H_k=\begin{bmatrix}\frac{m_x{(k-1)}}{m_x{(k)}}&0&0\\ 0&\frac{m_y{(k-1)}}{m_y{(k)}}&0\\ 0&0&1\end{bmatrix} Hk= mx(k)mx(k1)000my(k)my(k1)0001 。系统噪声 w k w_k wk和测量噪声 η k \eta_k ηk是不相关的零均值白噪声过程。因此,我们得到KF-1的状态空间模型如下: x g ( k ) = [ m x ( k − 1 ) m x ( k ) 0 0 0 m y ( k − 1 ) m y ( k ) 0 0 0 1 ] x g ( k − 1 ) + [ w x ( k − 1 ) w y ( k − 1 ) w z ( k − 1 ) ] x_g{(k)}=\begin{bmatrix}\frac{m_x{(k-1)}}{m_x{(k)}}&0&0\\ 0&\frac{m_y{(k-1)}}{m_y{(k)}}&0\\ 0&0&1\end{bmatrix}x_g{(k-1)+\begin{bmatrix}w_x{(k-1)}\\w_y{(k-1)}\\w_z{(k-1)}\end{bmatrix}} xg(k)= mx(k)mx(k1)000my(k)my(k1)0001 xg(k1)+ wx(k1)wy(k1)wz(k1)
z ( k ) = H k x g ( k ) + η ( k ) z(k)=H_kx_g{(k)}+\eta(k) z(k)=Hkxg(k)+η(k)

1.2 计算倾角和方位角

   钻柱旋转时,上述方程不适用。安装在钻柱中心的传感器,即x轴和y轴的测量信号,在旋转过程中呈现正弦波。
   通过KF-1,我们得到gx、gy和gz,它们分别定义为x、y和z轴上的重力加速度测量信号。然后定义系统的输入向量为 X ′ = [ g x g y g z m x m y m z ] X'=\begin{bmatrix}g_x&g_y&g_z\\ m_x&m_y&m_z\end{bmatrix} X=[gxmxgymygzmz]
   定义转速为R,若R为0,则用下式计算倾角和方位角。
在这里插入图片描述
   R ≠ 0 R\neq0 R=0表示钻柱在旋转,倾角I和方位角A的计算公式如下:
在这里插入图片描述
   式中 T M T_M TM为磁性工具面角: T M = t g − 1 ( − m y m x ) T_M=tg^-1(\frac{-m_y}{m_x}) TM=tg1(mxmy)。用于近直井的工具面角。磁性工具面是在垂直于井筒轴线的平面上,顺时针方向测量的井眼测量仪器在井筒内的角度或方位角;北、东、南、西方向的磁性工具面角分别为0°、90°、180°和270°。磁性工具面可以校正为参考栅格北或真北。
   虽然钻柱转速是判断钻柱是否旋转的一种方法,但这种方法的可靠性不高。相反,使用标准差统计方法来确定钻柱的运动更加可靠,因为它反映了群体中个体之间的分散程度。使用50个数据点作为时间窗口,假设它们是x1, x2,…, x49, x50,我们得到标准差 σ = 1 N ∑ i = 1 N ( x i − x ˉ ) 2 \sigma=\sqrt{\frac{1}{N}\sum_{i=1}^N(x_i-\bar{x})^2} σ=N1i=1N(xixˉ)2 ,当标准差σ接近于零时,钻柱可以认为是静止的。

1.3 KF-2的状态空间模型

   在钻孔过程中,运动状态为 σ = 0 \sigma=0 σ=0 σ ≠ 0 \sigma\neq0 σ=0交替出现。当 σ ≠ 0 \sigma\neq0 σ=0时,由于钻柱不旋转时振动较小,求解结果更为准确。我们开发了另一种卡尔曼滤波器(KF-2)来平滑钻井轨迹。(KF-2分别分为KF-2.1和KF2.2,如下图所示。)
在这里插入图片描述
   在正常的钻井过程中,为了测量倾角和方位角,钻井作业经常要停在测量站。然后根据数学假设计算两个测量站之间的井眼轨迹。例如,可以假设钻孔距离为直线、圆弧或多角线;每个都需要不同的计算方法。设实际钻井轨迹测量第N点的三维坐标为 ( x N , y N , z N ) (x_N,y_N,z_N) (xN,yN,zN),则测量(N +1)点为 ( x N + 1 , y N + 1 , z N + 1 ) (x_{N+1},y_{N+1},z_{N+1}) (xN+1,yN+1,zN+1),井深、垂深、倾角、方位角分别为 L N 、 H N 、 θ N 、 Ψ N L_N、H_N、\theta_N、\Psi_N LNHNθNΨN和$L_{N+1}、H_{N+1}、\theta_{N+1}、\Psi_{N+1}。两点之间的井眼轨迹定义如下:
在这里插入图片描述
   如果垂直深度H已知,则测量(N + 1)-th点的三维坐标可定义为
在这里插入图片描述
   我们可以通过递归计算得到每个点的空间坐标,从而得到整个钻井轨迹。如下图所示,我们可以使用井眼轨迹外推法建立相邻两个测点之间的递推关系。
在这里插入图片描述

   假设L是钻孔深度, γ \gamma γ是钻孔轨迹与其切线的夹角,然后: L ( k ) = L ( k − 1 ) + Δ L ( K ) L(k)=L(k-1)+\Delta L(K) L(k)=L(k1)+ΔL(K)
γ = a r c c o s [ c o s ( I k − 2 ) c o s ( I k − 1 ) + s i n ( I k − 2 ) s i n ( I k − 1 ) c o s ( A k − 1 − A k − 2 ) ] \gamma=arccos[cos(I_{k-2})cos(I_{k-1})+sin(I_{k-2})sin(I_{k-1})cos(A_{k-1}-A_{k-2})] γ=arccos[cos(Ik2)cos(Ik1)+sin(Ik2)sin(Ik1)cos(Ak1Ak2)]
γ ( k ) = γ Δ L ( k − 1 ) Δ L ( k ) \gamma(k)=\frac{\gamma}{\Delta L(k-1)}\Delta L(k) γ(k)=ΔL(k1)γΔL(k)
   从上面三个公式看出,我们可以使用两个点来估计下一个点,因此可以平滑井眼轨迹。我们可以使用Kalman Filter 2.2对钻井轨迹进行校准,其中将第k个测量点定义为 P ( k ) = [ I k A k ] P(k)=[I_k A_k] P(k)=[IkAk]。系统输入为P(K -2)和P(K -1), KF-2.2结合实测值和理论计算值估算P(K)。
   KF-2.2通过倾角和方位角作为输入,实现井眼轨迹平滑。假设状态向量为 x ( k ) = [ I k A k ] x(k)=\begin{bmatrix}I_k\\A_k\end{bmatrix} x(k)=[IkAk],则系统输出的测量向量为 x ( k ) = [ I k A k ] m x(k)=\begin{bmatrix}I_k\\A_k\end{bmatrix}_m x(k)=[IkAk]m H ( k ) = [ 1 0 0 1 ] H(k)=\begin{bmatrix}1&0\\0&1\end{bmatrix} H(k)=[1001]。得到KF-2.2的状态空间模型:
在这里插入图片描述
   如上所示,我们应该确定 Δ L \Delta L ΔL作为KF-2.2的输入。( Δ t \Delta t Δt期间钻头向前移动的距离)。我们可以通过测量z轴上的加速度来计算位移。 A z A_z Az是z轴上三轴加速度计的信号,它与重力加速度和振动加速度相结合。定义加速度在z轴时间序列上的测量为 a z ( k ) a_z(k) az(k)
   所以在计算位移 Δ L \Delta L ΔL之前,我们首先要排除重力的影响,如下所示: f g z ( k ) = a z ( k ) − G ⋅ c o s ( I k − 1 ) f_{g_z}(k)=a_z(k)-G·cos(I_{k-1}) fgz(k)=az(k)Gcos(Ik1)
   其中 f g z ( k ) f_{g_z}(k) fgz(k)是加速度时间序列函数,通过去掉 g z g_z gz的加速度,可以计算出 a z ( k ) a_z(k) az(k)和倾角 I k − 1 I_{k-1} Ik1对应的同一时间。
   然后我们可以使用z轴时间序列 f g z ( k ) f_{g_z}(k) fgz(k)上的加速度来计算钻探深度( Δ L \Delta L ΔL)。将状态向量定义为 Δ L ( k ) = [ Δ L ( k ) Δ L ˙ ( k ) Δ L ¨ ( k ) ] \Delta L(k)=\begin{bmatrix}\Delta L(k)\\\Delta \dot L(k)\\\Delta \ddot L(k)\end{bmatrix} ΔL(k)= ΔL(k)ΔL˙(k)ΔL¨(k) ,系统输出的测量向量为 z ( k ) = f g z ( k ) z(k)=f_{g_z}(k) z(k)=fgz(k),我们建立了KF-2.1的状态空间模型如下:
在这里插入图片描述
   通过对实测信号进行预处理,建立了一种基于底部旋转钻具动力学分析的动态方位和倾角求解算法。在理论模型的基础上,我们开发了一种卡尔曼滤波器来提高求解器的精度;在预测钻井轨迹的基础上,建立了卡尔曼滤波的状态空间方程。基于卡尔曼滤波的动态测量算法是一种能够大大降低求解误差的新模型。

二、 往期回顾

课题学习(一)----静态测量
课题学习(二)----倾角和方位角的动态测量方法(基于磁场的测量系统)
课题学习(三)----倾角和方位角的动态测量方法(基于陀螺仪的测量系统)
课题学习(四)----四元数解法
课题学习(五)----阅读论文《抗差自适应滤波的导向钻具动态姿态测量方法》
课题学习(六)----安装误差校准、实验方法
课题学习(七)----粘滑运动的动态算法

下面的代码运行中显示以下问题:科里奥利力矩阵以及更新机器人动力学 (真实状态)部分错误使用 * 用于矩阵乘法的维度不正确。请检查并确保第一个矩阵中的列数第二个矩阵中的行数匹配。要单独对矩阵的每个元素进行运算,请使用 TIMES (.*)执行按元素相乘。 classdef UnderwaterPositioningSystem < handle % 水下机器人位置保持系统 % 包含超声定位、EKF滤波和自适应控制器 properties % 超声探头参数 r = 0.5; % 探头安装半径(m) beta = deg2rad(60); % 探头倾角(60度) alpha = deg2rad([0, 120, 240]); % 探头方位角 % EKF参数 x_est; % 状态估计 [dx, dy, vx, vy, theta, phi, wx, wy]^T P_est; % 估计协方差矩阵 Q; % 过程声协方差 R; % 测量声协方差 % 控制器参数 Kp = 2.0; % 比例增益 Ki = 0.1; % 积分增益 Kd = 3.0; % 微分增益 integral_error = [0; 0]; % 积分误差项 prev_error = [0; 0]; % 上一次误差 % 机器人物理参数 mass = 50; % 质量(kg) inertia = [10, 0; 0, 10]; % 转动惯量(kg·m²) drag_coeff = 1.5; % 阻尼系数 thruster_gain = 0.8;% 推进器增益 % 流体动力学参数 added_mass = diag([15, 15, 20, 5, 5, 10]); % 附加质量矩阵 quadratic_drag = diag([0.8, 0.8, 1.2, 0.3, 0.3, 0.5]); % 二次阻尼 linear_drag = diag([1.5, 1.5, 2.0, 0.5, 0.5, 1.0]); % 线性阻尼 % 目标位置 target_position = [0; 0]; % [dx_target; dy_target] % 历史记录 position_history = []; true_position_history = []; control_history = []; time_history = []; thruster_history = []; % 推进器指令历史记录 error_history = []; % 位置误差历史记录 % 时间步长 dt = 0.1; % 控制周期(100ms) % 真实状态 (用于仿真) true_position = [0; 0]; true_velocity = [0; 0]; true_orientation = [0; 0]; true_ang_velocity = [0; 0]; end methods function obj = UnderwaterPositioningSystem() % 构造函数 - 初始化系统 % 初始化EKF状态 (dx, dy, vx, vy, theta, phi, wx, wy) obj.x_est = zeros(8, 1); % 初始化协方差矩阵 (增大不确定性) obj.P_est = diag([0.5, 0.5, 0.2, 0.2, 0.1, 0.1, 0.1, 0.1]); % 过程声协方差 (模型不确定性) obj.Q = diag([0.05, 0.05, 0.1, 0.1, 0.02, 0.02, 0.05, 0.05]); % 测量声协方差 (传感器) obj.R = diag([0.1, 0.1, 0.02, 0.02]); % dx, dy, theta, phi % 初始化历史记录 obj.position_history = zeros(2, 1000); obj.true_position_history = zeros(2, 1000); obj.control_history = zeros(4, 1000); obj.thruster_history = zeros(4, 1000); % 推进器指令历史 obj.time_history = zeros(1, 1000); obj.error_history = zeros(2, 1000); % 位置误差历史 end function [theta, phi, dx, dy] = ultrasound_positioning(obj, d, orientation) % 改进的超声定位解算,考虑姿态影响 % 输入:d - 三个探头的测量距离 [d1, d2, d3] % orientation - 当前姿态 [theta, phi] % 输出:姿态和水下位置 % 获取当前姿态 pitch = orientation(1); roll = orientation(2); % 旋转矩阵 (ZYX欧拉角) R = [cos(pitch)*cos(roll), -cos(pitch)*sin(roll), sin(pitch); sin(roll), cos(roll), 0; -sin(pitch)*cos(roll), sin(pitch)*sin(roll), cos(pitch)]; % 计算探头方向向量 (考虑姿态) v = @(i) R * [sin(obj.beta)*cos(obj.alpha(i)); sin(obj.beta)*sin(obj.alpha(i)); -cos(obj.beta)]; % 探头安装位置 p = @(i) [obj.r*cos(obj.alpha(i)); obj.r*sin(obj.alpha(i)); 0]; % 计算水底点坐标 P1 = p(1) + d(1) * v(1); P2 = p(2) + d(2) * v(2); P3 = p(3) + d(3) * v(3); % 拟合水底平面 vec1 = P2 - P1; vec2 = P3 - P1; n = cross(vec1, vec2); n_hat = n / norm(n); % 解算姿态角 theta = asin(-n_hat(1)); phi = atan2(n_hat(2), -n_hat(3)); % 计算水平偏移量 O_prime = dot(n_hat, P1) * n_hat; dx = O_prime(1); dy = O_prime(2); end function ekf_predict(obj, u) % EKF预测步骤 % u - 控制输入 [Fx; Fy; M_theta; M_phi] % 状态变量: x = [dx, dy, vx, vy, theta, phi, wx, wy]^T % 从状态向量中提取当前值 dx = obj.x_est(1); dy = obj.x_est(2); vx = obj.x_est(3); vy = obj.x_est(4); theta = obj.x_est(5); phi = obj.x_est(6); wx = obj.x_est(7); wy = obj.x_est(8); % 控制输入 Fx = u(1); Fy = u(2); M_theta = u(3); M_phi = u(4); % 状态转移函数 (非线性动力学模型) % 位置更新 dx_new = dx + vx * obj.dt; dy_new = dy + vy * obj.dt; % 速度更新 (考虑阻尼和水流) vx_new = vx + (Fx - obj.drag_coeff * vx)/obj.mass * obj.dt; vy_new = vy + (Fy - obj.drag_coeff * vy)/obj.mass * obj.dt; % 姿态更新 theta_new = theta + wx * obj.dt; phi_new = phi + wy * obj.dt; % 角速度更新 wx_new = wx + (M_theta - obj.drag_coeff * wx)/obj.inertia(1,1) * obj.dt; wy_new = wy + (M_phi - obj.drag_coeff * wy)/obj.inertia(2,2) * obj.dt; % 更新状态预测 obj.x_est = [dx_new; dy_new; vx_new; vy_new; ... theta_new; phi_new; wx_new; wy_new]; % 改进的雅可比矩阵计算 (包含耦合项) F = eye(8); F(1,3) = obj.dt; % ∂dx/∂vx F(2,4) = obj.dt; % ∂dy/∂vy F(3,3) = 1 - (obj.drag_coeff/obj.mass)*obj.dt; % ∂vx/∂vx F(4,4) = 1 - (obj.drag_coeff/obj.mass)*obj.dt; % ∂vy/∂vy F(5,7) = obj.dt; % ∂theta/∂wx F(6,8) = obj.dt; % ∂phi/∂wy F(7,7) = 1 - (obj.drag_coeff/obj.inertia(1,1))*obj.dt; % ∂wx/∂wx F(8,8) = 1 - (obj.drag_coeff/obj.inertia(2,2))*obj.dt; % ∂wy/∂wy % 添加耦合项 F(3,5) = -sin(theta)*obj.dt; % 俯仰角对X加速度的影响 F(4,6) = -sin(phi)*obj.dt; % 横滚角对Y加速度的影响 % 更新协方差矩阵 obj.P_est = F * obj.P_est * F' + obj.Q; end function ekf_update(obj, z) % EKF更新步骤 % z - 测量值 [dx_meas; dy_meas; theta_meas; phi_meas] % 观测矩阵 (测量是状态的一部分) H = [1 0 0 0 0 0 0 0; 0 1 0 0 0 0 0 0; 0 0 0 0 1 0 0 0; 0 0 0 0 0 1 0 0]; % 计算卡尔曼增益 S = H * obj.P_est * H' + obj.R; K = obj.P_est * H' / S; % 更新状态估计 innovation = z - H * obj.x_est; obj.x_est = obj.x_est + K * innovation; % 更新协方差矩阵 obj.P_est = (eye(8) - K * H) * obj.P_est; end function u = adaptive_controller(obj) % 改进的自适应PID控制器 % 返回控制输入 u = [Fx; Fy; M_theta; M_phi] % 从状态估计中提取位置和姿态 dx_est = obj.x_est(1); dy_est = obj.x_est(2); vx_est = obj.x_est(3); vy_est = obj.x_est(4); theta_est = obj.x_est(5); phi_est = obj.x_est(6); % 位置误差 error_pos = [dx_est; dy_est] - obj.target_position; obj.error_history(:, end+1) = error_pos; % 记录误差 % 姿态误差 (目标姿态为水平) error_att = [theta_est; phi_est]; % 计算误差范数和变化率 error_norm = norm(error_pos); error_rate = norm([vx_est; vy_est]); % 自适应增益调整 if error_norm > 0.5 % 大误差区域 - 高增益 Kp_adaptive = obj.Kp * 2.0; Ki_adaptive = obj.Ki * 0.5; % 限制积分 Kd_adaptive = obj.Kd * 1.5; elseif error_norm > 0.1 % 中等误差区域 Kp_adaptive = obj.Kp * 1.2; Ki_adaptive = obj.Ki * 0.8; Kd_adaptive = obj.Kd * 1.0; else % 小误差区域 - 低增益防震荡 Kp_adaptive = obj.Kp * 0.8; Ki_adaptive = obj.Ki * 1.0; Kd_adaptive = obj.Kd * 0.8; end % 添加误差变化率影响 if error_rate > 0.3 Kd_adaptive = Kd_adaptive * min(2.0, 1 + error_rate); end % 位置控制 (PID) obj.integral_error = obj.integral_error + error_pos * obj.dt; % 积分项限幅防饱和 max_integral = 5.0; % 积分上限 obj.integral_error(1) = sign(obj.integral_error(1)) * min(abs(obj.integral_error(1)), max_integral); obj.integral_error(2) = sign(obj.integral_error(2)) * min(abs(obj.integral_error(2)), max_integral); % 力计算 Fx = -Kp_adaptive * error_pos(1) - Ki_adaptive * obj.integral_error(1) - Kd_adaptive * vx_est; Fy = -Kp_adaptive * error_pos(2) - Ki_adaptive * obj.integral_error(2) - Kd_adaptive * vy_est; % 姿态-位置解耦补偿 Fx_comp = Fx * cos(theta_est) - Fy * sin(phi_est); Fy_comp = Fx * sin(theta_est) + Fy * cos(phi_est); % 姿态控制 (PD) M_theta = -Kp_adaptive * 0.5 * error_att(1) - Kd_adaptive * obj.x_est(7); M_phi = -Kp_adaptive * 0.5 * error_att(2) - Kd_adaptive * obj.x_est(8); % 限制最大控制量 max_force = 50; % N max_torque = 20; % N·m Fx_comp = sign(Fx_comp) * min(abs(Fx_comp), max_force); Fy_comp = sign(Fy_comp) * min(abs(Fy_comp), max_force); M_theta = sign(M_theta) * min(abs(M_theta), max_torque); M_phi = sign(M_phi) * min(abs(M_phi), max_torque); % 更新误差记录 obj.prev_error = error_pos; % 控制输出 u = [Fx_comp; Fy_comp; M_theta; M_phi]; end function thruster_commands = thruster_allocation(obj, u) % 推进器分配 % 输入:u = [Fx; Fy; M_theta; M_phi] % 输出:推进器指令 [T1, T2, T3, T4] % 假设有4个推进器:前、后、左、右 % 分配矩阵 % [ Fx ] [ 0 0 1 -1 ] [ T1 ] % [ Fy ] = [ 1 -1 0 0 ] [ T2 ] % [ M_theta] [ l l 0 0 ] [ T3 ] % [ M_phi ] [ 0 0 l l ] [ T4 ] l = 0.8; % 力臂长度(m) % 构建分配矩阵 A = [0, 0, 1, -1; 1, -1, 0, 0; l, l, 0, 0; 0, 0, l, l]; % 求解推进器指令 thruster_commands = A \ u; % 应用推进器增益和饱和限制 max_thrust = 30; % N thruster_commands = obj.thruster_gain * thruster_commands; thruster_commands = min(max(thruster_commands, -max_thrust), max_thrust); end function S = skew(~, v) % 计算叉乘的斜对称矩阵 S = [0, -v(3), v(2); v(3), 0, -v(1); -v(2), v(1), 0]; end function update_robot_dynamics(obj, u, flow_disturbance) % 改进的机器人动力学模型,包含流体动力学效应 % 输入:控制输入u,水流干扰flow_disturbance % 控制力 Fx = u(1); Fy = u(2); M_theta = u(3); M_phi = u(4); % 从状态中提取当前值 vel = [obj.true_velocity; 0]; % 添加Z方向速度(设为0) ang_vel = [0; 0; obj.true_ang_velocity]; % 添加偏航角速度(设为0) % 构建广义速度向量 (6自由度) nu = [vel; ang_vel]; % 科里奥利力矩阵 C = [zeros(3) -obj.mass*obj.skew(ang_vel); -obj.mass*obj.skew(ang_vel) -obj.skew(obj.added_mass(4:6,4:6)*ang_vel)]; % 阻尼力 (线性 + 二次) D = obj.linear_drag * abs(nu) + obj.quadratic_drag * (nu .* abs(nu)); % 广义力 (只考虑水平方向) F_generalized = [Fx; Fy; 0; 0; M_theta; M_phi] + [flow_disturbance; 0; 0; 0; 0]; % 质量矩阵 (包含附加质量) M = diag([obj.mass, obj.mass, obj.mass, ... obj.inertia(1,1), obj.inertia(2,2), obj.inertia(1,1)]) + obj.added_mass; % 加速度计算 (考虑附加质量和阻尼) acc = M \ (F_generalized - C*nu - D*nu); % 状态更新 (只更新水平方向) obj.true_velocity = obj.true_velocity + acc(1:2)*obj.dt; obj.true_ang_velocity = obj.true_ang_velocity + acc(5:6)*obj.dt; obj.true_position = obj.true_position + obj.true_velocity*obj.dt; obj.true_orientation = obj.true_orientation + obj.true_ang_velocity*obj.dt; end function simulate(obj, total_time, flow_conditions) % 系统仿真 % total_time: 总仿真时间(s) % flow_conditions: 水流干扰条件 % 初始化真实状态 obj.true_position = [0.5; -0.3]; % 初始偏移 obj.true_velocity = [0; 0]; obj.true_orientation = [deg2rad(1.5); deg2rad(-1.0)]; obj.true_ang_velocity = [0; 0]; % 初始化EKF状态 obj.x_est(1:2) = obj.true_position; % 初始位置偏移 obj.x_est(5:6) = obj.true_orientation; % 初始姿态 % 初始化控制器状态 obj.integral_error = [0; 0]; obj.prev_error = [0; 0]; % 主循环 num_steps = round(total_time / obj.dt); time = 0; for k = 1:num_steps % 记录时间 time = time + obj.dt; obj.time_history(k) = time; % 生成水流干扰 (随时间变化) flow_disturbance = [0.5 * sin(0.3*time) + flow_conditions(1); 0.3 * cos(0.2*time) + flow_conditions(2)]; % 控制器计算 u = obj.adaptive_controller(); % 推进器分配 thruster_cmds = obj.thruster_allocation(u); obj.thruster_history(:, k) = thruster_cmds; % 记录推进器指令 % 更新机器人动力学 (真实状态) obj.update_robot_dynamics(u, flow_disturbance); % 模拟超声测量 () measurement_noise = 0.02; % 测量声水平 % 计算真实距离 (考虑姿态影响) true_depth = 10; % 水深10米 true_distances = true_depth * ones(3,1); for i = 1:3 % 考虑姿态对测量的影响 attitude_factor = 1 + 0.1 * sin(obj.true_orientation(1)) * cos(obj.alpha(i)) ... + 0.1 * sin(obj.true_orientation(2)) * sin(obj.alpha(i)); true_distances(i) = true_depth / attitude_factor; end d_meas = true_distances + measurement_noise * randn(3,1); % 超声定位 (使用改进的算法) [theta_meas, phi_meas, dx_meas, dy_meas] = ... obj.ultrasound_positioning(d_meas, obj.true_orientation); % EKF更新 z = [dx_meas; dy_meas; theta_meas; phi_meas]; obj.ekf_update(z); % 记录位置和控制历史 obj.position_history(:, k) = obj.x_est(1:2); obj.true_position_history(:, k) = obj.true_position; obj.control_history(:, k) = u; % EKF预测 (使用上一个控制输入) obj.ekf_predict(u); % 显示状态 if mod(k, 10) == 0 fprintf('时间: %.1fs, 估计X偏移: %.3fm, 真实X偏移: %.3fm, 估计Y偏移: %.3fm, 真实Y偏移: %.3fm\n', ... time, obj.x_est(1), obj.true_position(1), obj.x_est(2), obj.true_position(2)); end end % 截断历史记录 obj.time_history = obj.time_history(1:num_steps); obj.position_history = obj.position_history(:, 1:num_steps); obj.true_position_history = obj.true_position_history(:, 1:num_steps); obj.control_history = obj.control_history(:, 1:num_steps); obj.thruster_history = obj.thruster_history(:, 1:num_steps); obj.error_history = obj.error_history(:, 1:num_steps); % 绘制结果 obj.plot_results(); % 分析性能 obj.analyze_performance(); end function plot_results(obj) % 绘制仿真结果 figure('Position', [100, 100, 1200, 800], 'Name', '水下机器人位置保持系统仿真结果'); % 位置轨迹 subplot(2, 2, 1); plot(obj.time_history, obj.position_history(1, :), 'b-', 'LineWidth', 1.5); hold on; plot(obj.time_history, obj.true_position_history(1, :), 'r--', 'LineWidth', 1.5); plot(obj.time_history, obj.position_history(2, :), 'g-', 'LineWidth', 1.5); plot(obj.time_history, obj.true_position_history(2, :), 'm--', 'LineWidth', 1.5); plot(obj.time_history, zeros(size(obj.time_history)), 'k--', 'LineWidth', 1.0); title('水下机器人位置偏移'); xlabel('时间 (s)'); ylabel('偏移量 (m)'); legend('估计X偏移', '真实X偏移', '估计Y偏移', '真实Y偏移', '目标位置'); grid on; % 控制力 subplot(2, 2, 2); plot(obj.time_history, obj.control_history(1, :), 'b', 'LineWidth', 1.5); hold on; plot(obj.time_history, obj.control_history(2, :), 'r', 'LineWidth', 1.5); title('控制力'); xlabel('时间 (s)'); ylabel('力 (N)'); legend('X方向控制力', 'Y方向控制力'); grid on; % 姿态角 subplot(2, 2, 3); plot(obj.time_history, rad2deg(obj.x_est(5) * ones(size(obj.time_history))), 'b', 'LineWidth', 1.5); hold on; plot(obj.time_history, rad2deg(obj.x_est(6) * ones(size(obj.time_history))), 'r', 'LineWidth', 1.5); plot(obj.time_history, rad2deg(obj.true_orientation(1)) * ones(size(obj.time_history)), 'b--', 'LineWidth', 1.0); plot(obj.time_history, rad2deg(obj.true_orientation(2)) * ones(size(obj.time_history)), 'r--', 'LineWidth', 1.0); plot(obj.time_history, zeros(size(obj.time_history)), 'k--', 'LineWidth', 1.0); title('机器人姿态'); xlabel('时间 (s)'); ylabel('角度 (°)'); legend('估计俯仰角', '估计横滚角', '真实俯仰角', '真实横滚角', '目标姿态'); grid on; % 位置轨迹图 subplot(2, 2, 4); plot(obj.position_history(1, :), obj.position_history(2, :), 'b', 'LineWidth', 1.5); hold on; plot(obj.true_position_history(1, :), obj.true_position_history(2, :), 'r', 'LineWidth', 1.5); plot(0, 0, 'go', 'MarkerSize', 10, 'LineWidth', 2); title('水下机器人运动轨迹'); xlabel('X位置 (m)'); ylabel('Y位置 (m)'); legend('估计轨迹', '真实轨迹', '目标位置'); axis equal; grid on; % 添加标题 sgtitle('水下机器人位置保持系统仿真结果', 'FontSize', 16, 'FontWeight', 'bold'); end function analyze_performance(obj) % 性能分析 % 计算稳态误差 (取后20%数据) steady_state_idx = round(0.8 * length(obj.time_history)); ss_error_x = rms(obj.position_history(1, steady_state_idx:end) - obj.true_position_history(1, steady_state_idx:end)); ss_error_y = rms(obj.position_history(2, steady_state_idx:end) - obj.true_position_history(2, steady_state_idx:end)); % 收敛时间 (误差<0.05m的时间) conv_threshold = 0.05; conv_time_x = find(abs(obj.position_history(1,:) - obj.true_position_history(1,:)) < conv_threshold, 1) * obj.dt; conv_time_y = find(abs(obj.position_history(2,:) - obj.true_position_history(2,:)) < conv_threshold, 1) * obj.dt; % 控制能量消耗 control_energy = sum(sum(obj.control_history.^2)) * obj.dt; fprintf('\n性能分析结果:\n'); fprintf('X方向估计误差: %.4f m\n', ss_error_x); fprintf('Y方向估计误差: %.4f m\n', ss_error_y); fprintf('X方向收敛时间: %.2f s\n', conv_time_x); fprintf('Y方向收敛时间: %.2f s\n', conv_time_y); fprintf('控制能量消耗: %.2f N²·s\n', control_energy); % 绘制估计误差 figure('Name', '位置估计误差分析'); subplot(2,1,1); plot(obj.time_history, obj.position_history(1,:) - obj.true_position_history(1,:), 'b'); hold on; plot(obj.time_history, obj.position_history(2,:) - obj.true_position_history(2,:), 'r'); title('位置估计误差'); xlabel('时间 (s)'); ylabel('误差 (m)'); legend('X方向误差', 'Y方向误差'); grid on; subplot(2,1,2); plot(obj.time_history, abs(obj.position_history(1,:) - obj.true_position_history(1,:)), 'b'); hold on; plot(obj.time_history, abs(obj.position_history(2,:) - obj.true_position_history(2,:)), 'r'); title('位置估计绝对误差'); xlabel('时间 (s)'); ylabel('绝对误差 (m)'); legend('X方向', 'Y方向'); grid on; end end end
最新发布
07-23
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

致虚守静~归根复命

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值