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

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

模型概述

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 all

pathFigure = "..\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
%-- parameters

rho=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:N
    xo(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:L
    n = 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:KTda
    kn=k-1; 

    %- ensemble forecast 集合预报
    kdan = kda; 
    kda = kda+1;
    x0da = xda(:,kdan);
    [tk,xk] = ode45(@lorenz3_dxdt,[tda(kn),tda(k)],x0da);            % 使用 ode45 对每个集合成员集体积分

    %- store background
    xda(:,kda)=xk(end,:);
%   [xda(:,kda) xda(:,kdan)]
    %- ensemble analysis by ETKF
    ko=ko+1;
    if(abs(tda(k)-tt(kt2da(ko)))>0.1*dTt); %  check timing
        fprintf('\n   inconsistency in timing (to,tda)=(%s,%s)',...
            num2str(tt(kt2da(ko))),num2str(tda(k)))
    else
        %%- prep
        %- actual obs
        yok=H*xo(:,kt2da(ko));
        %- ensemble mean and spread in x
        kdan=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 y
        ymnbk=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 计算分析增益 wa
        wa=Pa*Ybk'*Rinv*dmnobk;
        %- Wa
        Wa=sqrtm((Ken-1)*Pa);

        %- analysis 更新集合 xenak
        xenak=(xmnbk+Xbk*wa)*ones(1,Ken)+Xbk*Wa;
        %%- store analysis
        xda(:,kda)=xenak(:);
%        [xda(:,kda) xda(:,kdan)], pause
    end 
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:KTba
    xen = 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:N
    subplot(N,1,n); 
    hold on;box on;

    % 真值 True
    h(1) = plot(tt, xt(n,:), [clrt lnt], 'LineWidth',1.2);                        
    
    % 背景值/第一猜/先验 Background
    h(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
end

for l=1:L
    subplot(N,1,idy2x(l)); 
    % 观测值 Observation
    h(4) = plot( tt(kt2da),yo(l,:),[clro styo], 'MarkerSize', 6);          
end

subplot(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 Ken

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