【技术实现步骤摘要】
一种多变量频域同化方法及系统
[0001]本专利技术涉及一种多变量数据同化方法及系统，尤其涉及一种多变量频域同化方法及系统。
技术介绍
[0002]陆面过程模型依靠其内在的物理过程和动力学机制，可以对地表
‑
大气界面的地表温度、土壤湿度、土壤温度等与水、能循环相关的地表状态变量在时空上进行连续模拟，是研究大尺度和流域尺度水分和能量循环的重要方法，但陆面过程模型输入数据质量和模型参数化方案的非完备性使其模拟精度较低。
[0003]地面观测可以时间连续地对地表变量进行实测，但观测点位通常较少、空间代表性有限，无法实现地表状态变量大面积动态监测。陆面数据同化方法通过有机结合地面观测产品和陆面模型模拟各自的优势，能够有效提高陆面过程模型模拟地表状态变量的精度，为区域水、热通量准确连续模拟提供了新途径。
[0004]当前，国内外常用的同化方法主要有顺序同化和连续同化两类。顺序同化根据积分时刻观测产品和陆面过程模型模拟的不确定性，利用观测产品对陆面过程模型状态变量进行更新，从而依次获得陆面过程模型状态变量的最优后验估计。连续同化在同化窗口内构建最优状态变量与观测变量和初始状态变量之间偏差的目标函数，利用迭代优化算法求解最优状态变量。
[0005]顺序同化和连续同化都是在时间域内寻找陆面过程模型状态变量的最优解。但地表过程状态演化通常涉及到很多时间尺度的物理信息，空间均匀性较大的状态变量，例如风速，变化较慢，而空间均匀性较小的状态变量，例如地表温度，则变化较快，这导致不同地表状态变量时序信
【技术保护点】
【技术特征摘要】
1.一种多变量频域同化方法，其特征在于，包括以下步骤：步骤1：确定研究区域、同化窗口长度、校正向量及同化向量，准备陆面过程模型及其驱动数据、观测算子及其参数数据、同化向量的观测产品；步骤2：生成陆面过程模型的校正向量集合；步骤3：利用校正向量集合向前积分陆面过程模型和观测算子得到同化窗口内各积分时刻校正向量的预报集合和同化向量的预报集合；步骤4：利用校正向量的预报集合和同化向量的预报集合借助同化向量的观测产品在频率域实现校正向量的同化；步骤5：基于同化过的校正向量，按照步骤2
‑
步骤4执行后续时刻的同化。2.如权利要求1所述的一种多变量频域同化方法，其特征在于，所述步骤2利用校正向量对应同化时刻的历史均值向量和历史方差向量生成校正向量集合，具体步骤为：步骤2011：产生N
×
M个在(0，1)上均匀分布的随机变量u1
nm
和u2
nm
，1≤n≤N，1≤m≤M，其中，n，m表示随机变量的下标变量，N为校正向量的维数，M为校正向量集合的大小；步骤2012：利用N
×
M个u1
nm
和u2
nm
产生N
×
M个服从标准正态分布的随机变量U
nm
，公式如下：步骤2013：利用N
×
M个U
nm
产生一个大小为M的服从N维正态分布的随机校正向量集合{X
Nm
}
1≤m≤M
：X
Nm
＝[x
1m
,
…
,x
Nm
]'x
1m
＝U
1m
×
σ1+μ1,
…
,x
Nm
＝U
Nm
×
σ
N
+μ
N
式中，
′
表示向量或矩阵转置，[μ1,
…
,μ
N
]
′
和[σ1,
…
,σ
N
]
′
分别为校正向量X
N
＝[x1,
…
,x
n
,
…
,x
N
]
′
对应同化时刻的历史均值向量和历史方差向量。3.如权利要求1所述的一种多变量频域同化方法，其特征在于，所述步骤3的具体步骤为：步骤3011：将{X
Nm
}
1≤m≤M
每个成员X
Nm
输入陆面过程模型并在驱动数据驱动下向前积分生成同化窗口内各积分时刻校正向量的预报结果，积分公式为：式中，t表示当前同化时刻，w表示同化窗口内的积分步数，W表示同化窗口的长度，F
t
，F
t+1
，
…
，F
t+w
分别表示在积分时刻t，t+1，
…
，t+w时的驱动数据，М表示陆面过程模型，M(X
Nm
,F
t
)表示将X
Nm
输入М并在F
t
驱动下向前积分一步，X
t+wNm
＝[x
t+w1m
,
…
,x
t+wnm
,
…
,x
t+wNm
]
′
表示将X
Nm
输入М并在F
t
，F
t+1
，
…
，F
t+w
驱动下向前积分w步的预报结果，x
t+w1m
，
…
，x
t+wnm
，
…
，x
t+wNm
表示X
t+wNm
的第1
…
，n，
…
，N维校正变量的预报结果。步骤3012：由步骤3011中的X
t+wNm
构成校正向量预报集合X
TNM
：χ
t+wNM
＝{X
t+wN1
,
…
,X
t+wNm
,
…
,X
t+wNM
}X
TNM
＝{χ
tNM
,
…
,χ
t+wNM
…
,χ
t+WNM
}步骤3013：将X
TNM
中每个成员输入到观测算子并结合观测算子参数数据得到同化向量的预报集合，公式为：
Y
t+wLm
＝H(X
t+wNm
,P
t+w
)0≤w≤W 1≤m≤MΓ
t+wLM
＝{Y
t+wL1
,
…
,Y
t+wLM
}Y
TLM
＝{Γ
tLM
,
…
,Γ
t+WLM
}式中，H表示观测算子，L表示同化向量的维数，P
t+w
表示在t+w时观测算子的参数数据，H(X
t+wNm
,P
t+w
)表示将X
t+wNm
输入H并结合P
t+w
得到t+w时同化向量的预报结果Y
t+wLm
＝[y
t+w1m
,
…
,y
t+wlm
,
…
,y
t+wLm
]
′
，y
t+w1m
，
…
，y
t+wlm
，
…
，y
t+wLm
表示第1
…
，l，
…
，L维同化变量的预报结果，Γ
t+wLM
为Y
t+wL1
，
…
，Y
t+wLM
组成的预报集合，Y
TLM
为由Γ
tLM
，
…
，Γ
t+WLM
组成的同化向量预报集合的集合。4.如权利要求1所述的一种多变量频域同化方法，其特征在于，所述步骤4的具体步骤为：步骤4011：利用对应同化窗口内各积分时刻的同化向量观测产品的历史均值向量和历史方差向量生成同化向量观测产品集合，生成具体过程为：a)产生(W+1)
×
L
×
M个在(0，1)上均匀分布的随机变量v1
wlm
和v2
wlm
，w，l，m表示随机变量的下标变量；b)利用v1
wlm
和v2
wlm
产生(W+1)
×
L
×
M个服从标准正态分布的随机变量V
wlm
，公式如下：c)利用V
wlm
产生W+1个大小为M的服从L维正态分布的同化向量观测产品的集合O
TLM
：o
t+wlm
＝V
wlm
×
δ
t+wl
+ν
t+wl
式中，δ
t+wl
，ν
t+wl
分别为同化向量中第l维同化变量观测对应t+w时的历史方差和历史均值，o
t+wlm
为产生的第l维同化变量的第m个随机变量；步骤4012：利用校正向量预报集合、同化向量预报集合及观测产品集合计算同化窗口内各积分时刻的卡尔曼增益矩阵，计算具体过程为：a)利用X
TNM
和Y
TLM
计算同化窗口内各积分时刻校正向量与同化向量之间的协方差矩阵，计算公式为：计算公式为：计算公式为：计算公式为：计算公式为：
式中，P
t+wn
H
t+wi
为在t+w时第n维校正变量与第i维同化变量之间的预报协方差，为在t+w时N维校正向量与L维同化向量之间的预报协方差矩阵，P
TN
H
TL
为N维校正向量与L维同化向量之间的预报协方差矩阵的集合；b)利用Y
TLM
计算各积分时刻同化向量预报结果的协方差矩阵，计算公式为：计算各积分时刻同化向量预报结果的协方差矩阵，计算公式为：计算各积分时刻同化向量预报结果的协方差矩阵，计算公式为：计算各积分时刻同化向量预报结果的协方差矩阵，计算公式为：计算各积分时刻同化向量预报结果的协方差矩阵，计算公式为：式中，H
t+wi
H
t+wl
为在t+w时第i维与第l维同化变量之间的预报协方差，为在t+w时L个同化变量之间的预报协方差矩阵，H
TL
H
TL
为同化变量之间的预报协方差矩阵的预报集合；c)利用O
TLM
计算同化窗口内各积分时刻同化向量观测产品的协方差矩阵，计算公式为：计算同化窗口内各积分时刻同化向量观测产品的协方差矩阵，计算公式为：计算同化窗口内各积分时刻同化向量观测产品的协方差矩阵，计算公式为：计算同化窗口内各积分时刻同化向量观测产品的协方差矩阵，计算公式为：计算同化窗口内各积分时刻同化向量观测产品的协方差矩阵，计算公式为：式中，O
t+wi
O
t+wl
为在t+w时第i维与第l维同化变量之间的观测协方差，为在t+w时同化变量之间的观测协方差矩阵，O
TL
O
TL
为观测协方差矩阵的集合；d)利用P
TN
H
TL
、H
TL
H
TL
及O
TL
O
TL
计算同化窗口内各积分时刻的卡尔曼增益矩阵计算公式如下：O
t+wi
H
t+wl
＝O
t+wi
O
t+wl
+H
t+wi
H
t+wl
步骤4013：将同化窗口内各积分时刻同化向量预报集合的均值向量、同化向量观测产品串联成时序信号，串联过程如下：a)将Γ
tLM
，
…
，Γ
t+WLM
的预报均值向量串联成时序信号，串联方法为：Y
Wl
＝[Y
tl
,...,Y
t+wl
,...,Y
t+Wl
]式中，Y
tl
，
…
，Y
t+wl
，
…
，Y
t+Wl
分别为各积分时刻第l维同化变量的预报均值，Y
Wl
为同化窗口内第l维同化变量预报均值的时序信号；b)将各积分时刻同化向量观测产品O
tL
＝[o
t1
,
…
,o
tl
,
…
,o
tL
]
′
，
…
，O
t+wL
＝[o
t+w1
,
…
,o
t+wl
,
…
,o
t+wL
]
′
，
…
，O
t+WL
＝[o
t+W1
,
…
,o
t+Wl
,
…
,o
t+WL
]
′
串联成时序信号，串联方法为：o
Wl
＝[o
tl
,...,o
t+wl
,...,o
t+Wl
]1≤l≤LO
WL
＝[o
W1
,...,o
Wl
,...,o
WL
]'式中，o
Wl
为第l维同化变量的观测值所组成的时序信号，O
Wl
为同化向量观测产品的时序信号矩阵；步骤4014：将同化窗口内各积分时刻的卡尔曼增益矩阵串联成卡尔曼增益矩阵时序信号，串联方法为：g
Wnl
＝[G
tnl
,...,G
t+wnl
,...,G
t+Wnl
]1≤n≤N 1≤l≤L式中，g
Wnl
为由卡尔曼增益矩阵G
t+wnl
，t≤t+w≤t+W集合的第n行第l列元素组成的时序信号；步骤4015：将同化向量预报集合的均值向量的时序信号和同化变量观测产品的时序信号进行频率分解，得到各自的频率成份，频率分解过程如下：a)将Y
Wl
进行频率分解，得到频率成份，频率分解方法如下：式中，Ry
Wl
[j]为Y
Wl
中频率为j的余弦谐波系数值，Iy
Wl
[j]为Y
Wl
中频率为j的正弦谐波系数值；b)将o
Wl
进行频率分解，得到频率成份，频率分解方法如下：式中，Ro
Wl
[j]为o
Wl
中频率为j的余弦谐波系数值，Io
Wl
[j]为o
Wl
中频率为j的正弦谐波系数值；步骤4016：将g
Wnl
进行频率分解，得到卡尔曼增益矩阵各个矩阵元素时序信号的频率成份，频率分解方法如下：
式中，Rg
Wnl
[j]为g
Wnl
中频率为j的余弦谐波系数值，Ig
Wnl
[j]为g
Wnl
中频率为j的正弦谐波系数值。步骤4017：在不同频率成份上计算校正向量的校正量，并将所有频率成份上校正向量的校正量合成为同化窗口内各积分时刻的校正向量的校正量，合成过程如下：a)在不同频率成份上计算校正向量各维校正变量的校正量，计算方法为：a)在不同频率成份上计算校正向量各维校正变量的校正量，计算方法为：式中，Rκ
Wn
[j]为第n维校正变量时序信号中频率为j的余弦谐波的校正量，Iκ
Wn
[j]为第n维校正变量时序信号中频率为j的正弦谐波的校正量；b)将所有频率成份0≤j≤W/2上各维校正变量的校正量合成为同化窗口内各积分时刻各维校正变量的校正量，合成公式如下：式中，ξ
t+wn
为在t+w时第n维校正变量的校正量；步骤4018：计算同化窗口内各积分时刻的校正向量校正量的均值向量，计算方法为：Ξ
N
＝[ξ1,...,ξ
n
,...,ξ
N
]'式中，ξ
n
为所有积分时刻第n维校正变量的校正量的均值，Ξ
N
为所有N维校正变量的校正量的均值向量；步骤4019：利用校正向量校正量的均值向量对校正向量进行校正，具体为：X
tN
＝X
tN
+Ξ
N
式中，X
tN
为同化时刻t时校正向量的值及校正后的同化值。5.一种实现权利要求1
‑
4任一项一种多变量频域同化方法的系统，其特征在于，包括同化准...
【专利技术属性】
技术研发人员：陈少辉，
申请(专利权)人：中国科学院地理科学与资源研究所，
类型：发明
国别省市：
还没有人留言评论。发表了对其他浏览者有用的留言会获得科技券。