【数据同化案例1】ETKF求解 Lorenz-63 模型的同化系统（完整MATLAB实现）

本博客基于 集合变换卡尔曼滤波ETKF（Ensemble Transform Kalman Filter）求解 Lorenz-63 模型的同化系统。需要注意的是，本博客只同化状态变量。

模型概述

Lorenz-63 模型

Lorenz 1963 模型是描述大气对流的一个简化系统，具有混沌特性，是非线性动力系统研究的经典模型之一。该系统由三个常微分方程组成：

Lorenz-63 模型是一个低维混沌系统，具有著名的蝴蝶形吸引子。其特性包括：

• 对 初始条件 极度敏感（蝴蝶效应）
  

• 存在一个 非线性吸引子 ，即系统的长期行为会限制在某个特定的形状/范围内（而不是无界增长）
  

Lorenz-63 系统的吸引子

吸引子（Attractor）基本定义：

在动力系统中，吸引子是系统长期演化后，状态变量所趋近的集合。它可以是：

• 一个点（稳定平衡）
• 一个周期轨道（极限环）
• 一个分形结构的集合（典型的混沌吸引子，如 Lorenz 吸引子）

对于 Lorenz-63 模型，对应的吸引子是一个著名的 “蝴蝶形”混沌吸引子，是一个在三维相空间中复杂蜿蜒的轨道集合，不是一个点！

吸引子特点：

特征	描述
维度	在 3D 空间中演化（状态变量 x, y, z）
不是一个点	而是无数点组成的轨道集合
依赖参数	不同的 σ、ρ、β 会导致不同类型的吸引子（甚至非混沌）
对初值敏感	但长期都会落入吸引子（吸引性）

Lorenz 系统的行为高度依赖参数，例如：

参数组合	吸引子形态	案例图形
ρ < 1	Lorenz 系统在 𝜌<1 时，原点是唯一稳定点，所有轨道都会收敛到原点，系统不出现振荡或混沌。	轨道从初始点缓慢收敛到原点，轨迹几乎不动。（ρ=0.5<1）
ρ ≈ 10	出现两个稳定平衡点（非混沌），轨道会围绕其中一个平衡点收敛，表现为规则螺旋	轨道围绕两个点之间螺旋收敛，看不到明显的“蝴蝶形”混沌结构。
ρ = 28	出现混沌吸引子（蝴蝶形）	Lorenz 系统最经典的混沌参量，轨道在两个“翅膀”之间不可预测地跳跃，形成典型混沌吸引子。
ρ > 100	系统变得更混乱，吸引子形状更复杂	两个“翅膀”结构非常清晰，但比 𝜌=28 时范围更广、点分布更散，Z 方向振幅超过 200。

如果你从一个随机或不在吸引子上的初始点开始模拟：

• 初期行为不代表系统的真实长期特性；
• 整个模拟轨迹会表现出非典型行为，影响分析；
• 在数据同化中，这种初始点会导致算法性能下降，甚至发散。

Spin-up 初始条件

“Spin-up” 是指在进行正式模拟之前，先运行模型一段时间，使其状态从任意初始值演化到系统吸引子上的过程。这个过程也可以称为“系统热身”。

进行系统预热的代码如下：

x0sut = [1;1;1];   % 初始条件（随意选一个接近吸引子的点）
Tsu = 100;         % spin-up 时间（模拟 100 单位时间）
x0 = x0sut';       % 转置为行向量
Ken = 1;           % 只跑一个集合成员
rho = rhot; beta = betat; sigma = sigmat;[tsu,xsu] = ode45(@lorenz3_dxdt,[0,Tsu],x0);
x0t = xsu(end,:);

1、选择一个初始点 x0sut = [1; 1; 1]：

这个点不一定在吸引子上，但在吸引子附近。

2、设定 spin-up 时间 Tsu = 100：

让 Lorenz 模型自由演化 100 个时间单位，足够让系统落入吸引子区域。

3、使用 ODE45 进行积分：

数值积分微分方程，得到从 t=0 到 t=100 的轨道 xsu。

4、取最终状态作为新初始点：

x0t = xsu(end,:) 是 t=100 时的状态，它已经稳定在吸引子上，可以作为后续模拟或同化的起点。

🔄 为什么只跑一个集合成员？
这里只是为了生成一个“参考轨道”的起点，不需要多个 ensemble；后续数据同化中会构建 ensemble，但这一步只需要一个典型的吸引子上的点。

MATLAB实现代码

三维状态变量的 3D 轨迹图如下：

---

真实值与观测值的误差如下：（白噪声生成）

---

绘制状态变量的真值/估计值图形如下：(x1和x3变量有观测数据)

---

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

%% ETKF codes
% 基于 Lorenz 3 模型（Lorenz 1963 系统）的 Ensemble Transform Kalman Filter (ETKF) 同化过程%% initialization
close all
clear allpathFigure = "..\Figures\";
%% global 声明全局变量
% rho, beta, sigma: Lorenz 系统的参数
% N: 系统状态变量维度（Lorenz 3 系统中为 3）
% Ken: 集合成员数（ensemble size）global rho beta sigma
global N Ken%% sample etkf
fprintf('\n\n sample etkf using lorenz 3 variable: etkf v1\n')%% setup
%- input & intial spinup of the true system
etkf_v1_ip%% true 生成真实系统轨迹（Nature run）
%- setup
%-- ic is computed by spin up in input file
%-- parametersrho=rhot; 
beta=betat; 
sigma=sigmat;%-- time 时间设置
tt = 0:dTt:Tend;  
KTt = size(tt,2);     % 长度，即length(tt)%-- output
xt=zeros(N,KTt);
kn=1; xt(:,kn)=x0t;
for k=2:KTt;  kn=k-1; x0=xt(:,kn);% 使用 ode45 数值求解器，逐步积分 Lorenz 微分方程[tk,xk]=ode45(@lorenz3_dxdt,[tt(kn),tt(k)],x0);% xt: 保存真实状态向量 x 随时间变化的结果xt(:,k)=xk(end,:);
end%%  observation 生成观测数据
% 通过向真实值 xt 添加高斯噪声模拟观测
% epso(n): 第 n 个变量的观测误差标准差
% xo: 为含噪声的观测数据xo = randn(N,KTt);
for n=1:Nxo(n,:) = epso(n)*xo(n,:) + xt(n,:);
end%% etkf loop 
%- setup 初始化 ETKF 同化过程
%-- parameters
% rho=rhot; beta=betat; sigma=sigmat;% Ken: 集合成员数
Ken=Kda;     %Kda=3%-- precomputed LETKF covariance
% 预设的协方差逆（用于简化计算）
Pbinv=eye(Ken)*(Ken-1);        %Ken=3%-- obs operator and inv(R)
% 观测算子，从状态空间映射到观测空间
H=zeros(L,N);  % 观测误差协方差矩阵的逆
Rinv=zeros(L);     %N=3(numbers of variable), L=2(number of observation)
for l=1:Ln = idy2x(l); H(l,n) = 1; Rinv(l,l) = 1/(epsoda(n)^2); 
end%-- timing
kdT_t2da = max([1,round(dTda/dTt)]);  
dTda = kdT_t2da*dTt;
tda = [0:dTda:Tend]; 
KTda = size(tda,2);
kt2da = [1+kdT_t2da:kdT_t2da:KTt];
KTba = 2*KTda-1;%-- ic: perturbed around the mean
% 初始化集合成员：在真实初始值 x0t 附近加扰动，生成初始集合成员
x0da=zeros(N,Ken);  
for n=1:N% 第 n 个变量的 Ken 个集合初始值x0da(n,:) = randn(1,Ken)*epsda0(n)+x0t(n);
end
% 存储整个集合随时间变化的状态
xda=zeros(N*Ken,KTba);  
k=1; 
xda(:,k)=x0da(:); %-- counter setup
kda=1;   
ko=0;%-- ETKF ETKF 同化循环for k=2:KTdakn=k-1; %- ensemble forecast 集合预报kdan = kda; kda = kda+1;x0da = xda(:,kdan);[tk,xk] = ode45(@lorenz3_dxdt,[tda(kn),tda(k)],x0da);            % 使用 ode45 对每个集合成员集体积分%- store backgroundxda(:,kda)=xk(end,:);
%   [xda(:,kda) xda(:,kdan)]%- ensemble analysis by ETKFko=ko+1;if(abs(tda(k)-tt(kt2da(ko)))>0.1*dTt); %  check timingfprintf('\n   inconsistency in timing (to,tda)=(%s,%s)',...num2str(tt(kt2da(ko))),num2str(tda(k)))else%%- prep%- actual obsyok=H*xo(:,kt2da(ko));%- ensemble mean and spread in xkdan=kda; kda=kda+1; xenbk=reshape(xda(:,kdan),N,Ken);xmnbk=mean(xenbk,2); Xbk=(xenbk-xmnbk*ones(1,Ken));%- ensemble mean and spread in yymnbk=H*xmnbk; Ybk=H*Xbk;%- innovation 生成观测预测（H*X）和实际观测对比（Innovation）dmnobk=yok-ymnbk;%%- LETKF%- inverse by SVD 使用 SVD 解协方差[U,S,V]=svd(Pbinv/rhoda+Ybk'*Rinv*Ybk);Pa=V*inv(S)*V';%- wa 计算分析增益 wawa=Pa*Ybk'*Rinv*dmnobk;%- WaWa=sqrtm((Ken-1)*Pa);%- analysis 更新集合 xenakxenak=(xmnbk+Xbk*wa)*ones(1,Ken)+Xbk*Wa;%%- store analysisxda(:,kda)=xenak(:);
%        [xda(:,kda) xda(:,kdan)], pauseend 
end%% plot 结果分析与可视化%- timing of a and b 
% total number of output for b and a
kda2b = [1 2:2:KTba];   
kda2a = [3:2:KTba];   %- xda mean
% 计算集合平均作为分析结果
xmn = zeros(N,KTba);
for k=1:KTbaxen = reshape(xda(:,k),N,Ken);xmn(:,k) = mean(xen,2);
end%- yo 观测结果
yo = H*xo(:,kt2da);figureUnits = 'centimeters';
figureWidth = 18; 
figureHeight = 12;%% 绘图1：（x, y, z）三维轨迹演化图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, :);                 % 每个点的颜色close all;
% 绘图
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',sigmat , rhot, betat));
view(30, 20);
axis tight;
set(gca, 'FontSize', 12, 'FontName', 'Times New Roman');% 保存图像
str = strcat(pathFigure, "Fig.1 EKTF 3D 轨迹图", '.tiff');
print(gcf, '-dtiff', '-r600', str);%% 绘图2：同化结果图
%- plot
figureHandle = figure('Name','x^t');
set(gcf, 'Units', figureUnits, 'Position', [0 0 figureWidth figureHeight]);for n=1:Nsubplot(N,1,n); hold on;box on;% 真值 Trueh(1) = plot(tt, xt(n,:), [clrt lnt], 'LineWidth',1.2);                        % 背景值/第一猜/先验 Backgroundh(2) = plot(tda,           xmn(n,kda2b),[clrb styb], 'MarkerSize', 6);        % 分析值h(3) = plot(tda(2:end),xmn(n,kda2a),[clra stya], 'MarkerSize', 6);    set(gca, 'FontSize', 12, 'FontName', 'Times New Roman');ylabel(['x_' int2str(n)],'FontSize',14,'FontName','Times New Roman')if(n==N); xlabel('time','FontSize',14,'FontName','Times New Roman'); end
endfor l=1:Lsubplot(N,1,idy2x(l)); % 观测值 Observationh(4) = plot( tt(kt2da),yo(l,:),[clro styo], 'MarkerSize', 6);          
endsubplot(N,1,1); 
hl = legend(h([1  2  3  4 ]),"True", "Background", "Analyze","Observation");
set(hl, 'NumColumns', 4, 'box', 'off', 'FontName', 'Times New Roman', 'FontSize', 12,'Location','northoutside');% 保存图像
str = strcat(pathFigure, "Fig.1 Results", '.tiff');
print(gcf, '-dtiff', '-r600', str);

---

ETKF配置代码如下：

%% etkf_v1_ip
% 输入配置脚本，定义初始条件、参数（如 rho, beta, sigma）、扰动、观测误差等%% input for etkf_v1%% definition
N = 3;                 %  number of variable  状态变量数量（Lorenz-63 模型有 3 个变量）
Kda=3;               % ensemble size  in da 数据同化用的集合成员数量（ensemble size）%% model parameters%- true
% 真正的物理系统参数（用于生成真值轨道）
rhot = 28; 
betat = 8/3; 
sigmat = 10;%- da system (can be different from DA)
% 同化系统使用的参数（可以与真实值不同）
rhoda = rhot; 
betada = betat; 
sigmada = sigmat;%% Time 时间设置Tend = 20;      % 模拟总时长
dTt = 0.025;    % xt output  模型积分输出时间间隔（真值轨道）
dTda = 0.5;     % da window 同化周期（观测间隔）%% Observation 观测设置epso=[0.1;0.1;0.1];  % noise std in yo 观测误差标准差% selective obs
L=2;                 % number of obs 实际观测到的变量数量
idy2x=[1 3];     % 被观测变量的索引（观测 x1 和 x3）%% DA system 同化系统配置%--initial perturnation
epsda0=[1;1;1];         % 初始扰动%-- obs
epsoda=epso;           % 同化系统中的观测误差%-- inflation
rhoda=2;                   % 通胀因子（用于解决样本方差不足）%% graphical input 图形绘制设置
% 设置不同类型数据的绘图颜色与样式，用于后续可视化。% clrt='m';       
clrt=[211,47,47]/255;     
lnt='-';clro='g';        
styo='x';clrb='r';  
styb='*';clra='b'; 
stya='o';%% spinup to get IC on manifold
% Spin-up 初始条件（使系统稳定在吸引子上）x0sut = [1;1;1];   % 初始条件
Tsu = 100;          % spin-up 时间
x0 = x0sut';
Ken = 1;             % 仅用于 spin-up 的单个轨道
rho=rhot; beta=betat; sigma=sigmat;% 为了确保初始条件位于 Lorenz 系统的吸引子上，通过积分一段时间（Tsu=100）进行"spin-up"
[tsu,xsu] = ode45(@lorenz3_dxdt,[0,Tsu],x0);% 最终状态 x0t 作为后续模拟的初始状态
x0t = xsu(end,:);

---

Lorenz 3 模型定义如下：

%% Lorenz 3 模型定义
function [dxdt]=lorenz3_dxdt(t , x)
%%
% rho, beta, sigma: Lorenz 系统的参数
% N: 系统状态变量维度（Lorenz 3 系统中为 3）
% Ken: 集合成员数（ensemble size）global rho beta sigma
global N Kendxdt = zeros(N, Ken);
xen=reshape(x,N,Ken);dxdt(1,:)=sigma*(xen(2,:)-xen(1,:));
dxdt(2,:)=xen(1,:).*(rho-xen(3,:))-xen(2,:);
dxdt(3,:)=xen(1,:).*xen(2,:)-beta*xen(3,:);dxdt=dxdt(:);return

参考