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

WW、forever

已于 2025-07-14 19:22:38 修改

阅读量261

点赞数 7

CC 4.0 BY-SA版权

分类专栏： WRFDA原理及实现文章标签： matlab 数据结构开发语言

于 2025-07-14 19:19:29 首次发布

本文链接：https://blog.csdn.net/qq_44246618/article/details/149337257

WRFDA原理及实现专栏收录该内容

23 篇文章

订阅专栏

模型概述

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

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

MATLAB实现代码

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

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

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

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

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

%% ETKF codes
% 参数-状态联合估计的同化系统，使用ETKF（Ensemble Transform Kalman Filter）方法对 Lorenz-63 模型的状态变量和参数进行估计
clear all
close all

pathFigure = "..\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 步，生成一个真实轨道 xt
  x = 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:mem
  x = randn(nx,1);
  for i=1:1500
    x = tinteg_rk4(modelname,pa(:,m,1),x,dt);
  end
  xa(:,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 );
  end
  xt(:,it+1) = x;

  % observations 生成观测值
  % 对真实状态加上观测误差（ober为观测误差方差），生成模拟观测
  y = h * xt(:,it+1) + ober * randn(ny,1);

  % first guess 集合预测（forward prediction）
  % 每个成员根据自己的参数预测未来状态
  for m=1:mem
    x = xa(:,m,it);
    for i=1:itassim
      x = 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);
  % 扰动变换矩阵 T
  T = sqrt(mem-1)*V*sqrt(inv(D))*transpose(V);

  % 更新扩展状态集合 za
  za = 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);         % k4

xp = 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