扩散方程反问题的正则化方法
Regularization Method for Inverse Problem of Diffusion Equation
DOI: 10.12677/pm.2024.144115, PDF, HTML, XML, 下载: 39  浏览: 111 
作者: 王清艳:兰州交通大学数理学院,甘肃 兰州
关键词: 扩散方程反问题辐射系数磨光法正则化Diffusion Equation Inverse Problem Radiation Coefficient Mollification Method Regularization
摘要: 研究了从终端观测数据重构扩散方程中辐射系数问题的一种正则化方法。对于带有噪声数据的终端观测值,运用磨光化的正则化方法,得出重构问题近似解的误差估计以及收敛速率。
Abstract: A regularization method for reconstructing the radiation coefficient problem in the diffusion equation from the terminal observation data is studied. For the terminal observations with noise data, the error estimation and convergence rate of the approximate solution of the reconstruction problem are obtained by using the polished regularization method.
文章引用:王清艳. 扩散方程反问题的正则化方法[J]. 理论数学, 2024, 14(4): 98-106. https://doi.org/10.12677/pm.2024.144115

1. 引言

扩散方程反问题在环境保护、生物医学、材料科学和化工等多个领域中具有重要的应用价值,可以用于污染物的追踪、地下水污染评估、药物输送模拟等领域的研究。通过反问题的求解,可以更准确地预测和控制扩散过程,优化扩散模型,从而制定理想的方案,解决物理问题,反问题的解决需要观测数据作为附加条件,但是由于观测数据有限、噪声干扰和不完全信息等因素的存在,反问题往往会面临不稳定性和非唯一性等困难。这意味着一个小的扰动或噪声可能导致问题解的巨大变化,或者存在多个解都能够很好地拟合观测数据。反问题正则化的主要目的就是抑制解的不稳定性和非唯一性,从而使得解更加稳定有效。通过合适的正则化选择和参数调整,得出更加合理有效的问题解,降低不确定性,并提供对解的有效性评估。因此,反问题正则化在物理、地球科学、医学成像等领域的问题研究中都有广泛的应用。

扩散方程反问题包括扩散系数,辐射系数等系数的识别研究。其中,辐射系数是描述介质对辐射传输的影响的重要参数,对于许多领域的研究和应用具有重要意义。然而,由于测量数据的不完备性和噪声的存在,直接从实验数据中准确地重构辐射系数是一项具有挑战性的任务。为了解决这个问题,许多正则化方法被提出并成功应用于辐射系数的重构。Deng在文献 [1] 中利用Tikhonov正则化方法研究了基于终端观测数据重构热传导方程

u t Δ u + q ( x ) u = 0

中辐射系数 q ( x ) 的反问题,得到了稳定的数值解。在文献 [2] 中研究了发展型反问题

u t Δ u + q ( x , t ) u = 0

中与时间和空间均相关的辐射系数 q ( x , t ) 的重构。Zhang等 [3] 根据扩散模型构造了一个单调算子重构辐射系数,并且证明了重构问题解的唯一性、稳定性,运用Galerkin有限元方法进行数值重构。Ivanova [4] 利用最小二乘方法研究了非线性热方程中热导率的识别问题。Yuksek [5] 利用变分法研究了扩散方程中空间相关扩散系数识别的反问题,并证明了反问题解的唯一性、稳定性。除此之外,还有贝叶斯方法 [6] ,Fourier正则化方法 [7] ,Landweber正则化方法 [8] 也可以解决不适定的反问题。

本文主要在Zhang等 [3] 研究出了分数阶扩散方程中辐射系数重构问题解的唯一性与稳定性的基础之上,采用磨光化 [9] 的正则化方法,解决重构辐射系数问题中的不适定性。其基本思想是对不准确的观测数据光滑化,化为一个近似光滑函数的求导问题,这种光滑化是通过与光滑核函数的卷积来实现的。这种方法不仅适用于Hilbert空间的问题,也适用于Banach空间中的问题。相对于Tikhonov正则化方法涉及的额外罚项、Fourier正则化方法需要考虑的周期性的变换、Landweber迭代正则化方法收敛速度较慢的问题,磨光化正则化方法直接作用于观测数据,对其进行磨光后再求导得到近似导数,保证其稳定性。

考虑一个有界域 D = [ 0 , L ] ,从初始时间0到终端时刻 T > 0 的时间区间内的扩散过程 α ( 0 , 1 ]

{ t α u ( x , t ) u x x ( x , t ) + q ( x ) u ( x , t ) = f ( x ) , ( x , t ) D × ( 0 , T ] , u ( 0 , t ) = b ( 0 ) , u ( L , t ) = b ( L ) , t ( 0 , T ] , u ( x , 0 ) = v ( x ) , x Ω , (1)

其中 f ( x ) 是与空间相关的源项, v ( x ) 是初始数据,函数 q ( x ) 是指标准抛物线情况下的辐射系数或反应系数或势函数。自始至终,假设辐射系数 q ( x ) 是空间相关的。

假设给出终端观测值作为附加条件

u ( x , T ) = g ( x ) , x D .

在实际测量中会不可避免地带有噪声影响,因此,取噪声水平 δ ,考虑如下形式的观测值

g δ ( x ) = g ( x ) [ 1 + δ × r a n d o m ( x ) ] (2)

2. 预备知识

根据扩散模型(1),定义一个算子K:

K q ( x ) = f ( x ) t α u ( x , T ; q ) + g ( x ) g ( x ) , q Q (3)

假设精确的辐射系数满足

q D ( K ) = { q C ( D ¯ ) : 0 q f ( x ) + Δ g ( x ) g ( x ) }

在文献 [3] 中,证明了算子K的单调性,以及迭代序列 q n = K q n 1 , n + 的收敛性,给出了重构势函数问题的唯一性与稳定性。本文在此基础上运用一种正则化方法克服反问题的不适定性。

假设1 假设问题(1)中的初始数据v,边界数据b以及源项f满足下列条件:

v H 2 ( D ) , v N > 0 , x D ,且对所有的 x D ,有 v ( x ) = b ( x )

b H 2 ( D ) , b N > 0 , x D

f W 1 , p ( D ) C ( D ¯ ) ( p > 2 ) , f 0

引理1 假设 q Q ,且函数 f , b , v 满足假设1中的条件,则问题(1)的解 u ( x , t ; q ) 满足以下性质

① 对所有 t > 0 u ( x , t ; q ) H 2 ( D ) ,存在一个与q无关的常数C,使得

max t [ 0 , T ] u ( x , t ; q ) L ( D ) C

t α u ( x , t ; q ) H 2 ( D ) H 0 1 ( D ) Δ u ( x , t ; q ) C ( D ¯ ) ,且对所有的 ( x , t ) D ¯ × [ 0 , T ]

t α u ( x , t ; q ) 0 , u ( x , t ; q ) M

引理2 假设 q 1 , q 2 Q ,且假设1成立,则对 t > t 0 以及任意正参数 ε < min ( 1 , 2 d 2 ) , d = 1 , 2 , 3 ,有

t α ( u ( q 1 ) u t ( q 2 ) ) ( t ) H 2 ( Ω ) c max ( t α , t ( 1 ε ) α ) q 2 q 1 L 2 ( Ω )

其中c与 q 1 , q 2 和t无关。

3. 正则化方法

首先,定义一个Gaussian核函数

ρ ε ( x ) = 1 ε π exp ( x 2 ε 2 ) C ( ) (4)

其中, ε > 0 表示平均滤波器的磨光半径或正则化参数。

然后,利用卷积的方法对噪声数据进行磨光化,根据上述核函数可以定义一个磨光算子 J ε

( J ε g ) ( x ) = ( ρ ε g ) ( x ) = ρ ε ( x s ) g ( s ) d s x a x + a ρ ε ( x s ) g ( s ) d s

它是 上的 函数,且 J ε g C [ a , L + a ] ,满足

J ε g C [ a , L + a ] g C ( D ¯ ) (5)

假设终端观测数据 g ( x ) , g δ ( x ) C 2 ( D ¯ ) g ( x ) D ¯ = [ 0 , L ] Lipschitz连续,且

g ( x ) g δ ( x ) C ( D ¯ ) δ

为了使正则化方法能够处理整个区间 D ¯ = [ 0 , L ] 上点的求导问题,先将函数 g ( x ) , g δ ( x ) 的定义域向区间两边作小的延拓。记 a = 3 ε , D ε = [ a , L + a ] ,对任意 g ( x ) C ( D ¯ ) ,可以在整个区间 D ε 上定义一个光滑的延拓函数 g ˜ ( x ) ,使得 g ˜ ( x ) , g ˜ δ ( x ) C 2 [ α , L + α ] ,且

g ˜ ( x ) g ˜ δ ( x ) C [ a , L + a ] C g ( x ) g δ ( x ) C ( D ¯ ) C δ (6)

延拓函数 g ˜ ( x ) 可以定义为

g ˜ ( x ) = { 0 , x ( , a ] [ L + a , ) , g ( 0 ) exp [ x 2 / ( a 2 x 2 ) ] , x ( a , 0 ) , g ( x ) , x D , g ( L ) exp [ ( x L ) 2 / ( ( x L ) 2 a 2 ) ] , x [ L , L + a ] .

g ( x ) 延拓到 上, g δ ( x ) 可以同样的方法延拓。

注意到,磨光算子 J ε 总是为正,并且在以原点为中心和半径为 3 ε 的区间外为零,即对 | x | 3 ε ρ ε ( x ) = 0 且满足 ρ ε ( x ) d x = 1 。从而噪声数据 g δ ( x ) 的导数在区间 [ 3 ε , L + 3 ε ] 上有比较好的结果。虽然 g δ ( x ) 是非光滑的,但其磨光后的函数 J ε ( g δ ) ( x ) 是一个 函数,因此可微。磨光后的导数可以由下式计算

在磨光算子 J ε 的作用下,可以用磨光化的数据来近似噪声数据,也即可以用 J ε ( g δ ) ( x ) 代替 g δ ( x )

接下来,给出正则化的终端时刻数据以及迭代算子的定义

g δ , a ( x ) = J ε ( g ˜ δ ) ( x ) (7)

K δ , a q = g δ , a t α u ( x , T ; q ) + f ( x ) g δ , a ( x ) (8)

迭代序列

q δ , a n = K δ , a q δ , a n 1 q δ , a n = K δ , a n q δ , a 0 (9)

其中,迭代 K δ , α n q δ , α 0 的迭代初值是

q δ , a 0 = ( g δ , a + f ) / g δ , a (10)

通过不断迭代趋于 q δ , α n ,接近于精确的辐射系数 q

4. 收敛性分析

首先给出上述延拓和磨光的一些性质。

引理3 设 g ( x ) D ¯ Lipschitz连续,则

g J ε g ˜ L ( D ) C a

其中,C只与gLipschitz常数有关。

证明 对 x D ,由表达式

J ε g ˜ ( x ) = ( ρ ε g ˜ ) ( x ) = x a x + a ρ ε ( x s ) g ˜ ( s ) d s

g ( x ) = x a x + a ρ ε ( x s ) g ( x ) d s

以及Lipschitz连续,有

| g ( x ) J ε g ˜ ( x ) | = | x a x + a ρ ε ( x s ) ( g ( x ) g ˜ ( s ) ) d s | C x a x + a ρ ε ( x s ) | x s | d s = C a

引理证毕。

引理4 设 g ( x ) , g δ ( x ) C 2 ( D ¯ ) g ( x ) D ¯ Lipschitz连续,且 g ( x ) g δ ( x ) C ( D ¯ ) δ 。则下式成立

g g δ , a L ( D ) C ε + C δ ε π

其中 g δ , a = J ε ( J ε ( g ˜ δ ) ) C ε δ 无关。

证明 易知 g δ , a = ( J ε g ˜ δ ) = J ε ( g ˜ δ )

g g δ , a L ( D ) g J ε g ˜ L ( D ) + J ε g ˜ J ε g ˜ δ L ( D )

x D ,有

J ε g ˜ ( x ) g ( x ) = ρ ε ( x s ) [ g ˜ ( s ) g ( x ) ] d s

进而可得

| J ε g ˜ ( x ) g ( x ) | ρ ε ( x s ) C | s x | d s x a x + a ρ ε ( x s ) C | s x | d s C α x a x + a ρ ε ( x s ) d s C ε

另一方面,有估计

| J ε g ˜ ( x ) J ε g ˜ δ ( x ) | = | ( J ε ( g ˜ g ˜ δ ) ) ( x ) | = ρ ε ( x s ) ( g ˜ ( s ) g ˜ δ ( s ) ) d s | ρ ε ( x s ) | | g ˜ ( s ) g ˜ δ ( s ) | d s C δ | ρ ε ( x s ) | d s = 2 C δ 0 ρ ε ( x ) d x = C δ ε π

综上得

g g δ , a L ( D ) C ε + C δ ε π

引理证毕。

引理5 设 g ( x ) , g δ ( x ) C 2 ( D ¯ ) D ¯ Lipschitz连续,且 g ( x ) g δ ( x ) C ( D ¯ ) δ 。则下式成立

g g δ , a L ( D ) C ε + C δ ε 2 π

其中C与 ε δ 无关。

证明 易知 g δ , a = ( J ε g ˜ δ ) = J ε ( g ˜ δ )

g g δ , a L ( D ) = g J ε g ˜ L ( D ) + J ε g ˜ J ε g ˜ δ L ( D )

首先,根据 g ( x ) D ¯ Lipschitz连续有

| g J ε g ˜ | ρ ε ( x s ) ( g ( x ) g ˜ ( s ) ) d s x a x + a ρ ε ( x s ) C | x s | d s a C

x D

I ( x ) : = ( J ε g ˜ ) ( x ) ( J ε g ˜ δ ) ( x ) = ( J ε ( g ˜ g ˜ δ ) ) ( x ) = x a x + a ρ ε ( x s ) ( g ˜ ( s ) g ˜ δ ( s ) ) d s

分部积分

I ( x ) = ρ ε ( a ) [ g ˜ ( x + a ) g ˜ δ ( x + a ) ] ρ ε ( a ) [ g ˜ ( x a ) g ˜ δ ( x a ) ] + x a x + a ρ ε ( x s ) [ g ˜ δ ( s ) g ˜ ( s ) ] d s

又已知 ρ ε ( a ) = ρ ε ( a ) = 0 ,从而

I ( x ) = x a x + a ρ ε ( x s ) [ g ˜ δ ( s ) g ˜ ( s ) ] d s

| I ( x ) | x a x + a | 2 ( x s ) ε 3 π | e ( x s ) 2 ε 2 [ g ˜ δ ( s ) g ˜ ( s ) ] d s 2 α ε 2 x a x + a 1 ε π e ( x s ) 2 ε 2 [ g ˜ δ ( s ) g ˜ ( s ) ] d s = 2 α ε 2 x a x + a ρ ε ( x s ) [ g ˜ δ ( s ) g ˜ ( s ) ] d s C δ 2 a ε 2 x a x + a | ρ ε ( x s ) | d s = C δ 2 a ε 2 2 0 a ρ ε ( x ) d x = C δ ε 2 π

综上有

g g δ , a L ( D ) C ε + C δ ε 2 π

引理证毕。

以下两个引理为 q q δ , a L 2 ( D ) 的误差提供了一个界限。

引理6 设 q , q δ , a 如(3)和(8)式定义,则对足够大的T有下式成立

q q δ , a L 2 ( D ) E 1 C max ( T α , T ( 1 ε ) α )

其中

E : = g g g δ , a g δ , a L 2 ( D ) + 1 g 1 g δ , a L 2 ( D ) [ t α u ( x , T ; 0 ) L ( D ) + f L ( D ) ] (11)

证明 因为 q 0 ( x ) > 0 ,根据文献 [3] 中引理2.5的证明有,对任意的

q 1 q 2 t α u ( x , t ; q 1 ) t α u ( x , t ; q 2 )

从而有

t α u ( x , t ; 0 ) t α u ( x , t ; q 0 ) 0

由初值 q 0 = ( g + f ) / g q δ , a 0 = ( g δ , a + f ) / g δ , a ,以及三角不等式有

q 0 q δ , a 0 L 2 ( D ) = g g g δ , a g δ , a L 2 ( D ) + 1 g 1 g δ , a L 2 ( D ) f L ( D ) E

然后对势函数的第一次迭代值有

q 1 q δ , α 1 L 2 ( D ) g g g δ , a g δ , a L 2 ( D ) + t α u ( x , T ; q 0 ) + f g t α u ( x , T ; q δ , a 0 ) + f g δ , a L 2 ( D ) g g g δ , a g δ , a L 2 ( D ) + 1 g 1 g δ , a L 2 ( D ) f L ( D ) + 1 g 1 g δ , a L 2 ( D ) t α u ( x , T ; q 0 ) L ( D ) + 1 g δ , a [ t α u ( x , T ; q 0 ) t α u ( x , T ; q δ , a 0 ) ] L 2 ( D ) E + 1 g δ , a [ t α u ( x , T ; q 0 ) t α u ( x , T ; q δ , a 0 ) ] L 2 ( D )

又已知 g δ , a 有界,根据引理2

1 g δ , a [ t α u ( x , T ; q 0 ) t α u ( x , T ; q δ , a 0 ) ] L 2 ( D ) C max ( T α , T ( 1 ε ) α ) q 0 q δ , a 0 L 2 ( D )

从而

q 1 q δ , a 1 L 2 ( D ) ( 1 + C max ( T α , T ( 1 ε ) α ) ) E

同理可得

q 2 q δ , a 2 L 2 ( D ) E + 1 g δ , a [ t α u ( x , T ; q 1 ) t α u ( x , T ; q δ , a 1 ) ] L 2 ( D ) E + C max ( T α , T ( 1 ε ) α ) q 1 q δ , a 1 L 2 ( D ) E + C max ( T α , T ( 1 ε ) α ) [ ( 1 + C max ( T α , T ( 1 ε ) α ) ) E ] ( 1 + C max ( T α , T ( 1 ε ) α ) + ( C max ( T α , T ( 1 ε ) α ) ) 2 ) E

继续这一论证可得

q q δ , a L 2 ( D ) lim n q n q δ , a n L 2 ( D ) lim n ( k = 0 n C k ( max ( T α , T ( 1 ε ) α ) ) k ) E E 1 C max ( T α , T ( 1 ε ) α )

引理证毕。

引理7 设 g ( x ) , g δ ( x ) C ( D ¯ ) g ( x ) D ¯ Lipschitz连续,且 g ( x ) g δ ( x ) C ( D ¯ ) δ 。则由(11)式定义的E,满足下式

E C ( ε + δ ) + C δ ε π + C δ ε 2 π

证明 已知

E = g g g δ , a g δ , a + 1 g 1 g δ , a [ t α u ( x , T ; 0 ) L ( D ) + f L ( D ) ]

对上式第一项有

g g g δ , a g δ , a L 2 ( D ) g ( g δ , a g ) g g δ , a L 2 ( D ) + g g δ , a g δ , a L 2 ( D )

根据引理1知, g ( x ) 有界,从而

E C [ g g δ , a L 2 ( D ) + g g δ , a L 2 ( D ) ]

其中

g g δ , a L ( D ) g J ε g ˜ L ( D ) + J ε ( g ˜ J ε ( g ˜ ) ) L ( D ) + J ε ( J ε ( g ˜ g ˜ δ ) ) L ( D )

根据(5)式有

g g δ , a L ( D ) C ( ε + δ )

根据引理4和引理5有

E C ( ε + δ ) + C δ ε π + C δ ε 2 π

引理证毕。

定理1 设 g ( x ) , g δ ( x ) C 2 ( D ¯ ) g ( x ) D ¯ Lipschitz连续,且 g ( x ) g δ ( x ) C ( D ¯ ) δ 。则对足够大的T,有下式成立

q q δ , a L 2 ( D ) C 1 C max ( T α , T ( 1 ε ) α ) ( ( ε + δ ) + δ ε π + δ ε 2 π )

这一定理是上述引理的直接结果。

推论1 设 g ( x ) , g δ ( x ) C 2 ( D ¯ ) g ( x ) D ¯ Lipschitz连续,且 g ( x ) g δ ( x ) C ( D ¯ ) δ 。T足够大,则选取 ε = δ 1 / 3 ,可以得到最优误差

q q δ , a L 2 ( D ) C 1 C max ( T α , T ( 1 ε ) α ) δ 1 / 3

5. 总结

文中研究了不适定问题的磨光化正则化方法,并且通过一系列的性质研究,给出了重构辐射系数问题的误差估计,以及 O ( δ 1 / 3 ) 的收敛速度。本文是在前人研究的理论基础之上,提出的解决线性扩散系统不适定问题的正则化方法,在后续的研究中可以考虑磨光化正则化在非线性系统中的应用。

参考文献

[1] Deng, Z.C., Yang, L. and Yu, J.N. (2009) Identifying the Radiative Coefficient of Heat Conduction Equations from Discrete Measurement Data. Applied Mathematics Letters, 22, 495-500.
https://doi.org/10.1016/j.aml.2008.06.023
[2] Deng, Z.C., Yang, L., Yu, J.N., et al. (2010) Identifying the Radiative Coefficient of an Evolutional Type Heat Conduction Equation by Optimization Method. Journal of Mathematical Analysis and Applications, 362, 210-223.
https://doi.org/10.1016/j.jmaa.2009.08.042
[3] Zhang, Z., Zhang, Z. and Zhou, Z. (2022) Identification of Potential in Diffusion Equations from Terminal Observation: Analysis and Discrete Approximation. SIAM Journal on Numerical Analysis, 60, 2834-2865.
https://doi.org/10.1137/21M1446708
[4] Ivanova, A., Migorski, S., Wyczolkowski, R., et al. (2020) Numerical Identification of Temperature Dependent Thermal Conductivity Using Least Squares Method. International Journal of Numerical Methods for Heat & Fluid Flow, 30, 3083-3099.
https://doi.org/10.1108/HFF-10-2018-0589
[5] Yüksek, K., Koca, Y. and Sadikoglu, H. (2009) Solution of Counter Diffusion Problem with Position Dependent Diffusion Coefficent by Using Variational Methods. Journal of Computational and Applied Mathematics, 232, 285-294.
https://doi.org/10.1016/j.cam.2009.06.009
[6] Latz, J. (2023) Bayesian Inverse Problems Are Usually Well-Posed. SIAM Review, 65, 831-865.
https://doi.org/10.1137/23M1556435
[7] Kerimov, N.B. and Ismailov, M.I. (2012) An Inverse Coefficient Problem for the Heat Equation in the Case of Nonlocal Boundary Conditions. Journal of Mathematical Analysis and Applications, 396, 546-554.
https://doi.org/10.1016/j.jmaa.2012.06.046
[8] Yang, F., Wang, N. and Li, X.X. (2020) Landweber Iterative Method for an Inverse Source Problem of Time-Fractional Diffusion-Wave Equation on Spherically Symmetric Domain. Journal of Applied Analysis & Computation, 10, 514-529.
https://doi.org/10.11948/20180279
[9] 刘继军. 不适定问题的正则化方法及应用[M]. 北京: 科学出版社, 2005.