当前位置：首页 > ds >正文

【数据同化案例1】ETKF求解参数-状态联合估计的同化系统（完整MATLAB实现）

ds 2025/7/15 5:23:54

模型概述

Lorenz 1963 模型的详细介绍可参见另一博客-常见混沌系统：Lorenz 1963 模型。

在这里插入图片描述
此案例采用Lorenz 1963 模型作为状态变量。

MATLAB实现代码

三维状态变量的 3D 轨迹图如下所示：
在这里插入图片描述

参数估计结果如下：（随同化次数增多，参数趋于稳定）
在这里插入图片描述

相应的状态估计过程如下：（实线表示估计值，虚线表示真实值；后续估计值与真实值几乎一致）
在这里插入图片描述

需要注意，先验数据随机生成，因而每次的结果不尽相同。

完整主函数如下：（包括图形绘制代码）

%% ETKF codes
% 参数-状态联合估计的同化系统，使用ETKF（Ensemble Transform Kalman Filter）方法对 Lorenz-63 模型的状态变量和参数进行估计
clear all
close allpathFigure = "..\Figures\";%% =============
%%% 1 . Settings
%%%=============% 使用 Lorenz-63 模型，有三个参数（a, r, b）
modelname='lorenz63';% true model parameters
np = 3;
a = 10;
r = 28;
b = 8/3;
pt = zeros(np,1);
pt(1) = a;
pt(2) = r;
pt(3) = b;% model settings
nx = 3;           % 状态变量个数（x, y, z）
dt = 0.01;       % 积分时间步长% observation settings  观测设置
ny = nx;
h = eye(ny);                             % 每个状态变量都被观测（h=I）：单位矩阵
ober = 0.1;                               % 观测误差标准差为 0.1
Rinv = eye(ny)/ober/ober;       % 观测误差协方差矩阵的逆% assimilation settings  同化设置
itassim = 5;            % 每次同化间隔进行 5 步积分
infl = 1.04;              % infl 是膨胀因子inflation coefficient%% 初始化（first guess）% experimental period (assimilation steps)
nt = 1000;              % nt=1000 是同化总步数% nature run: xt   生成真实轨道（nature run）
xt = zeros(nx , nt);
x = randn(nx , 1);
for i=1:1500% 多次调用：连续积分 1500 步，生成一个真实轨道 xtx = tinteg_rk4(modelname, pt , x, dt);
end
xt(:,1) = x;                  % xt 存储真实状态变量% initial parameter ensemble 初始化参数和状态集合（ensemble）
mem = 5;                                % 集合成员个数
pa = zeros( np, mem, nt);       % pa 是分析后的参数集合
pf = zeros( np, mem, nt);        % pf 是预测的参数集合
for i=1:np% 初始参数集合设置为真实值加上偏移和噪声（人为设偏差增加挑战）% pa(i,:,1) = pt(i)*ones(1,mem,1) + pt(i)/5*randn(1,mem,1);pa(i,:,1) = pt(i)*ones(1,mem,1)+ 3*ones(1,mem,1)+3*randn(1,mem,1);
end% initial condition ensemble
% 每个成员都使用自己的参数跑 1500 步，初始化状态集合
xa = zeros(nx, mem, nt);
xf = zeros(nx, mem, nt);
for m=1:memx = randn(nx,1);for i=1:1500x = tinteg_rk4(modelname,pa(:,m,1),x,dt);endxa(:,m,1) = x;
end% 构造扩展状态向量（Augmented State Vectors）
% 使用 Kalman 滤波方法 同时估计参数和状态
zf  = zeros(nx+np,mem);           % 预测时刻的扩展状态集合，每列是一个集合成员的 [x; p] 组合
zfm = zeros(nx+np,1);               % zfm 是 zf 的集合均值，所有成员扩展状态的平均值（用于计算扰动矩阵）
dzf = zeros(nx+np,mem);          % 扰动矩阵（deviation from mean），表示每个成员与集合均值的偏离
za  = zeros(nx+np,mem);          % 分析时刻的扩展状态集合，ETKF 更新后计算得到，每列是更新后的 [x; p]
hz  = zeros(ny,nx+np);               % 扩展观测矩阵/扩展观测算子
hz(:, 1:nx) = h;            % 对状态使用 h
hz(:, nx+1:end) = 0;   % 对参数无观测（观测不涉及参数部分）%% ===========================
%%% 2 . Main assimilation loop 同化主循环
%%%===========================% nt=1000 是同化总步数（共同化1000-1次）
for it=1:nt-1% nature run 真实轨道前向积分x = xt(:,it);for i=1:itassim% 使用真实参数积分 5 步，生成下一个真实状态x = tinteg_rk4( modelname, pt , x , dt );endxt(:,it+1) = x;% observations 生成观测值% 对真实状态加上观测误差（ober为观测误差方差），生成模拟观测y = h * xt(:,it+1) + ober * randn(ny,1);% first guess 集合预测（forward prediction）% 每个成员根据自己的参数预测未来状态for m=1:memx = xa(:,m,it);for i=1:itassimx = tinteg_rk4(modelname,pa(:,m,it),x,dt);end% 存储预每个集合预测成员的状态和参数xf(:,m,it+1) = x;pf(:,m,it+1) = pa(:,m,it);end% ETKF assimilation 数据同化zf(1   :nx   ,:) = xf(:,:,it+1);           % 预测时刻的扩展状态集合zf(nx+1:nx+np,:) = pf(:,:,it+1);   % 预测时刻的扩展状态集合zfm = mean(zf,2);                      % zfm 是平均扩展状态dzf = zf - zfm*ones(1,mem);     % dzf 是扰动矩阵% 计算协方差矩阵并进行特征值分解（用于ETKF分析增益计算）% 加上 infl（膨胀因子（inflation coefficient）） 是为了防止滤波器崩溃[V ,D ] = eig((mem-1)*eye(mem)/infl + transpose(hz*dzf)*Rinv*hz*dzf);% ETKF 更新公式：更新后的扩展状态集合 za% 均值偏移 δdelta = V*inv(D)*transpose(V)*transpose(hz*dzf)*Rinv*(y-hz*zfm)*ones(1,mem);% 扰动变换矩阵 TT = sqrt(mem-1)*V*sqrt(inv(D))*transpose(V);% 更新扩展状态集合 zaza = zfm*ones(1,mem) + dzf * ( delta + T);% 包含状态和参数的更新xa(:,:,it+1) = za(1:nx,:);                     % 状态更新pa(:,:,it+1) = za(nx+1:nx+np,:);        % 参数更新
end%% ============
%%% 3 . plot
%%%============%% 绘制 3D 轨迹图T = size(xt, 2);                   % 获取轨迹长度
colormapName = hot;       % 可以改为 jet、hot、cool、turbo 等% 创建颜色映射（T 行的 RGB，每一行是一个颜色）
cmap = colormap(colormapName);
nColors = size(cmap, 1);
colorIdx = round(linspace(1, nColors, T));  % 将时间映射到颜色索引
colors = cmap(colorIdx, :);                 % 每个点的颜色% 绘图
figure('Name','Lorenz Attractor','Color','w');
scatter3(xt(1,:), xt(2,:), xt(3,:), 8, colors, 'filled'); % 使用颜色渐变
grid on;xlabel('x(t)', 'FontSize', 14, 'FontName', 'Times New Roman');
ylabel('y(t)', 'FontSize', 14, 'FontName', 'Times New Roman');
zlabel('z(t)', 'FontSize', 14, 'FontName', 'Times New Roman');
title(sprintf('Lorenz Attractor: σ=%.1f, ρ=%.1f, β=%.3f', pt(1), pt(2), pt(3)));
view(30, 20);
axis tight;
set(gca, 'FontSize', 12, 'FontName', 'Times New Roman');% 保存图像
str = strcat(pathFigure, "Fig.1 3D 轨迹图", '.tiff');
print(gcf, '-dtiff', '-r600', str);%% 绘制 参数估计图
anal=squeeze(mean(pa(:,:,:),2));
true=pt(:)*ones(1,nt);figure
hold on;box on;% 估计值
plot(anal(1,:),'b', 'LineWidth',1.2)
plot(anal(2,:),'g', 'LineWidth',1.2)
plot(anal(3,:),'r', 'LineWidth',1.2)% 真实值
plot(true(1,:),'b--', 'LineWidth',1.5)
plot(true(2,:),'g--', 'LineWidth',1.5)
plot(true(3,:),'r--', 'LineWidth',1.5)
set(gca, 'FontSize', 12, 'FontName', 'Times New Roman');title('Parameter estimation')
xlabel('Assimilation steps','FontSize',14,'FontName','Times New Roman')
ylabel('Values','FontSize',14,'FontName','Times New Roman')
hl = legend('a','r','b');
set(hl, 'NumColumns', 3, 'box', 'off', 'FontName', 'Times New Roman', 'FontSize', 15,'Location','north');
set(gca,'Layer','top');str= strcat(pathFigure, "Fig.2 parameter estimation", '.tiff');
print(gcf, '-dtiff', '-r600', str);%% 绘制 状态同化图figureUnits = 'centimeters';
figureWidth = 27; 
figureHeight = 12;figureHandle = figure;
set(gcf, 'Units', figureUnits, 'Position', [0 0 figureWidth figureHeight]);anal=squeeze(mean(xa(:,:,:),2));
true=xt;
plot(anal(1,:),'b')
hold on
plot(anal(2,:),'g')
plot(anal(3,:),'r')
plot(true(1,:),'b--')
plot(true(2,:),'g--')
plot(true(3,:),'r--')
set(gca, 'FontSize', 12, 'FontName', 'Times New Roman');title('State estimation','FontSize',14,'FontName','Times New Roman')
xlabel('Assimilation steps','FontSize',14,'FontName','Times New Roman')
ylabel('Values','FontSize',14,'FontName','Times New Roman')
hl = legend('x','y','z');
set(hl, 'NumColumns', 3, 'box', 'off', 'FontName', 'Times New Roman', 'FontSize', 15,'Location','north');
set(gca,'Layer','top');str= strcat(pathFigure, "Fig.3 state estimation", '.tiff');
print(gcf, '-dtiff', '-r600', str);

四阶龙格-库塔（RK4）数值积分方法的实现函数如下：
函数 tinteg_rk4 实现了一个 四阶龙格-库塔（RK4） 数值积分方法，用于数值求解微分方程（ODE）。它是整个程序中对动力学模型（如 Lorenz-63）进行时间积分的核心函数。

function xp = tinteg_rk4(modelname, p, x, dt)
%% 四阶龙格-库塔（RK4）数值积分方法，用于数值求解微分方程（ODE）
% modelname：模型名称字符串（如 'lorenz63'），用于调用模型函数；
% p：模型参数（如 [a, r, b]）；
% x：当前状态变量（如 [x; y; z]）；
% dt：时间步长；
% 输出：xp 是在时间步 dt 后的状态。q1 = dt * feval(modelname, p, x);                  % k1
q2 = dt * feval(modelname, p, x + q1*0.5);   % k2
q3 = dt * feval(modelname, p, x + q2*0.5);   % k3
q4 = dt * feval(modelname, p, x + q3);         % k4xp = x + (q1 + 2*q2 + 2*q3 + q4)/6;             % x_{n+1}end

RK4 方法是一种经典的显式积分方法，适用于常微分方程（ODE）：
在这里插入图片描述

Loren模型定义如下：

function xp = lorenz63(p,x)
% Lorenz-63 模型的右端项，用于积分a = p(1);
r = p(2);
b = p(3);xdot = a*(x(2) - x(1));
ydot = -x(1)*x(3) + r*x(1) - x(2);
zdot = x(1)*x(2) - b*x(3);xp = [xdot;ydot;zdot];end