# 定义模式方程和积分格式
import numpy as np
def Lorenz63(state,*args):       #Lorenz63模式
    sigma = args[0]
    beta = args[1]
    rho = args[2]
    x, y, z = state 
    f = np.zeros(3) 
    f[0] = sigma * (y - x)
    f[1] = x * (rho - z) - y
    f[2] = x * y - beta * z
    return f 
def RK4(rhs,state,dt,*args):    # Runge-Kutta格式，输入的rhs代表模式右端方程
    k1 = rhs(state,*args)
    k2 = rhs(state+k1*dt/2,*args)
    k3 = rhs(state+k2*dt/2,*args)
    k4 = rhs(state+k3*dt,*args) # 输出新的一步的状态
    new_state = state + (dt/6)*(k1+2*k2+2*k3+k4)
    return new_state

# 之后使用 RK4(Lorenz63,x_in,dt,sigma,beta,rho)可以进行步长为dt的一步积分

# Lorenz63模式真值试验和观测构造
sigma = 10.0; beta = 8.0/3.0; rho = 28.0          # 模式参数值   
dt = 0.01                                         # 模式积分步长
n = 3                                             # 模式的状态维数，即有多少个状态变量
m = 3                                             # 模式的观测维数，即多少个变量可以被观测到，这里先假定所有变量都能被观测
tm = 10                                           # 同化试验时间窗口
nt = int(tm/dt)                                   # 总积分步数
t = np.linspace(0,tm,nt+1)                        # 模式时间网格
x0True = np.array([1.508870, -1.531271, 25.46091])# 真实值的初值
np.random.seed(seed=1)              # 设置随机种子，由于电脑的随机数是伪随机，记录了随机种子之后，每次运行这个脚本产生的“随机”的观测误差都是一样的。
sig_m= 0.5                                        # 观测误差标准差
R = sig_m**2*np.eye(n)                            # 观测误差协方差矩阵，设为对角阵使得不同变量的误差互不相干
dt_m = 0.2                                        # 观测之间的时间间隔（即每20模式步观测一次）
tm_m = 10                                         # 最大观测时间（即多少时间之后停止同化，可小于同化试验时间窗口）
nt_m = int(tm_m/dt_m)                             # 进行同化的总次数
ind_m = (np.linspace(int(dt_m/dt),int(tm_m/dt),nt_m)).astype(int)   # 这是有观测的时刻在总时间网格中的位置指标
t_m = t[ind_m]                                                      # 观测网格
def h(x):                                         # 定义观测算子：观测算子用于构建模式变量和观测之间的关系。当所有变量被观测时，使用单位阵。
    H = np.eye(n)                                 # 观测矩阵为单位阵。
    yo = H@x                                      # 单位阵乘以状态变量，即所有变量直接被观测。
    return yo
Xtrue = np.zeros([n,nt+1])                        # 真值保存在xTrue变量中
Xtrue[:,0] = x0True                               # 初始化真值
km = 0                                            # 观测计数
yo = np.zeros([3,nt_m])                           # 观测保存在yo变量中
for k in range(nt):                               # 按模式时间网格开展模式积分循环
    Xtrue[:,k+1] = RK4(Lorenz63,Xtrue[:,k],dt,sigma,beta,rho)       # 真实值积分
    if (km<nt_m) and (k+1==ind_m[km]):                              # 用指标判断是否进行观测
        yo[:,km] = h(Xtrue[:,k+1]) + np.random.normal(0,sig_m,[3,]) # 通过判断，在观测时间取出真值作为观测值，同时叠加高斯分布随机噪声
        km = km+1                                                   # 观测计数，用于循环控制
        
# 这个脚本产生了同化实验中的观测数据yo和用于对比结果的真值Xtrue。

def Lin3dvar(xb,w,H,R,B):                  # 三维变分同化算法
    A = R + H@B@(H.T)
    b = (w-H@xb)
    xa = xb + B@(H.T)@np.linalg.solve(A,b) # 求解线性方程组
    return xa
# 由于我并不打算过于详细地介绍3DVar的各种算法细节，我只使用最简单的一个版本展开介绍

# 同化实验流程
x0b = np.array([1,-1,20])                # 背景积分的初值
Xb = np.zeros([3,nt+1]); Xb[:,0] = x0b   # 控制试验结果存在xb中
# --------------- 背景积分实验 ------------------------
for k in range(nt):                      # 模式积分循环
    Xb[:,k+1] = RK4(Lorenz63,Xb[:,k],dt,sigma,beta,rho)   # 不加同化的背景积分结果，后面和同化结果进行对比
# --------------- 数据同化实验 ------------------------
sig_b= 1.0                               # 设定初始的背景误差
B = sig_b**2*np.eye(3)                   # 设定初始背景误差协方差矩阵，B矩阵的取值对于变分同化比较重要，这个简单模式可以使用简单的对角阵
Xa = np.zeros([3,nt+1]); Xa[:,0] = x0b   # 同化试验结果存在Xa中，第一步没同化，所以数值也是x0b
km = 0                                   # 同化次数计数
H = np.eye(3)                            # 如前述，观测算子是单位阵
for k in range(nt):                      # 模式积分循环
    Xa[:,k+1] = RK4(Lorenz63,Xa[:,k],dt,sigma,beta,rho)   # 没有遇到观测的时候，就是正常积分模式
    if (km<nt_m) and (k+1==ind_m[km]):   # 当有观测时，使用3dvar同化
        Xa[:,k+1] = Lin3dvar(Xa[:,k+1],yo[:,km],H,R,B)  # 调用3dvar，更新状态和协方差
        km = km+1

RMSEb = np.sqrt(np.mean((Xb-Xtrue)**2,0))
RMSEa = np.sqrt(np.mean((Xa-Xtrue)**2,0))
# 结果画图
import matplotlib.pyplot as plt
plt.figure(figsize=(8,6))
lbs = ['x','y','z']
for j in range(3):
    plt.subplot(4,1,j+1)
    plt.plot(t,Xtrue[j],'b-',lw=2,label='True')
    plt.plot(t,Xb[j],'--',color='orange',lw=2,label='Background')
    plt.plot(t_m,yo[j],'go',ms=8,markerfacecolor='white',label='Observation')
    plt.plot(t,Xa[j],'-.',color='red',lw=2,label='Analysis')
    plt.ylabel(lbs[j],fontsize=16)
    plt.xticks(np.arange(0,10.1,2),[]);plt.yticks(fontsize=16)
    if j==0:
        plt.legend(ncol=2, loc='upper right',fontsize=12)
        plt.title("3DVar",fontsize=16)
plt.subplot(4,1,4)
plt.plot(t,RMSEb,'--',color='orange',label='Background')
plt.plot(t,RMSEa,color='red',label='Analysis')
plt.text(0,20,'mean RMSEb = %0.2f' %np.mean(RMSEb))
plt.text(0,25,'mean RMSEa = %0.2f' %np.mean(RMSEa))
plt.ylabel('RMSE',fontsize=16)
plt.xlabel('Time',fontsize=16)
plt.xticks(fontsize=16);plt.yticks(fontsize=16)
plt.legend(ncol=2, loc='upper right',fontsize=12)

<matplotlib.legend.Legend at 0x1067d8510>

#%% Lorenz63模式的切线性模式
def JLorenz63(state,*args):         # Lorenz63 方程雅克比矩阵
    sigma = args[0]
    beta = args[1]
    rho = args[2]
    x, y, z = state
    
    df = np.zeros([3,3])           # 以下是切线性矩阵的9个元素配置
    df[0,0] = sigma * (-1)
    df[0,1] = sigma * (1)
    df[0,2] = sigma * (0)
    df[1,0] = 1 * (rho - z) 
    df[1,1] = -1
    df[1,2] = x * (-1)
    df[2,0] = 1 * y 
    df[2,1] = x * 1 
    df[2,2] = - beta
    return df 

def JRK4(rhs,Jrhs,state,dt,*args):   # 切线性模式的积分格式    
    n = len(state)
    k1 = rhs(state,*args)
    k2 = rhs(state+k1*dt/2,*args)
    k3 = rhs(state+k2*dt/2,*args)
    # 以下是对矩阵的Runge-Kutta格式
    dk1 = Jrhs(state,*args)
    dk2 = Jrhs(state+k1*dt/2,*args) @ (np.eye(n)+dk1*dt/2) 
    dk3 = Jrhs(state+k2*dt/2,*args) @ (np.eye(n)+dk2*dt/2) 
    dk4 = Jrhs(state+k3*dt,*args) @ (np.eye(n)+dk3*dt)  
    DM = np.eye(n) + (dt/6) * (dk1+2*dk2+2*dk3+dk4)
    return DM

def EKF(xb,yo,ObsOp,Dh,R,B):
# 输入的变量分别为:xb 预报、yo 观测、ObsOp 观测算子、Dh 切线观测矩阵，R 观测误差协方差， # B 背景误差协方差。
    n = xb.shape[0]                      # 状态空间维数
    D = Dh@B@Dh.T + R
    K = B @ Dh.T @ np.linalg.inv(D)      # 卡尔曼增益矩阵
    xa = xb + K @ (yo-ObsOp(xb))         # 更新状态
    P = (np.eye(n) - K@Dh) @ B           # 更新误差协方差矩阵
    return xa, P                         # 输出分析状态场和分析误差协方差矩阵

# 同化实验流程
x0b = np.array([1,-1,20])                # 背景积分的初值
Xb = np.zeros([3,nt+1]); Xb[:,0] = x0b   # 控制试验结果存在xb中
# --------------- 背景积分实验 ------------------------
for k in range(nt):                      # 模式积分循环
    Xb[:,k+1] = RK4(Lorenz63,Xb[:,k],dt,sigma,beta,rho)   # 不加同化的背景积分结果，后面和同化结果进行对比
# --------------- 数据同化实验 ------------------------
sig_b= 1.0                               # 设定初始的背景误差
Q = 0.0*np.eye(3)                        # !!!!!!!!!!!!!!!!!!!设置模式误差（若假设完美模式则取0）
B = sig_b**2*np.eye(3)                   # 设定初始背景误差协方差矩阵，B矩阵的取值对于变分同化比较重要，这个简单模式可以使用简单的对角阵
Xa = np.zeros([3,nt+1]); Xa[:,0] = x0b   # 同化试验结果存在Xa中，第一步没同化，所以数值也是x0b
km = 0                                   # 同化次数计数
H = np.eye(3)                            # 如前述，观测算子是单位阵
for k in range(nt):                      # 模式积分循环
    Xa[:,k+1] = RK4(Lorenz63,Xa[:,k],dt,sigma,beta,rho)      # 没有遇到观测的时候，就是正常积分模式
    DM = JRK4(Lorenz63,JLorenz63,Xa[:,k],dt,sigma,beta,rho)  # !!!!!!!!!!!!!!!!!!!!!!!!使用切线性模式积分
    B = DM @ B @ DM.T + Q                                    # !!!!!!!!!!!!!!!!!!!!!!!!积分过程协方差更新
    if (km<nt_m) and (k+1==ind_m[km]):                       # 当有观测时
        Xa[:,k+1],B = EKF(Xa[:,k+1],yo[:,km],h,H,R,B)        # !!!!!!!!!!!!!!!!!!调用EKF，更新状态和协方差
        km = km+1

RMSEb = np.sqrt(np.mean((Xb-Xtrue)**2,0))
RMSEa = np.sqrt(np.mean((Xa-Xtrue)**2,0))
# 结果画图
import matplotlib.pyplot as plt
plt.figure(figsize=(8,6))
lbs = ['x','y','z']
for j in range(3):
    plt.subplot(4,1,j+1)
    plt.plot(t,Xtrue[j],'b-',lw=2,label='True')
    plt.plot(t,Xb[j],'--',color='orange',lw=2,label='Background')
    plt.plot(t_m,yo[j],'go',ms=8,markerfacecolor='white',label='Observation')
    plt.plot(t,Xa[j],'-.',color='red',lw=2,label='Analysis')
    plt.ylabel(lbs[j],fontsize=16)
    plt.xticks(np.arange(0,10.1,2),[]);plt.yticks(fontsize=16)
    if j==0:
        plt.legend(ncol=2, loc='upper right',fontsize=12)
        plt.title("EKF",fontsize=16)     #!!!!!!!!!!!!!!!!!!!!!!!!!!!
plt.subplot(4,1,4)
plt.plot(t,RMSEb,'--',color='orange',label='Background')
plt.plot(t,RMSEa,color='red',label='Analysis')
plt.text(0,20,'mean RMSEb = %0.2f' %np.mean(RMSEb))
plt.text(0,25,'mean RMSEa = %0.2f' %np.mean(RMSEa))
plt.ylabel('RMSE',fontsize=16)
plt.xlabel('Time',fontsize=16)
plt.xticks(fontsize=16);plt.yticks(fontsize=16)
plt.legend(ncol=2, loc='upper right',fontsize=12)

<matplotlib.legend.Legend at 0x1066d0510>

plt.figure(figsize=(5,5))
plt.imshow(B/np.max(B),cmap='Reds')
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x106afc0d0>

def EnKF(xbi,yo,ObsOp,Dh,R):
# 输入变量依次为预报集合、观测、观测算子、切线观测矩阵、观测误差协方差
    n,N = xbi.shape                        # n-状态维数，N-集合成员数
    m = yo.shape[0]                        # m-观测维数
    xb = np.mean(xbi,1)                    # 计算预报集合平均                      
    B = (1/(N-1)) * (xbi - xb.reshape(-1,1)) @ (xbi - xb.reshape(-1,1)).T

    #背景误差协方差
    D = Dh@B@Dh.T + R
    K = B @ Dh.T @ np.linalg.inv(D) #计算卡尔曼增益矩阵
    yoi = np.zeros([m,N])                                                # 预分配空间，保存扰动后的观测集合 xai = np.zeros([n,N]) 
    for i in range(N):
        yoi[:,i] = yo + np.random.multivariate_normal(np.zeros(m), R)        # 随机扰动观测
        xai[:,i] = xbi[:,i] + K @ (yoi[:,i]-ObsOp(xbi[:,i]))                 # 更新每个成员
    return xai                            # 输入和输出都是一整个集合。不同于 EKF，不需要输出分析误差协方差

# EnKF同化实验流程
x0b = np.array([1,-1,20])                # 背景积分的初值
Xb = np.zeros([3,nt+1]); Xb[:,0] = x0b   # 控制试验结果存在xb中
# --------------- 背景积分实验 ------------------------
for k in range(nt):                      # 模式积分循环
    Xb[:,k+1] = RK4(Lorenz63,Xb[:,k],dt,sigma,beta,rho)   # 不加同化的背景积分结果，后面和同化结果进行对比
# --------------- 数据同化实验 ------------------------
sig_b= 1.0                               # 设定初始的背景误差
B = sig_b**2*np.eye(3)                   # 设定初始背景误差协方差矩阵，B矩阵的取值对于变分同化比较重要，这个简单模式可以使用简单的对角阵
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
N = 20                                   # 设定集合成员数
xai = np.zeros([3,N])                    # 设定集合，保存在 xai 中
for i in range(N):
    xai[:,i] = x0b + np.random.multivariate_normal(np.zeros(n), B)    # 通过对预报初值进行随机扰动构造初始集合
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Xa = np.zeros([3,nt+1]); Xa[:,0] = x0b   # 同化试验结果存在Xa中，第一步没同化，所以数值也是x0b
km = 0                                   # 同化次数计数
H = np.eye(3)                            # 如前述，观测算子是单位阵
for k in range(nt):                      # 模式积分循环
    for i in range(N):                   # !!!!!!!!!!!!!!!!!!!!!!!!!对每个集合成员积分
        xai[:,i] = RK4(Lorenz63,xai[:,i],dt,sigma,beta,rho)     #!!!!!!!!!!!!!!!!!
    if (km<nt_m) and (k+1==ind_m[km]):   # 当有观测时，使用3dvar同化
        xai = EnKF(xai,yo[:,km],h,H,R)  # !!!!!!!!!!!!!!!!!!!!!!!!调用 EnKF 同化
        km = km+1
    Xa[:,k+1] = np.mean(xai,1)           # !!!!!!!!!!!!!!!!!!!!!!!!非同化时刻求预报平均，同化时刻求分析平均

RMSEb = np.sqrt(np.mean((Xb-Xtrue)**2,0))
RMSEa = np.sqrt(np.mean((Xa-Xtrue)**2,0))
# 结果画图
import matplotlib.pyplot as plt
plt.figure(figsize=(8,6))
lbs = ['x','y','z']
for j in range(3):
    plt.subplot(4,1,j+1)
    plt.plot(t,Xtrue[j],'b-',lw=2,label='True')
    plt.plot(t,Xb[j],'--',color='orange',lw=2,label='Background')
    plt.plot(t_m,yo[j],'go',ms=8,markerfacecolor='white',label='Observation')
    plt.plot(t,Xa[j],'-.',color='red',lw=2,label='Analysis')
    plt.ylabel(lbs[j],fontsize=16)
    plt.xticks(np.arange(0,10.1,2),[]);plt.yticks(fontsize=16)
    if j==0:
        plt.legend(ncol=2, loc='upper right',fontsize=12)
        plt.title("EnKF",fontsize=16)     #!!!!!!!!!!!!!!!!!!!!!!!!!!!
plt.subplot(4,1,4)
plt.plot(t,RMSEb,'--',color='orange',label='Background')
plt.plot(t,RMSEa,color='red',label='Analysis')
plt.text(0,20,'mean RMSEb = %0.2f' %np.mean(RMSEb))
plt.text(0,25,'mean RMSEa = %0.2f' %np.mean(RMSEa))
plt.ylabel('RMSE',fontsize=16)
plt.xlabel('Time',fontsize=16)
plt.xticks(fontsize=16);plt.yticks(fontsize=16)
plt.legend(ncol=2, loc='upper right',fontsize=12)

<matplotlib.legend.Legend at 0x106ca71d0>

#%% 集合平方根滤波器（直接法）
def EnSRF(xbi,yo,ObsOp,Dh,R):
    from scipy.linalg import sqrtm
    n,N = xbi.shape                      # n-状态维数，N-集合成员数
    m = yo.shape[0]                      # m-观测维数
    xb = np.mean(xbi,1)                   # 预报集合平均 
    B = (1/(N-1)) * (xbi - xb.reshape(-1,1)) @ (xbi - xb.reshape(-1,1)).T
    D = Dh@B@Dh.T + R
    K = B @ Dh.T @ np.linalg.inv(D)         # !!!以上与EnKF一致
    xa = xb + K @ (yo-ObsOp(xb))          # 用确定性格式更新集合平均

    A = xbi - xb.reshape(-1,1)              # 集合异常
    Z = A/np.sqrt(N-1)                    # 标准化集合异常值
    Y = np.linalg.inv(D)@Dh@Z           
    X = sqrtm(np.eye(N)-(Dh@Z).T@Y)      # 矩阵平方根
    X = np.real(X)                         # 保证矩阵平方根为实数
    
    Z = Z@X                             # 更新集合异常值
    A = Z*np.sqrt(N-1)
    xai = xa.reshape(-1,1)+A               # 用集合平均和集合异常计算集合成员
    return xai

#%% 集合平方根卡尔曼滤波器(串行格式)
def sEnSRF(xbi,yo,ObsOp,Dh,R):
    n,N = xbi.shape                     # n-状态维数，N-集合成员数
    m = yo.shape[0]                     # m-观测维数
    xb = np.mean(xbi,1)                 # 预报集合平均 
    B = (1/(N-1)) * (xbi - xb.reshape(-1,1)) @ (xbi - xb.reshape(-1,1)).T
    D = Dh@B@Dh.T + R
    K = B @ Dh.T @ np.linalg.inv(D)       # !!!以上与EnKF一致
    xa = xb + K @ (yo-ObsOp(xb))        # 用确定性格式更新集合平均

    A = xbi - xb.reshape(-1,1)            # 集合异常
    Z = A/np.sqrt(N-1)                  # 标准化集合异常值
    V = (Dh@Z).T
    for j in range(m):                       # 根据每个观测循环
        Dj = V[:,j].T @ V[:,j] + R[j,j] 
        betaj = 1/(Dj+np.sqrt(R[j,j]*Dj))
        Z = Z@(np.eye(N)-betaj*V[:,j]@V[:,j].T)     # 集合异常更新公式
    A = Z*np.sqrt(N-1)
    xai = xa.reshape(-1,1)+A               # 用集合平均和集合异常计算集合成员
    return xai

#%% 集合转换卡尔曼滤波器
def ETKF(xbi,yo,ObsOp,Dh,R):
    n,N = xbi.shape                      # n-状态维数，N-集合成员数
    m = yo.shape[0]                     # m-观测维数
    xb = np.mean(xbi,1)                  # 预报集合平均 
    B = (1/(N-1)) * (xbi - xb.reshape(-1,1)) @ (xbi - xb.reshape(-1,1)).T
    D = Dh@B@Dh.T + R
    K = B @ Dh.T @ np.linalg.inv(D)        # !!!以上与EnKF一致
    xa = xb + K @ (yo-ObsOp(xb))         # 用确定性格式更新集合平均

    A = xbi - xb.reshape(-1,1)             # 集合异常
    Z = A/np.sqrt(N-1)                   # 标准化集合异常值
    V = (Dh@Z).T
    CTC = V@np.linalg.inv(R)@V.T
    Gamma, C = np.linalg.eig(CTC)
    Gamma = np.real(Gamma);C = np.real(C)    
    Z = Z@C@np.diag((Gamma+1)**(-0.5))  # 集合异常更新公式
    A = Z*np.sqrt(N-1)
    xai = xa.reshape(-1,1)+A              # 用集合平均和集合异常计算集合成员
    return xai

#%% 集合调整卡尔曼滤波器
def EAKF(xbi,yo,ObsOp,Dh,R):
    n,N = xbi.shape                       # n-状态维数，N-集合成员数
    m = yo.shape[0]                      # m-观测维数
    xb = np.mean(xbi,1)                   # 预报集合平均 
    B = (1/(N-1)) * (xbi - xb.reshape(-1,1)) @ (xbi - xb.reshape(-1,1)).T
    D = Dh@B@Dh.T + R
    K = B @ Dh.T @ np.linalg.inv(D)          # !!!以上与EnKF一致
    xa = xb + K @ (yo-ObsOp(xb))           # 用确定性格式更新集合平均

    A = xbi - xb.reshape(-1,1)               # 集合异常
    Z = A/np.sqrt(N-1)                         # 标准化集合异常值
    V = (Dh@Z).T
    CTC = V@np.linalg.inv(R)@V.T
    Gamma, C = np.linalg.eig(CTC)
    Gamma = np.real(Gamma);C = np.real(C)   
    F,G,U = np.linalg.svd(Z)
    IG2 = np.diag((Gamma+1)**(-0.5))
    Gtilde = np.concatenate([np.diag(1/G),np.zeros([N-m,m])],0)   # !!!公式（6-1-29）
    Adj = Z@C@IG2@Gtilde@F.T                             # !!!公式（6-1-30）
    Z = Adj@Z                                              # 公式（6-1-27）                    
    A = Z*np.sqrt(N-1)
    xai = xa.reshape(-1,1)+A            # 用集合平均和集合异常计算集合成员
    return xai

#%% 确定性集合卡尔曼滤波器
def DEnKF(xbi,yo,ObsOp,Dh,R):
    n,N = xbi.shape         # n-状态维数，N-集合成员数
    m = yo.shape[0]         # m-观测维数
    xb = np.mean(xbi,1)     # 预报集合平均 
    B = (1/(N-1)) * (xbi - xb.reshape(-1,1)) @ (xbi - xb.reshape(-1,1)).T
    D = Dh@B@Dh.T + R
    K = B @ Dh.T @ np.linalg.inv(D)     # !!!以上与EnKF一致
    xa = xb + K @ (yo-ObsOp(xb))        # 更新集合成员
    
    A = xbi - xb.reshape(-1,1)          # 集合异常
    A = A - (1/2) * K @ Dh @ A          # !! DEnKF关键步骤：补偿离散度的操作
    xai = xa.reshape(-1,1) + A
          
    return xai

#%% （串行）集合调整卡尔曼滤波器
def obs_increment_eakf(ensemble, observation, obs_error_var): # 1：计算观测空间增量
    prior_mean = np.mean(ensemble);
    prior_var = np.var(ensemble);
    if prior_var >1e-6:                                                 # 用于避免退化的先验集合造成错误更新
        post_var = 1.0 / (1.0 / prior_var + 1.0 / obs_error_var);       # 公式（2）
        post_mean = post_var * (prior_mean / prior_var + observation / obs_error_var);  # 公式（3）
    else:
        post_var = prior_var; post_mean = prior_mean;    
    updated_ensemble = ensemble - prior_mean + post_mean;
    var_ratio = post_var / prior_var;
    updated_ensemble = np.sqrt(var_ratio) * (updated_ensemble - post_mean) + post_mean; # 公式（4）
    obs_increments = updated_ensemble - ensemble;
    return obs_increments
def get_state_increments(state_ens, obs_ens, obs_incs): # 2将观测增量回归到状态增量
    covar = np.cov(state_ens, obs_ens);
    state_incs = obs_incs * covar[0,1]/covar[1,1];
    return state_incs
def sEAKF(xai,yo,ObsOp, R, RhoM):
    n,N = xai.shape;                  # 状态维数
    m = yo.shape[0];                  # 观测数
    Loc = ObsOp(RhoM)                 # 观测空间局地化
    for i in range(m):                # 针对每个标量观测的循环
        hx = ObsOp(xai);              # 投影到观测空间 
        hxi = hx[i];                  # 投影到对应的矢量观测
        obs_inc = obs_increment_eakf(hxi,yo[i],R[i,i]);
        for j in range(n):            # 针对状态变量的每个元素的循环
            state_inc = get_state_increments(xai[j], hxi,obs_inc)   # 获取状态增量
            cov_factor=Loc[i,j]       # 使用局地化矩阵的相应元素
            if cov_factor>1e-6:       # 在局地化范围内加增量
                xai[j]=xai[j]+cov_factor*state_inc; 
    return xai

数据同化的算法演变：3DVar、EKF、EnKF¶

关于同化流程的代码回顾¶

三维变分¶

卡尔曼滤波器和扩展卡尔曼滤波器¶

集合卡尔曼滤波器¶

EnKF的衍生算法¶

结语¶