Page 103 - 2023年第54卷第10期

{ １ ，ａ ≤θ ｉ ≤ｂ
ｉ
ｉ
） ＝ ｂ － ａ （７）
Ｐ（ θ ｉ ｉ ｉ
０，其他
先验分布的下限和上限。
式中 ａ、ｂ分别为参数 θ ｉ
ｉ ｉ
如果各参数独立，总先验分布 Ｐ（ θ ）可定义为：
ｎ
Ｐ（ θ ） ＝ ∏ Ｐ（ θ ｉ ） （８）
ｉ ＝１
， ）服从 ｎ
根据水头观测值 ｄ与模拟值 Ｆ（ θ ）之间的总误差（ ε ＝ ｄ － Ｆ（ θ ）），假设 ε ＝ （ ε １ ε ２ ，…，ε ｎ
维高斯分布，得出 ε 的条件概率密度函数（似然函数）ｐ（ｄ "θ ）为：
１ １
ｐ（ｄ "θ ） ＝ １?２ [ Ｔ － １ ] （９）
ｅｘｐ － （ｄ － Ｆ（ θ ））Ｃ （ｄ － Ｆ（ θ ））
ｄ
（２ π ） ｎ?２ Ｃ ｄ ２
－ １
式中：Ｃ 为总误差协方差矩阵； Ｃ 为矩阵 Ｃ 的行列式；Ｃ 为矩阵 Ｃ 的逆。
ｄ
ｄ
ｄ
ｄ
ｄ
由于 ｐ（ｄ "θ ）计算中需要调用正演模型 Ｆ（ θ ），当参数维度较高时，一般无法得到 ｐ（ θ"ｄ）的解析
解。可采用基于 ＭＣＭＣ的概率抽样方法和 ＥＳ － ＭＤＡ数据同化方法，推求 ｐ（ θ"ｄ）的统计量（如均值、
标准差等），当作参数反演问题的解。
２．４．２ 基于 ＭＣＭＣ的参数后验分布求解 ＭＣＭＣ通过随机采样构建一个马尔科夫链，当马尔科夫链收敛
后的采样就相当于在后验分布 ｐ（ θ" ｄ）上抽取样本，利用样本的统计量就可以得到参数均值和标准差。
常用实现 ＭＣＭＣ方式是 Ｍｅｔｒｏｐｏｌｉｓ － Ｈａｓｔｉｎｇｓ（Ｍ－ Ｈ）算法 ［２１］ ，其基本思想是通过构造一个接受 － 拒
绝准则，以保证采到的样本服从参数后验分布 ｐ（ θ"ｄ），具体步骤如下：
，ｉ ＝ ０ ；
（ １）初始化参数。采用蒙特卡罗方法，随机生成一组参数初值 θ ｉ
（２）生成参数的候选值。从建议分布（如正态分布）ｑ（ θ ）中随机生成一个候选值 θ 。

（３）计算接受概率：
{ ｌｏｇ（ｐ（ｄ "θ ））ｑ（ θ ｉ ，θ ） }


，θ ） ＝ｍｉｎ１， （１０）
α （ θ ｉ 

ｌｏｇ（ｐ（ｄ "θ ｉ ））ｑ（ θ ，θ ｉ ）
）。
式中：ｌｏｇ（ｐ（ｄ "θ ））为对数似然函数；在建议分布为正态分布时，ｑ（ θ ｉ ，θ ） ＝ｑ（ θ ，θ ｉ


， ＝
（４）接受或拒绝候选状态。从均匀分布 Ｕ（０，１）中抽取随机数 ｕ，若 ｕ ≤α （ θ ｉ θ ），则令 θ ｉ ＋ １

＝ 。

θ ，否则令 θ ｉ ＋ １ θ ｉ
（５）令 ｉ ＝ ｉ ＋ １ ，返回步骤（２），重复多次，直到采样得到足够的参数样本。
，
２．４．３ 基于 ＥＳ － ＭＤＡ的参数后验分布求解 ＥＳ － ＭＤＡ在每次迭代中对观测误差引进一个膨胀系数 α ｉ
以较低的计算代价解决高维非线性的参数反演问题，便于实现多次数据同化。具体步骤如下：
满足：
ａ
（ １）给定数据同化次数 Ｎ，集合样本 Ｎ，第 ｉ次数据同化时的膨胀系数 α ｉ ，其中 α ｉ
Ｎ ａ
１
∑ ＝１，ｉ ＝１，２，…，Ｎ ａ （１１）
ｉ ＝１α ｉ
：
（２）根据各参数的先验分布 ｐ（ θ ）抽取 Ｎ个样本，组成初始参数集合 θ １
］ （１２）
θ １ ＝ ［ θ １，１ ，θ １，２ ，…，θ １，Ｎ
（３）初始 ｉ ＝ １ ，对集合中每个样本进行正演模型计算，得到模拟值集合：
），ｊ ＝ １，２，…，Ｎ （１３）
ｄ ＝ Ｆ（ θ ｉ，ｊ
ｉ，ｊ
（ ４）对样本参数集合进行更新：
１?２
ｄ ＝ ｄ ＋ α ｉ Ｒ ｚ （１４）
ｕｃ，ｊ ｏｂｓ，ｊ 槡
－ １
＝ ＋ Ｃ （Ｃ ＋ Ｒ） （ｄ － ｄ ） （１５）
θ ｉ ＋ １ ，ｊ θ ｉ，ｊ ＭＤ ＤＤ ｕｃ， ｊ ｉ，ｊ
扰动后的观测值；ｚ为修正引入服从高斯分布的随机误差；Ｃ 为模
ｕｃ
式中：ｄ 为观测值；ｄ 为添加 α ｉ ＭＤ
ｏｂｓ
与模拟值 ｄ之间的协方差矩阵；Ｃ 与 Ｒ分别为模拟值和观测误差的自协方差矩阵。
型参数 θ ｉ ｉ ＤＤ
，
＝ ［ θ Ｎ ａ ，１
ａ
（ ５）ｉ ＝ ｉ ＋ １，重复执行步骤（３）和（４）直到完成 Ｎ 次数据 同 化，得到 最终 样 本 集合 θ Ｎ ａ
２
— １ ３ ９ —