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

sigma = 10.0; beta = 8.0/3.0; rho = 28.0        # 参数值   
dt = 0.01                                       # 模式步长

x0 = np.array([1.508870, -1.531271, 25.46091])  # 我们也用洛伦茨开展他的实验用的初值，这里的矢量x0代表了(x,y,z)的三个应变量，写成矢量主要为方便编程
Xt = np.zeros([3000,3]);Xt[0]=x0    # 我们积分3000步看看结果，这里涉及python的数组保存方式，第一维代表时间步数，第二维代表x,y,z的3个变量
                                    # 二维数值Xt使用Xt[0]的方式调用时，相当于Xt[0,:]，意味着第一维取0，第二维用全体，所以时间维放前面代码简单
for k in range(2999):             # 从0积分到2999
    Xt[k+1] = RK4(Lorenz63,Xt[k],dt,sigma,beta,rho) 

#%% 画图    
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(Xt[:,0], Xt[:,1], Xt[:,2],'r', label='Lorenz 63 model')
ax.legend()
plt.xlabel('x');plt.ylabel('y');
ax.set_zlabel('z')

Text(0.5, 0, 'z')

Text(0.5, 0, 'z')

# 模式积分
x0p = np.array([1.508, -1.531, 25.460])           # 这一次，我们只使用3位小数
Xp = np.zeros([3000,3]);Xp[0]=x0p
for k in range(2999):                             # 从0积分到2999
    Xp[k+1] = RK4(Lorenz63,Xp[k],dt,sigma,beta,rho) 
# 这一次，我们使用完全相同的流程积分模式，只不过初值忽略了小于千分之一的量，即初始误差小于千分之一。
#%% 画图
plt.figure(figsize=(8,4))
Ylabels = ['x','y','z']
for j in range(3):
    plt.subplot(3,1,j+1)
    plt.plot(dt*np.arange(0,3000), Xt[:,j],'k:', label='original')
    plt.plot(dt*np.arange(0,3000), Xp[:,j],'b',label='restart')
    plt.xlim(0,30);plt.ylabel(Ylabels[j],fontsize=14)
    plt.grid(axis='x')
    if j<2:
        plt.xticks(dt*np.arange(0,3001,300),[])
    if j==2:
        plt.legend(ncol=2);plt.xlabel('Time Units')
        plt.xticks(dt*np.arange(0,3001,300))

#%% 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.50                                       # 观测误差标准差
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                                                   # 观测计数，用于循环控制
# 画图
## 以下提供真值和观测画图的参考脚本
import matplotlib.pyplot as plt
plt.figure(figsize=(10,6))
lbs = ['x','y','z']
for j in range(3):
    plt.subplot(3,1,j+1)
    plt.plot(t,Xtrue[j],'b-',lw=2,label='True')
    plt.plot(t_m,yo[j],'go',ms=8,markerfacecolor='white',label='Observation')
    plt.ylabel(lbs[j],fontsize=16)
    plt.xticks(fontsize=16);plt.yticks(fontsize=16)
    if j==0:
        plt.legend(ncol=2, fontsize=12)
        plt.ylim(-21,21)
        plt.title('The True States and Observations',fontsize=16)
    if j==2:
        plt.xlabel('Time',fontsize=16)

def Lin3dvar(xb,w,H,R,B,opt):
# 三维变分的算法
    if opt == 1:                        # 模式空间算法
        Bi = np.linalg.inv(B)
        Ri = np.linalg.inv(R)
        A = Bi + (H.T)@Ri@H
        b = Bi@xb + (H.T)@Ri@w
        xa = np.linalg.solve(A,b)       # 求解线性方程组
    elif opt == 2:                      # 模式空间增量算法
        Bi = np.linalg.inv(B)
        Ri = np.linalg.inv(R)
        A = Bi + (H.T)@Ri@H
        b = (H.T)@Ri@(w-H@xb)
        xa = xb + np.linalg.solve(A,b)  # 求解线性方程组
    elif opt == 3:                      # 观测空间增量方法
        A = R + H@B@(H.T)
        b = (w-H@xb)
        xa = xb + B@(H.T)@np.linalg.solve(A,b) #solve a linear system
    return xa

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                                 # 设定初始的背景误差
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,1)  # 调用3dvar，更新状态和协方差
        km = km+1

# 3DVar结果画图
plt.figure(figsize=(10,8))
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(fontsize=16);plt.yticks(fontsize=16)
    if j==0:
        plt.legend(ncol=2, loc='upper right',fontsize=12)
        plt.title("3DVar",fontsize=16)
RMSEb = np.sqrt(np.mean((Xb-Xtrue)**2,0))
RMSEa = np.sqrt(np.mean((Xa-Xtrue)**2,0))
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 0x10c528ed0>

<matplotlib.legend.Legend at 0x10c528ed0>

用Lorenz1963模式展示的数据同化实验¶

关于实验和模型的一些说明¶

关于编程环境配置¶

关于数值模式¶

Lorenz63模式的数值解¶

解对于初值的依赖性¶

同化实验¶

讨论¶