单眼平滑追踪

[Reza-111]

为什么图片显示全蓝和全红

[donghaiwang]

还需要贴出绘图的脚本和说明绘制的数据是什么

[Reza-111]

%% 基于RNN+眼球肌肉控制的2D平滑追踪（目标定中央凹）
% 功能：整合论文FORCE-RNN预测模型+眼球肌肉动力学+2D追踪+肌肉疲劳
% 修正：去掉sat函数，用max/min内联饱和逻辑（兼容所有MATLAB版本）

clear; clc; close all;
rng(0); % 固定随机种子，保证结果可复现

%% ====================== 1. 核心参数初始化 ======================
% ---------------------- 1.1 时间参数（与论文一致） ----------------------
dt = 0.016; % 积分步长(16ms，论文FORCE学习步长)
t_total = 20; % 总模拟时间(s)
t = 0:dt:t_total; % 时间向量
Nt = length(t); % 时间步数
delay_vis = 0.08; % 视觉处理延迟(80ms，论文值)
delay_step = round(delay_vis / dt); % 延迟步数

% ---------------------- 1.2 目标轨迹参数（2D正弦运动） ----------------------
freq_h = 0.5; % 水平方向频率(Hz)
amp_h = 10; % 水平运动幅度(deg，中央凹±10deg)
freq_v = 0.3; % 垂直方向频率(Hz)
amp_v = 8; % 垂直运动幅度(deg)

% ---------------------- 1.3 眼球肌肉动力学参数（生理贴合） ----------------------
% 肌肉最大转动范围(deg)：水平(内/外直肌)、垂直(上/下直肌)
K_m_h = 45; K_m_v = 40;
% 肌肉时间常数(s)：激活延迟/惯性/阻尼
tau_act = 0.05; tau_inert = 0.02; zeta = 0.7;
% 肌肉疲劳参数：疲劳系数（随时间增加，肌肉激活效率下降）
fatigue_rate = 0.001; % 疲劳增长速率
fatigue = ones(1, Nt);% 肌肉疲劳因子（初始=1，无疲劳）

% ---------------------- 1.4 PID控制器参数（反馈校正） ----------------------
Kp = 18; Kd = 2.5; Ki = 0.6; % 比例/微分/积分增益

% ---------------------- 1.5 FORCE-RNN参数（论文核心） ----------------------
N = 500; % 神经元数量（论文值）
g = 1.5; % 递归权重增益（论文值，保证混沌）
tau_rnn = 0.16; % RNN时间常数(160ms，论文值)
alpha = 100; % FORCE学习正则化常数（论文值）
K = 1; % 反馈增益（论文值）
tau_lpf = 0.128; % 漏积分器时间常数(128ms，论文值)
K_lpf = 1; % 漏积分器增益（论文值）

%% ====================== 2. 生成2D目标轨迹（位置+速度） ======================
% 水平方向（x）：目标位置+速度
P_Tx = amp_h * sin(2pifreq_h*t);
v_Tx = diff(P_Tx)/dt; v_Tx = [v_Tx, v_Tx(end)]; % 补全长度

% 垂直方向（y）：目标位置+速度
P_Ty = amp_v * sin(2pifreq_v*t + pi/4); % 相位偏移，增加2D复杂性
v_Ty = diff(P_Ty)/dt; v_Ty = [v_Ty, v_Ty(end)];

% 合并2D目标速度（用于RNN学习）
v_T = [v_Tx; v_Ty]; % 2×Nt矩阵

%% ====================== 3. FORCE-RNN模型实现（论文核心：预测眼速度） ======================
% 初始化RNN状态
x = zeros(N, 1); % 神经元状态
r = tanh(x); % 神经元激活（tanh）
W = g * randn(N, N) / sqrt(N); % 递归权重（论文初始化方式）
W_out = zeros(2, N); % 读出权重（2维：x/y方向）
W_feedback = randn(N, 2); % 反馈权重（论文：-1~1均匀分布，此处简化为正态）
P = eye(N) / alpha; % FORCE学习的逆相关矩阵

% 初始化RNN输出（预测眼速度）
v_E_rnn = zeros(2, Nt); % 2×Nt：RNN预测的眼速度（x/y）
v_E_lpf = zeros(2, Nt); % 漏积分器后的平滑速度

% FORCE学习迭代（在线学习目标速度）
for k = 2:Nt
% 1. 计算视网膜滑移（RS=目标速度-当前预测速度，论文误差信号）
rs = v_T(:, k) - v_E_lpf(:, k-1);

% 2. 更新RNN神经元状态（论文动力学方程）
dx = (-x + W * r + K * W_feedback * v_E_lpf(:, k-1)) / tau_rnn;
x = x + dx * dt;
r = tanh(x);

% 3. FORCE学习更新权重（论文公式7/8）
k_force = P * r;
rP = r' * P;
P = P - k_force * rP / (1 + rP * r);
W_out = W_out + rs * k_force';

% 4. 读出层输出（预测眼速度）
v_E_rnn(:, k) = W_out * r;

% 5. 漏积分器（低通滤波，论文公式6）
dv_lpf = (-v_E_lpf(:, k-1) + K_lpf * v_E_rnn(:, k)) / tau_lpf;
v_E_lpf(:, k) = v_E_lpf(:, k-1) + dv_lpf * dt;

end

% 积分得到RNN预测的眼球位置指令（前馈用）
P_E_cmd_x = cumsum(v_E_lpf(1, :)) * dt;
P_E_cmd_y = cumsum(v_E_lpf(2, :)) * dt;

%% ====================== 4. 初始化眼球与肌肉状态变量 ======================
% 眼球状态：位置(deg)、速度(deg/s)、加速度(deg/s²)（x/y方向）
P_E = zeros(2, Nt); v_E = zeros(2, Nt); a_E = zeros(2, Nt);

% 肌肉状态：激活水平(0~1)、肌肉力(N)（x/y方向）
u_muscle = zeros(2, Nt); F_muscle = zeros(2, Nt);

% 误差变量：位置误差（目标-眼球）、延迟误差、积分误差
e_P = zeros(2, Nt); e_P_delay = zeros(2, Nt); int_eP = zeros(2, 1);

%% ====================== 5. 眼球肌肉控制迭代（闭环反馈） ======================
for k = 2:Nt
% ---------------------- 步骤1：更新肌肉疲劳因子 ----------------------
fatigue(k) = fatigue(k-1) * exp(-fatigue_rate * dt); % 指数衰减，模拟疲劳
if fatigue(k) < 0.5 % 疲劳下限：肌肉效率不低于50%
fatigue(k) = 0.5;
end

% ---------------------- 步骤2：计算视觉延迟后的位置误差 ----------------------
% 实时位置误差（目标-当前眼球位置）
e_P(1, k) = P_Tx(k) - P_E(1, k-1); % 水平x
e_P(2, k) = P_Ty(k) - P_E(2, k-1); % 垂直y

% 80ms视觉延迟（延迟后的误差用于控制）
if k > delay_step
    e_P_delay(1, k) = e_P(1, k - delay_step);
    e_P_delay(2, k) = e_P(2, k - delay_step);
else
    e_P_delay(:, k) = 0;
end

% ---------------------- 步骤3：PID反馈控制（调整肌肉激活） ----------------------
for dim = 1:2 % 遍历x/y方向
    int_eP(dim) = int_eP(dim) + e_P_delay(dim, k) * dt; % 积分项
    dedt = (e_P_delay(dim, k) - e_P_delay(dim, k-1)) / dt; % 微分项
    u_pid = Kp * e_P_delay(dim, k) + Kd * dedt + Ki * int_eP(dim);
    
    % 前馈项：RNN位置指令→肌肉激活
    if dim == 1
        u_ff = P_E_cmd_x(k) / K_m_h;
    else
        u_ff = P_E_cmd_y(k) / K_m_v;
    end
    
    % 总肌肉激活（含疲劳+饱和限制：用max/min替代sat函数）
    u_muscle(dim, k) = max(min((u_pid + u_ff) * fatigue(k), 1), 0);
end

% ---------------------- 步骤4：肌肉动力学→眼球运动 ----------------------
for dim = 1:2
    % 肌肉激活→肌肉力（F_max=10N，生理最大力）
    F_muscle(dim, k) = 10 * u_muscle(dim, k);
    
    % 肌肉力→眼球加速度（逆动力学：克服惯性、阻尼、弹性）
    a_E(dim, k) = (F_muscle(dim, k) - 2*zeta*tau_inert*v_E(dim, k-1) - P_E(dim, k-1)/tau_inert) / (tau_inert^2);
    
    % 积分更新速度和位置（欧拉法）
    v_E(dim, k) = v_E(dim, k-1) + a_E(dim, k) * dt;
    P_E(dim, k) = P_E(dim, k-1) + v_E(dim, k) * dt;
    
    % 限制眼球转动范围（生理极限：用max/min替代sat函数）
    if dim == 1
        P_E(dim, k) = max(min(P_E(dim, k), K_m_h), -K_m_h);
    else
        P_E(dim, k) = max(min(P_E(dim, k), K_m_v), -K_m_v);
    end
end

end

%% ====================== 6. 可视化追踪效果（核心结果展示） ======================
figure('Position', [100, 100, 1400, 900], 'Color', 'white');

% 子图1：2D目标与眼球位置轨迹
subplot(2,3,1);
plot(P_Tx, P_Ty, 'k--', 'LineWidth', 2, 'DisplayName', '目标轨迹');
hold on;
plot(P_E(1, :), P_E(2, :), 'r-', 'LineWidth', 2, 'DisplayName', '眼球轨迹');
plot(0, 0, 'g+', 'MarkerSize', 12, 'LineWidth', 2, 'DisplayName', '中央凹中心');
xlabel('水平位置 (deg)'); ylabel('垂直位置 (deg)');
title('2D目标与眼球轨迹'); legend; grid on; axis equal;

% 子图2：水平方向位置追踪（x）
subplot(2,3,2);
plot(t, P_Tx, 'k--', 'LineWidth', 2, 'DisplayName', '目标x');
hold on;
plot(t, P_E(1, :), 'r-', 'LineWidth', 2, 'DisplayName', '眼球x');
xlabel('时间 (s)'); ylabel('位置 (deg)');
title('水平方向追踪'); legend; grid on;

% 子图3：垂直方向位置追踪（y）
subplot(2,3,3);
plot(t, P_Ty, 'k--', 'LineWidth', 2, 'DisplayName', '目标y');
hold on;
plot(t, P_E(2, :), 'r-', 'LineWidth', 2, 'DisplayName', '眼球y');
xlabel('时间 (s)'); ylabel('位置 (deg)');
title('垂直方向追踪'); legend; grid on;

% 子图4：中央凹对准误差（x+y）
subplot(2,3,4);
e_total = sqrt(e_P(1, :).^2 + e_P(2, :).^2);
plot(t, e_total, 'b-', 'LineWidth', 1.5);
xlabel('时间 (s)'); ylabel('总对准误差 (deg)');
title('中央凹对准误差（理想=0）');
grid on; ylim([0, max(e_total)+0.5]);

% 子图5：肌肉激活强度（x+y）
subplot(2,3,5);
plot(t, u_muscle(1, :), 'm-', 'LineWidth', 1.5, 'DisplayName', '水平肌肉');
hold on;
plot(t, u_muscle(2, :), 'c-', 'LineWidth', 1.5, 'DisplayName', '垂直肌肉');
xlabel('时间 (s)'); ylabel('肌肉激活水平 (0~1)');
title('眼球肌肉激活强度'); legend; grid on;

% 子图6：动态追踪动画（直观展示中央凹对准）
subplot(2,3,6);
h_target = plot(P_Tx(1), P_Ty(1), 'ko', 'MarkerSize', 8, 'MarkerFaceColor', 'k');
h_eye = plot(P_E(1,1), P_E(2,1), 'ro', 'MarkerSize', 15, 'MarkerFaceColor', 'r');
h_fovea = plot(0, 0, 'g+', 'MarkerSize', 12, 'LineWidth', 2);
xlabel('水平位置 (deg)'); ylabel('垂直位置 (deg)');
title('动态追踪效果（黑=目标，红=中央凹）');
xlim([-amp_h-5, amp_h+5]); ylim([-amp_v-5, amp_v+5]);
grid on; axis equal;

% 动画更新（放慢速度，便于观察）
for k = 1:Nt
set(h_target, 'XData', P_Tx(k), 'YData', P_Ty(k));
set(h_eye, 'XData', P_E(1,k), 'YData', P_E(2,k));
drawnow;
pause(dt/20); % 动画速度：dt/20为慢放系数
end

%% ====================== 7. 量化评估追踪效果（关键指标） ======================
e_x = abs(e_P(1, :)); e_y = abs(e_P(2, :));
e_total = sqrt(e_P(1, :).^2 + e_P(2, :).^2);

fprintf('====================== 追踪效果量化评估 ======================\n');
fprintf('水平方向平均对准误差：%.3f deg（理想<0.5deg）\n', mean(e_x));
fprintf('垂直方向平均对准误差：%.3f deg（理想<0.5deg）\n', mean(e_y));
fprintf('总平均对准误差：%.3f deg\n', mean(e_total));
fprintf('最大对准误差：%.3f deg\n', max(e_total));
fprintf('肌肉平均疲劳程度：%.3f（1=无疲劳，<0.5=重度疲劳）\n', mean(fatigue));