CN117035095A - 一种含拉普拉斯噪声的偏微分方程反问题的推断方法 - Google Patents

Abstract

本发明公开了一种含拉普拉斯噪声的偏微分方程反问题的推断方法，包括：假设观测噪声服从拉普拉斯分布，将噪声的概率密度函数表示为高斯尺度的混合；获取似然函数；利用马尔科夫随机场表示未知变量的先验分布、多层建模把尺度参数建模为随机变量；通过拉普拉斯分布的高斯尺度混合表示将反问题的后验概率密度函数重新建模为多层贝叶斯模型；基于Kullback‑Leibler散度，采用变分贝叶斯推断方法将后验概率密度函数的计算问题转化为概率分布的变分问题；施加平均场假设，在概率分布族中寻找真实后验概率分布的近似估计，反演出未知变量、噪声方差参数和尺度参数的近似分布；选取适当的逼近概率分布族，得到非退化方法和退化方法。

Description

一种含拉普拉斯噪声的偏微分方程反问题的推断方法

技术领域

本发明涉及反问题推断技术领域，特别是涉及一种含拉普拉斯噪声的偏微分方程反问题的推断方法。

背景技术

随着科学技术水平的进步和数学研究的不断深入，反问题已经被广泛应用于地质与环境科学、生命科学与医学、地球物理科学、工程控制、信号处理等诸多科学技术领域。在数学物理研究领域，反问题用于描述时空域逆时针的物理变化，可通过直接观测的信息来探求事物内部规律。例如在地球物理勘探领域，接收地震波信号，通过走时成像得到地震波在不同深度的传播速度，进一步利用声波方程或单程波方程偏移成像方法，得到反射界面的位置和形状等地层地貌形态，这类问题就是数学物理反问题。

反问题的求解相对困难，本质在于反问题通常是不适定的，即解的存在性、唯一性和稳定性中的一个或多个不满足。针对不适定问题，研究学者提出了两大类有效解决方法：正则化方法和贝叶斯方法。求解反问题的传统正则化方法的出发点是测量数据的误差特征己知或误差很小可忽略不计。然而在处理实际问题时，难以精确获得测量数据的误差特征，或难以获得测量数据的具体误差大小。此外，从统计方面讲，正则化方法求解得到的是一个确定的值，相当于贝叶斯统计方法的“点估计”，它给出的信息较少，忽略了误差的随机性和模型的不确定性信息，因此正则化方法存在一定的局限性。为了克服上述缺点，研究学者使用随机方法来量化模型的不确定性并求解不适定问题，其中最常用的是贝叶斯统计推断方法，简称为贝叶斯方法。贝叶斯方法将反问题转化为统计推断问题，利用先验概率刻画待估计函数的已知信息，借助贝叶斯公式有效地将已知信息和由测量数据获得的信息相结合，从而得到能够对待估计函数不确定性进行完整刻画的后验概率分布。相比于传统正则化方法，贝叶斯方法提供了不确定性量化分析。

贝叶斯反演方法中一个重要问题是后验信息的获取，以马尔科夫链蒙特卡罗(MarkovchainMonteCarlo，MCMC)采样为代表的采样算法在统计反演中获得了广泛关注，但在涉及海量数据或者模型较复杂的情况下，其计算成本过高。此时，被广泛应用于机器学习的变分推断法就是一种很好的替代方法，它在与点估计计算量相当的情况下能够给出待估计函数的不确定性分析，因此对于偏微分方程反问题，该方法具有广泛的适用性。

变分贝叶斯推断为解决各种反问题提供了一个灵活性框架，首先，它能够产生一组与观测数据一致的解的集合，从而实现解的不确定性量化，其次，通过分层建模解决了超参数选择的难题，从而提供了一种灵活的正则化方法。在现有的反问题贝叶斯推断研究中，噪声先验分布大多假设为高斯型分布。虽然高斯型先验分布使得计算方便并易于探索后验状态空间，但是高斯型分布的局限性在于它对观测数据中的异常值缺乏鲁棒性，一个异常数据点就可以显著影响模型中的所有参数的估计值，因此严重损害估计的准确性。另外，并不是所有的数据噪声都能被高斯模型充分描述，有时在获取信号时可能会产生拉普拉斯噪声。

发明内容

为解决上述背景技术中存在的至少一种问题，本发明提供一种含拉普拉斯噪声的偏微分方程反问题的推断方法。

本发明实施例提供的一种拉普拉斯噪声的偏微分方程反问题的推断方法，包括如下八个步骤：

步骤一：假设观测噪声服从拉普拉斯分布，将噪声的概率密度函数表示为高斯尺度的混合；

步骤二：获取似然函数；

步骤三：利用马尔科夫随机场表示未知变量的先验分布；

步骤四：利用多层建模把尺度参数建模为随机变量；

步骤五：通过拉普拉斯分布的高斯尺度混合表示将反问题的后验概率密度函数重新建模为多层贝叶斯模型；

步骤六：基于Kullback-Leibler散度，采用变分贝叶斯推断方法将后验概率密度函数的计算问题转化为关于概率分布的变分问题；

步骤七：施加平均场假设，在概率分布族中寻找真实后验概率分布的近似估计，反演出未知变量、噪声方差参数和尺度参数的近似分布；

步骤八：选取适当的逼近概率分布族，得到能够同时刻画未知变量、噪声方差参数和尺度参数的不确定性的非退化方法，以及仅刻画噪声方差参数和尺度参数的不确定性的退化方法。

进一步地，上述一种含拉普拉斯噪声的偏微分方程反问题的推断方法中，步骤一中假设观测噪声服从拉普拉斯分布，将噪声的概率密度函数表示为高斯尺度的混合，包括如下两个步骤：

步骤(1)建立含拉普拉斯噪声反问题的变分贝叶斯模型，对于有限维线性反问题：Hm＝d，其中分别代表系统矩阵，未知变量，观测数据；设d＝d⁺+ζ，其中d⁺表示无噪声数据，ζ表示噪声,ζ＝[ζ₁,ζ₂,···ζ_n]^T，ζ的每个分量ζ_i是独立同分布且服从期望为0，方差为的拉普拉斯分布，记为其概率密度函数表示为：

步骤(2)将噪声ζ的概率密度函数表示为高斯尺度的混合，具体形式为：

其中，ζ_i服从期望为0、方差为w_i的高斯分布，记为w_i是隐变量，w_i服从指数分布w_i～Exponential(τ)，其概率密度函数为w＝[w₁,w₂,···w_n]^T，w的概率密度函数

进一步地，上述一种含拉普拉斯噪声的偏微分方程反问题的推断方法中，步骤二中获取似然函数包括：

在考虑加性噪声模型时，似然函数是噪声ζ概率密度函数的平移，可表示为：即已知m和w时，d的概率密度函数，其中，m是未知变量，d是观测数据，w是服从指数分布的随机变量；d服从期望为Hm、协方差矩阵为diag(w₁,w₂,···w_n)的n维高斯分布，w_i是隐变量，w_i服从指数分布w_i～Exponential(τ)。

进一步地，上述一种含拉普拉斯噪声的偏微分方程反问题的推断方法中，步骤三中利用马尔科夫随机场表示未知变量的先验分布，包括：

采用先验建模工具马尔科夫随机场表示未知变量m的先验分布p(m|λ)，其形式为：其中，λ是描述交互作用强度的尺度参数，矩阵L表征相邻点之间交互作用的结构，s为矩阵L的秩；

步骤四中利用多层建模把尺度参数建模为随机变量，包括：选取伽马分布G(λ；α,β)为尺度参数λ的先验分布：其中，Γ(·)为标准伽马函数，α和β是非负常数。

进一步地，上述一种含拉普拉斯噪声的偏微分方程反问题的推断方法中，步骤五中通过拉普拉斯分布的高斯尺度混合表示将反问题的后验概率密度函数重新建模为多层贝叶斯模型，包括：利用多层贝叶斯建模将后验概率密度函数表示为：

其中，w是噪声方差参数，W是对角元素为w的对角矩阵，范数定义为(α,β)为伽马分布的参数对，m是未知变量，d是观测数据，λ是描述交互作用强度的尺度参数，矩阵L表征相邻点之间交互作用的结构，s为矩阵L的秩，α和β是非负常数，d服从期望为Hm、协方差矩阵为diag(w₁,w₂,···w_n)的n维高斯分布，w_i是隐变量，w_i服从指数分布w_i～Exponential(τ)；

步骤六中基于Kullback-Leibler散度，采用变分贝叶斯推断方法将后验概率密度函数的计算问题转化为关于概率分布的变分问题，包括：

将后验概率密度函数转化为求解优化问题在易计算的概率分布族中寻找近似分布q(m,w,λ)以逼近真实的后验概率密度函数p(m,w,λ|d)。

进一步地，上述一种含拉普拉斯噪声的偏微分方程反问题的推断方法中，步骤八中选取适当的逼近概率分布族，得到能够同时刻画未知变量、噪声方差参数和尺度参数的不确定性的非退化方法，包括：对未知变量m、噪声方差参数w和尺度参数λ施加平均场假设，近似分布q(m,w,λ)可分离为q(m)q(w)q(λ)，平均场假设同时刻画未知变量m、噪声方差参数w和尺度参数λ的不确定性，此条件下的近似称为非退化近似；

步骤八中选取适当的逼近概率分布族，得到能够仅刻画噪声方差参数w和尺度参数λ的不确定性的退化方法，包括：增加额外假设以降低计算复杂性，q(m)退化的情况下，即q(m)＝1时，仅能刻画噪声方差参数w和尺度参数λ不确定性，此条件下的近似称为退化近似。

进一步地，上述一种含拉普拉斯噪声的偏微分方程反问题的推断方法中，非退化近似，包括以下步骤：

(1)计算近似分布的解为q(m,w,λ)＝q^*(m)q^*(w)q^*(λ)，其中

(2)计算后验概率密度函数的对数函数ln p(m,w,λ|d)：

(3)反演未知变量m：

将(2)中公式代入(1)中第一个式子有：

其中，X(w,λ)是与m无关的项；

令则(3)前两项简化为：

再令m^*＝(H^TH+η^*L^TL)^-1H^Td，则上式为：

其中，最后一项对归一化常数有贡献，m服从均值为m^*，方差为(H^TW^*H+λ^*L^TL)^-1的高斯分布，

即：

(4)反演噪声方差参数w：

将(2)中公式代入(1)中第二个式子有：

其中，T(m,λ)是ln p(m,w,λ)中与w无关的项；

作变换z服从均值为μ_z，方差为λ_z的反高斯分布，即：其中，

(5)反演尺度参数λ：将(2)中公式代入(1)中第三个式子有：

其中，F(m,w)是与λ无关的项，λ服从均值为协方差为的伽马分布，即：

(6)对参数τ进行假设，参数τ与隐变量w_i直接相关，具体应用时仔细调整，以获得合理的结果；

(7)通过步骤(3)～步骤(5)参数反演，得到近似分布q(m,w,λ)的各变量边缘分布的非显式解析解；在分别求解q(m)、q(w)和q(λ)时，采用交替迭代法解决，即首先初始化q(m)、q(w)、q(λ)和q(m)、q(w)、q(λ)中涉及到的参数，然后从q(m)、q(w)、q(λ)中其中一个以及涉及到的参数开始交替迭代计算，直到算法收敛；交替迭代型算法对q(m)、q(w)和q(λ)估计；

其中，w是服从指数分布的随机变量；d服从期望为Hm、协方差矩阵为diag(w₁,w₂,···w_n)的n维高斯分布，w_i是隐变量，w_i服从指数分布w_i～Exponential(τ)，W是对角元素为w的对角矩阵，范数定义为(α,β)为伽马分布的参数对，m是未知变量，d是观测数据，λ是描述交互作用强度的尺度参数，矩阵L表征相邻点之间交互作用的结构，s为矩阵L的秩，α和β是非负常数。

进一步地，上述一种含拉普拉斯噪声的偏微分方程反问题的推断方法中，退化近似，包括以下步骤：

(1)进一步降低计算复杂性,在对m、w和λ施加条件独立性假设的基础上增加q(m)退化的假设，即q(m)＝1，对于固定的m，考虑最小化KL散度可以转化为公式：引入KL散度的拉格朗日函数定义如下所示：

通过计算，求解得到：

m^*＝(H^TH+λ^*L^TL)^-1H^Td

类似非退化近似中的计算，得到：

(2)具体估计噪声方差参数w，

令得到z服从均值为μ_z，方差为λ_z的反高斯分布，即

(3)具体估计尺度参数λ，

λ服从均值为协方差为的伽马分布，

即其中m^*＝(H^TH+λ^*L^TL)^-1H^Td；

(4)采用交替迭代型算法对q(w)和q(λ)进行退化近似估计；

其中，w是服从指数分布的随机变量；d服从期望为Hm、协方差矩阵为diag(w₁,w₂,···w_n)的n维高斯分布，w_i是隐变量，w_i服从指数分布w_i～Exponential(τ)，W是对角元素为w的对角矩阵，范数定义为(α,β)为伽马分布的参数对，m是未知变量，d是观测数据，λ是描述交互作用强度的尺度参数，矩阵L表征相邻点之间交互作用的结构，s为矩阵L的秩，α和β是非负常数。

本发明有益效果在于：本发明通过假设噪声服从拉普拉斯分布，使用拉普拉斯分布的高斯尺度混合表示将反问题重新表示为多层贝叶斯模型，同时反演出解以及噪声方差参数和尺度参数的近似分布；通过选取不同的逼近概率分布族，给出能够同时够刻画解与超参数不确定性的非退化方法以及仅刻画超参数不确定性的退化方法。迭代循环求解，实现在更高估计精度和准确性的前提下有更快的估计速度和更强的鲁棒性，充分克服高斯噪声假设下，对异常数据敏感从而鲁棒性较差的问题，更适合于有实时性要求的工程领域应用。

附图说明

为了更清楚地说明本发明实施例或传统技术中的技术方案，下面将对实施例或传统技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提供的一种含拉普拉斯噪声的偏微分方程反问题的推断方法流程图；

图2为计算区域Ω；

图3是当θ(x)＝sin(πx)e^π+x+1，噪声方差参数为0.5的高斯噪声时，高斯模型和拉普拉斯模型下非退化和退化模型算法的实现效果；

图4是当θ(x)＝sin(πx)e^π+x+1，噪声方差参数为5的拉普拉斯噪声时，高斯模型和拉普拉斯模型下非退化和退化模型算法的实现效果；

图5是当噪声方差参数为0.5的高斯噪声时，高斯模型和拉普拉斯模型下非退化和退化模型算法的实现效果；

图6是当噪声方差参数为5的拉普拉斯噪声时，高斯模型和拉普拉斯模型下非退化和退化模型算法的实现效果。

具体实施方式

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图对本发明的具体实施方式做详细的说明。在下面的描述中阐述了很多具体细节以便于充分理解本发明。但是本发明能够以很多不同于在此描述的其它方式来实施，本领域技术人员可以在不违背本发明内涵时做类似改进，因此本发明不受下面公开的具体实施的限制。

除非另有定义，本文所使用的所有的技术和科学术语与属于本发明的技术领域的技术人员通常理解的含义相同。本文中在本发明的说明书中所使用的术语只是为了描述具体的实施例的目的，不是旨在于限制本发明。本文所使用的术语“及/或”包括一个或多个相关的所列项目的任意的和所有的组合。

图1为本发明实施例提供的一种含拉普拉斯噪声的偏微分方程反问题的推断方法流程图。

本发明实施例提供的一种拉普拉斯噪声的偏微分方程反问题的推断方法，结合图1，包括S101至S108八个步骤：

S101：假设观测噪声服从拉普拉斯分布，将噪声的概率密度函数表示为高斯尺度的混合。

S102：获取似然函数；

S103：利用马尔科夫随机场表示未知变量的先验分布；

S104：利用多层建模把尺度参数建模为随机变量；

S105：通过拉普拉斯分布的高斯尺度混合表示将反问题的后验概率密度函数重新建模为多层贝叶斯模型；

S106：基于Kullback-Leibler散度，采用变分贝叶斯推断方法将后验概率密度函数的计算问题转化为关于概率分布的变分问题；

S107：施加平均场假设，在概率分布族中寻找真实后验概率分布的近似估计，反演出未知变量、噪声方差参数和尺度参数的近似分布；

S108：选取适当的逼近概率分布族，得到能够同时刻画未知变量、噪声方差参数和尺度参数的不确定性的非退化方法，以及仅刻画噪声方差参数和尺度参数的不确定性的退化方法。

进一步地，上述一种含拉普拉斯噪声的偏微分方程反问题的推断方法中，步骤S101中假设观测噪声服从拉普拉斯分布，将噪声的概率密度函数表示为高斯尺度的混合，包括如下两个步骤：

步骤(1)建立含拉普拉斯噪声反问题的变分贝叶斯模型，对于有限维线性反问题：Hm＝d，其中分别代表系统矩阵，未知变量，观测数据；设d＝d⁺+ζ，其中d⁺表示无噪声数据，ζ表示噪声,ζ＝[ζ₁,ζ₂,···ζ_n]^T，ζ的每个分量ζ_i是独立同分布且服从期望为0，方差为的拉普拉斯分布，记为其概率密度函数表示为：

步骤(2)将噪声ζ的概率密度函数表示为高斯尺度的混合，具体形式为：

其中，ζ_i服从期望为0、方差为w_i的高斯分布，记为w_i是隐变量，w_i服从指数分布w_i～Exponential(τ)，其概率密度函数为w＝[w₁,w₂,···w_n]^T，w的概率密度函数

进一步地，上述一种含拉普拉斯噪声的偏微分方程反问题的推断方法中，步骤S102中获取似然函数包括：

在考虑加性噪声模型时，似然函数是噪声ζ概率密度函数的平移，可表示为：即已知m和w时，d的概率密度函数，其中，m是未知变量，d是观测数据，w是服从指数分布的随机变量；d服从期望为Hm、协方差矩阵为diag(w₁,w₂,···w_n)的n维高斯分布，w_i是隐变量，w_i服从指数分布w_i～Exponential(τ)。

进一步地，上述一种含拉普拉斯噪声的偏微分方程反问题的推断方法中，步骤S103中利用马尔科夫随机场表示未知变量的先验分布，包括：

采用先验建模工具马尔科夫随机场表示未知变量m的先验分布p(m|λ)，其形式为：其中，λ是描述交互作用强度的尺度参数，矩阵L表征相邻点之间交互作用的结构，s为矩阵L的秩；

步骤S104中利用多层建模把尺度参数建模为随机变量，包括：选取伽马分布G(λ；α,β)为尺度参数λ的先验分布：其中，Γ(·)为标准伽马函数，α和β是非负常数。

进一步地，上述一种含拉普拉斯噪声的偏微分方程反问题的推断方法中，步骤S105中通过拉普拉斯分布的高斯尺度混合表示将反问题的后验概率密度函数重新建模为多层贝叶斯模型，包括：利用多层贝叶斯建模将后验概率密度函数表示为：

其中，w是噪声方差参数，W是对角元素为w的对角矩阵，范数定义为为伽马分布的参数对，m是未知变量，d是观测数据，λ是描述交互作用强度的尺度参数，矩阵L表征相邻点之间交互作用的结构，s为矩阵L的秩，α和β是非负常数，d服从期望为Hm、协方差矩阵为diag(w₁,w₂,···w_n)的n维高斯分布，w_i是隐变量，w_i服从指数分布w_i～Exponential(τ)；

步骤S106中基于Kullback-Leibler散度，采用变分贝叶斯推断方法将后验概率密度函数的计算问题转化为关于概率分布的变分问题，包括：

将后验概率密度函数转化为求解优化问题在易计算的概率分布族中寻找近似分布q(m,w,λ)以逼近真实的后验概率密度函数p(m,w,λ|d)。

进一步地，上述一种含拉普拉斯噪声的偏微分方程反问题的推断方法中，步骤S108中选取适当的逼近概率分布族，得到能够同时刻画未知变量、噪声方差参数和尺度参数的不确定性的非退化方法，包括：对未知变量m、噪声方差参数w和尺度参数λ施加平均场假设，近似分布q(m,w,λ)可分离为q(m)q(w)q(λ)，平均场假设同时刻画未知变量m、噪声方差参数w和尺度参数λ的不确定性，此条件下的近似称为非退化近似；

步骤S108中选取适当的逼近概率分布族，得到能够仅刻画噪声方差参数w和尺度参数λ的不确定性的退化方法，包括：增加额外假设以降低计算复杂性，q(m)退化的情况下，即q(m)＝1时，仅能刻画噪声方差参数w和尺度参数λ不确定性，此条件下的近似称为退化近似。

进一步地，上述一种含拉普拉斯噪声的偏微分方程反问题的推断方法中，非退化近似，包括以下步骤：

(1)计算近似分布的解为q(m,w,λ)＝q^*(m)q^*(w)q^*(λ)，其中

(2)计算后验概率密度函数的对数函数lnp(m,w,λd)：

(3)反演未知变量m：

将(2)中公式代入(1)中第一个式子有：

其中，X(w,λ)是与m无关的项；

令则(3)前两项简化为：

再令m^*＝(H^TH+η^*L^TL)^-1H^Td，则上式为：

其中，最后一项对归一化常数有贡献，m服从均值为m^*，方差为(H^TW^*H+λ^*L^TL)^-1的高斯分布，

即：

(4)反演噪声方差参数w：

将(2)中公式代入(1)中第二个式子有：

其中，T(m,λ)是lnp(m,w,λ)中与w无关的项；

作变换z服从均值为μ_z，方差为λ_z的反高斯分布，即：其中，

(5)反演尺度参数λ：将(2)中公式代入(1)中第三个式子有：

其中，F(m,w)是与λ无关的项，λ服从均值为协方差为的伽马分布，即：

(6)对参数τ进行假设，参数τ与隐变量w_i直接相关，具体应用时仔细调整，以获得合理的结果；

(7)通过步骤(3)～步骤(5)参数反演，得到近似分布q(m,w,λ)的各变量边缘分布的非显式解析解；在分别求解q(m)、q(w)和q(λ)时，采用交替迭代法解决，即首先初始化q(m)、q(w)、q(λ)和q(m)、q(w)、q(λ)中涉及到的参数，然后从q(m)、q(w)、q(λ)中其中一个以及涉及到的参数开始交替迭代计算，直到算法收敛；交替迭代型算法对q(m)、q(w)和q(λ)估计；

具体的，本发明实施例中，交替迭代型算法对q(m)、q(w)和q(λ)估计的步骤具体如下：

(7a)初始化q(m)、q(w)、q(λ)和q(m)、q(w)、q(λ)中涉及到的参数，输入初始值q⁰(w)和q⁰(λ)，令k＝0；

(7b)令k＝k+1；

(7c)计算q^k(m)

其中m^k＝(H^TW^kH+λ^kL^TL)^-1H^TW^kd，

(7d)计算q^k(w)

其中

(7e)计算q^k(λ)

(7f)设置停止准则；

(7g)输出q^k(m)q^k(w)q^k(λ)；

其中，w是服从指数分布的随机变量；d服从期望为Hm、协方差矩阵为diag(w₁,w₂,···w_n)的n维高斯分布，w_i是隐变量，w_i服从指数分布w_i～Exponential(τ)，W是对角元素为w的对角矩阵，范数定义为(α,β)为伽马分布的参数对，m是未知变量，d是观测数据，λ是描述交互作用强度的尺度参数，矩阵L表征相邻点之间交互作用的结构，s为矩阵L的秩，α和β是非负常数。

进一步地，上述一种含拉普拉斯噪声的偏微分方程反问题的推断方法中，退化近似，包括以下步骤：

(1)进一步降低计算复杂性,在对m、w和λ施加条件独立性假设的基础上增加q(m)退化的假设，即q(m)＝1，对于固定的m，考虑最小化KL散度可以转化为公式：引入KL散度的拉格朗日函数定义如下所示：

通过计算，求解得到：

m^*＝(H^TH+λ^*L^TL)^-1H^Td

类似非退化近似中的计算，得到：

(2)具体估计噪声方差参数w，

令得到z服从均值为μ_z，方差为λ_z的反高斯分布，即

m^*＝(H^TH+λ^*L^TL)^-1H^Td；

(3)具体估计尺度参数λ，

λ服从均值为协方差为的伽马分布，

即其中m^*＝(H^TH+λ^*L^TL)^-1H^Td；

(4)采用交替迭代型算法对q(w)和q(λ)进行退化近似估计；

具体的，本发明实施例中，采用交替迭代型算法对q(w)和q(λ)进行退化近似估计的步骤如下：

(4a)初始化q(m)、q(w)、q(λ)和q(m)、q(w)、q(λ)中涉及到的参数，输入初始值q⁰(w)和q⁰(λ)，令k＝0；

(4b)令k＝k+1；

(4c)计算m^k；

m^k＝(H^TH+λ^*L^TL)^-1H^Td

(4d)计算q^k(w)

其中

(4e)计算q^k(λ)

(4f)设置停止准则；

(4g)输出q^k(w)q^k(λ)；

其中，w是服从指数分布的随机变量；d服从期望为Hm、协方差矩阵为diag(w₁,w₂,···w_n)的n维高斯分布，w_i是隐变量，w_i服从指数分布w_i～Exponential(τ)，W是对角元素为w的对角矩阵，范数定义为(α,β)为伽马分布的参数对，m是未知变量，d是观测数据，λ是描述交互作用强度的尺度参数，矩阵L表征相邻点之间交互作用的结构，s为矩阵L的秩，α和β是非负常数。

以下结合一个Cauchy型问题的具体例子来说明含拉普拉斯噪声反问题的变分贝叶斯推断方法的性能。在热传导过程中，通常会出现用剩余部分的超定边界数据确定边界部分上未知边界条件的Cauchy型问题。例如，在重入式航天飞机的研究中，需要根据对内表面温度的测量来推断航天飞机外表面的温度。在淬火过程中，需要通过测量边界可接触部分的温度和流量来估计通过边界的对流热传递系数。

稳态热传导问题可用如下公式表示：

其中α(x)是电导率，f(x)是源项。

上式的限制如下：

其中，n是向外垂直于Γ_c的单位向量。反问题为寻找边界Γ_i上未知温度u(x)。

考虑Laplace方程，即令上式中α＝1,f＝0。图2为计算区域Ω，Ω定义为方形的单位开域(0,1)×(0,1)，实验测量不可接触的边界部分Γ_i和可接触的边界部分Γ_c，分别定义为：Γ_i＝[0,1]×{1}，Γ_c＝Γ\Γ_i，边界Γ₀定义为Γ₀＝{0,1}×(0,1)。利用3200个三角形有限元对该问题进行离散，在Γ₀上取80个测量值。

考虑重建光滑解的情况，令θ(x)＝sin(πx)e^π+x+1。当噪声为高斯噪声时，通过对比在不同噪声方差参数下假设为高斯噪声与拉普拉斯噪声的非退化和退化近似对真实解的重建，发现差别并不明显，都能实现较好的重建效果。细微差别在于两种噪声假设下的退化近似均比非退化近似算法更精确、更稳定。当噪声方差参数为0.5时，高斯模型和拉普拉斯模型下的实现效果如图3所示，消耗时间和重建误差对比如表1所示。图3中，黑色实线1为真实解，黑色实线2为假设噪声为高斯噪声的非退化近似，黑色虚线1为假设噪声为高斯噪声的退化近似，黑色虚线2为假设噪声为拉普拉斯噪声的非退化近似，黑色虚线3为假设噪声为拉普拉斯噪声的退化近似。在下图中，除非特殊说明，图中各线段含义相同。从图5和表1可得出四种近似均能以相对较小的误差快速重建真解，相比而言退化的近似更加精确，更加稳定。

表1

当噪声为拉普拉斯噪声时，通过对比在不同噪声方差参数下假设为高斯噪声与拉普拉斯噪声的非退化和退化近似对真实解的重建，发现假设噪声为高斯噪声的模型已不能很好的重建真解，而本发明中提出的拉普拉斯模型能以很高的精度重建出真解。当噪声方差参数为5时，高斯模型和拉普拉斯模型下的实现效果如图4所示，消耗时间和重建误差对比如表2所示。从图4和表2可得出噪声为拉普拉斯噪声时，高斯噪声模型的退化近似算法已完全不能重建真解，非退化近似算法的重建误差较大，而本发明中提出的两种拉普拉斯模型近似算法均以很高的精度和速度实现了真解的重建。

表2

考虑重建非光滑解的情况，令由于方程无法求得解析解，为了避免反演过失，将精细网格上正问题的数值解作为精确数据，在粗网格上进行反演。当噪声为高斯噪声时，通过对比在不同噪声方差参数下假设为高斯噪声与拉普拉斯噪声的非退化和退化近似对真实解的重建，发现差别并不明显，都能实现较好的重建效果。当噪声方差参数为0.5时，高斯模型和拉普拉斯模型下的实现效果如图5所示，消耗时间和重建误差对比如表3所示。从图5和表3可得出四种近似均能以相对较小的误差快速重建真解。

表3

当噪声为拉普拉斯噪声时，通过对比在不同噪声方差参数下假设为高斯噪声与拉普拉斯噪声的非退化和退化近似对真实解的重建，发现假设噪声为高斯噪声的模型已不能重建真解，而本发明中提出的拉普拉斯模型能以很高的精度重建出真解。当噪声方差参数为5时，高斯模型和拉普拉斯模型下的实现效果如图6所示，消耗时间和重建误差对比如表4所示。从图6和表4可得出噪声为拉普拉斯噪声时，高斯噪声模型的两种近似算法已完全不能重建真解，而本发明中提出的两种拉普拉斯模型近似算法均以很高的精度和速度实现了真解的重建。

表4

经实践，应用本发明，可充分克服高斯噪声假设下，模型对异常数据敏感和鲁棒性较差的问题，有效提高计算结果的准确性和鲁棒性。

本领域的技术人员能够理解，尽管在此所述的一些实施例包括其它实施例中所包括的某些特征而不是其它特征，但是不同实施例的特征的组合意味着处于本发明的范围之内并且形成不同的实施例。

本领域的技术人员能够理解，对各个实施例的描述都各有侧重，某个实施例中没有详述的部分，可以参见其他实施例的相关描述。

以上，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到各种等效的修改或替换，这些修改或替换都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。

Claims (8)

1.一种含拉普拉斯噪声的偏微分方程反问题的推断方法，其特征在于，包括如下八个步骤：

步骤一：假设观测噪声服从拉普拉斯分布，将噪声的概率密度函数表示为高斯尺度的混合；

步骤二：获取似然函数；

步骤三：利用马尔科夫随机场表示未知变量的先验分布；

步骤四：利用多层建模把尺度参数建模为随机变量；

步骤五：通过拉普拉斯分布的高斯尺度混合表示将反问题的后验概率密度函数重新建模为多层贝叶斯模型；

步骤六：基于Kullback-Leibler散度，采用变分贝叶斯推断方法将后验概率密度函数的计算问题转化为关于概率分布的变分问题；

步骤七：施加平均场假设，在概率分布族中寻找真实后验概率分布的近似估计，反演出未知变量、噪声方差参数和尺度参数的近似分布；

步骤八：选取适当的逼近概率分布族，得到能够同时刻画未知变量、噪声方差参数和尺度参数的不确定性的非退化方法，以及仅刻画噪声方差参数和尺度参数的不确定性的退化方法。

2.根据权利要求1所述的一种含拉普拉斯噪声的偏微分方程反问题的推断方法，其特征在于，所述步骤一中假设观测噪声服从拉普拉斯分布，将噪声的概率密度函数表示为高斯尺度的混合，包括如下两个步骤：

步骤(1)建立含拉普拉斯噪声反问题的变分贝叶斯模型，对于有限维线性反问题：Hm＝d，其中 分别代表系统矩阵，未知变量，观测数据；设d＝d⁺+ζ，其中d⁺表示无噪声数据，ζ表示噪声,ζ＝[ζ₁,ζ₂,···ζ_n]^T，ζ的每个分量ζ_i是独立同分布且服从期望为0，方差为的拉普拉斯分布，记为其概率密度函数表示为：

步骤(2)将噪声ζ的概率密度函数表示为高斯尺度的混合，具体形式为：

其中，ζ_i服从期望为0、方差为w_i的高斯分布，记为w_i是隐变量，w_i服从指数分布w_i～Exponential(τ)，其概率密度函数为w＝[w₁,w₂,···w_n]^T，w的概率密度函数

3.根据权利要求1所述的一种含拉普拉斯噪声的偏微分方程反问题的推断方法，其特征在于，所述步骤二中获取似然函数包括：

在考虑加性噪声模型时，似然函数是噪声ζ概率密度函数的平移，可表示为：即已知m和w时，d的概率密度函数，其中，m是未知变量，d是观测数据，w是服从指数分布的随机变量；d服从期望为Hm、协方差矩阵为diag(w₁,w₂,···w_n)的n维高斯分布，w_i是隐变量，w_i服从指数分布w_i～Exponential(τ)。

4.根据权利要求1所述的一种含拉普拉斯噪声的偏微分方程反问题的推断方法，其特征在于，所述步骤三中利用马尔科夫随机场表示未知变量的先验分布，包括：

采用先验建模工具马尔科夫随机场表示未知变量m的先验分布p(m|λ)，其形式为：其中，λ是描述交互作用强度的尺度参数，矩阵L表征相邻点之间交互作用的结构，s为矩阵L的秩；

所述步骤四中利用多层建模把尺度参数建模为随机变量，包括：选取伽马分布G(λ；α,β)为尺度参数λ的先验分布：其中，Γ(·)为标准伽马函数，α和β是非负常数。

5.根据权利要求1所述的一种含拉普拉斯噪声的偏微分方程反问题的推断方法，其特征在于，所述步骤五中通过拉普拉斯分布的高斯尺度混合表示将反问题的后验概率密度函数重新建模为多层贝叶斯模型，包括：利用多层贝叶斯建模将后验概率密度函数表示为：

其中，w是噪声方差参数，W是对角元素为w的对角矩阵，范数定义为(α,β)为伽马分布的参数对，m是未知变量，d是观测数据，λ是描述交互作用强度的尺度参数，矩阵L表征相邻点之间交互作用的结构，s为矩阵L的秩，α和β是非负常数，d服从期望为Hm、协方差矩阵为diag(w₁,w₂,···w_n)的n维高斯分布，w_i是隐变量，w_i服从指数分布w_i～Exponential(τ)；

所述步骤六中基于Kullback-Leibler散度，采用变分贝叶斯推断方法将后验概率密度函数的计算问题转化为关于概率分布的变分问题，包括：

将后验概率密度函数转化为求解优化问题在易计算的概率分布族中寻找近似分布q(m,w,λ)以逼近真实的后验概率密度函数p(m,w,λ|d)。

6.根据权利要求1所述的一种含拉普拉斯噪声的偏微分方程反问题的推断方法，其特征在于，所述步骤八中选取适当的逼近概率分布族，得到能够同时刻画未知变量、噪声方差参数和尺度参数的不确定性的非退化方法，包括：对未知变量m、噪声方差参数w和尺度参数λ施加平均场假设，近似分布q(m,w,λ)可分离为q(m)q(w)q(λ)，平均场假设同时刻画未知变量m、噪声方差参数w和尺度参数λ的不确定性，此条件下的近似称为非退化近似；

步骤八中选取适当的逼近概率分布族，得到能够仅刻画噪声方差参数w和尺度参数λ的不确定性的退化方法，包括：增加额外假设以降低计算复杂性，q(m)退化的情况下，即q(m)＝1时，仅能刻画噪声方差参数w和尺度参数λ不确定性，此条件下的近似称为退化近似。

7.根据权利要求6所述的一种含拉普拉斯噪声的偏微分方程反问题的推断方法，其特征在于，非退化近似，包括以下步骤：

(1)计算近似分布的解为q(m,w,λ)＝q^*(m)q^*(w)q^*(λ)，其中

(2)计算后验概率密度函数的对数函数ln p(m,w,λ|d)：

(3)反演未知变量m：

将(2)中公式代入(1)中第一个式子有：

其中，X(w,λ)是与m无关的项；

令则(3)前两项简化为：

再令m^*＝(H^TH+η^*L^TL)^-1H^Td，则上式为：

其中，最后一项对归一化常数有贡献，m服从均值为m^*，方差为(H^TW^*H+λ^*L^TL)^-1的高斯分布，

即：

(4)反演噪声方差参数w：

将(2)中公式代入(1)中第二个式子有：

其中，T(m,λ)是lnp(m,w,λ)中与w无关的项；

作变换z服从均值为μ_z，方差为λ_z的反高斯分布，即：其中，

(5)反演尺度参数λ：将(2)中公式代入(1)中第三个式子有：

其中，F(m,w)是与λ无关的项，λ服从均值为协方差为的伽马分布，即：

(6)对参数τ进行假设，参数τ与隐变量w_i直接相关，具体应用时仔细调整，以获得合理的结果；

(7)通过步骤(3)～步骤(5)参数反演，得到近似分布q(m,w,λ)的各变量边缘分布的非显式解析解；在分别求解q(m)、q(w)和q(λ)时，采用交替迭代法解决，即首先初始化q(m)、q(w)、q(λ)和q(m)、q(w)、q(λ)中涉及到的参数，然后从q(m)、q(w)、q(λ)中其中一个以及涉及到的参数开始交替迭代计算，直到算法收敛；交替迭代型算法对q(m)、q(w)和q(λ)估计；

其中，w是服从指数分布的随机变量；d服从期望为Hm、协方差矩阵为diag(w₁,w₂,···w_n)的n维高斯分布，w_i是隐变量，w_i服从指数分布w_i～Exponential(τ)，W是对角元素为w的对角矩阵，范数定义为(α,β)为伽马分布的参数对，m是未知变量，d是观测数据，λ是描述交互作用强度的尺度参数，矩阵L表征相邻点之间交互作用的结构，s为矩阵L的秩，α和β是非负常数。

8.根据权利要求6所述的一种含拉普拉斯噪声的偏微分方程反问题的推断方法，其特征在于，退化近似，包括以下步骤：

(1)进一步降低计算复杂性,在对m、w和λ施加条件独立性假设的基础上增加q(m)退化的假设，即q(m)＝1，对于固定的m，考虑最小化KL散度可以转化为公式：

引入KL散度的拉格朗日函数定义如下所示：

通过计算，求解得到：

m^*＝(H^TH+λ^*L^TL)^-1H^Td

类似非退化近似中的计算，得到：

(2)具体估计噪声方差参数w，

令得到z服从均值为μ_z，方差为λ_z的反高斯分布，即

(3)具体估计尺度参数λ，

λ服从均值为协方差为的伽马分布，

即其中m^*＝(H^TH+λ^*L^TL)^-1H^Td；

(4)采用交替迭代型算法对q(w)和q(λ)进行退化近似估计；

其中，w是服从指数分布的随机变量；d服从期望为Hm、协方差矩阵为diag(w₁,w₂,···w_n)的n维高斯分布，w_i是隐变量，w_i服从指数分布w_i～Exponential(τ)，W是对角元素为w的对角矩阵，范数定义为(α,β)为伽马分布的参数对，m是未知变量，d是观测数据，λ是描述交互作用强度的尺度参数，矩阵L表征相邻点之间交互作用的结构，s为矩阵L的秩，α和β是非负常数。