近几十年以来, 由于科学技术的突飞猛进, 实验室里越来越多地发现自然界存在众多非局部扩散现象. 基于分数阶导数自身具有的记忆性与遗传性等非局部特征, 分数阶模型被广泛地应用于科学与工程等诸多关键问题, 特别是在处理具有分数阶特性的信号问题时, 分数阶数值微分已成为处理该类问题的一项重要工具, 被常用于如生物医学、地震及卫星等信号的去噪、特征提取等问题.

由于分数阶数值微分有着重要的实际应用价值, 而数值微分 (整数阶/分数阶) 为典型的不适定问题^[1,2], 需要采取正则化方法进行计算, 目前已有众多研究者围绕分数阶数值微分正则化方法开展了数值和理论上的广泛研究. Murio^[3] 利用磨光正则化方法研究了扰动数据的分数阶数值微分, 且推导了相应磨光正则化数值微分的收敛阶; Zhang 等^[4]基于对扰动数据的磨光, 利用 Jacobi 谱方法对磨光数据进行分数阶数值微分研究, 在精确数据满足一定光滑性假设的前提下推得了相应正则化数值微分的收敛率. Qian 等^[5]利用傅里叶谱截断正则化方法研究了扰动数据的分数阶数值微分, 且导出了该正则分数阶数值微分 Hölder 型的收敛率结果; Li 等^[6]将 Caputo 导数离散化为加权和形式, 通过显式权重设计和依据噪声水平自适应选择步长的方式有效正则化了 Caputo 导数计算中的奇异性与不适定性问题, 且平衡了计算的复杂度和精度. Liu 等^[7,8,9]分别基于雅可比正交多项式函数逼近技术和 B- 样条函数逼近技术, 将分数阶导数计算问题转为分数次幂函数求和问题, 数值计算取得了较好的结果; Li 等^[10]基于对扰动数据的 Gaussian 卷积磨光化, 将扰动数据的数值微分运算转为对扰动数据与 Gaussian 核分数阶导数的卷积运算, 数值上取得了不错的计算效果. 最近 Duc 等^[11]基于对一阶导数的 Lavrentiev 正则化和分数阶导数定义, 构造了数值求解 Caputo 分数阶导数的一类正则化方法, 且在 Banach 空间 $L^{\infty}[T]$ 下分析了正则化分数阶导数的收敛率; 在文献[12] 中, Duc 等基于 Tikhonov 正则化方法, 构造了正则化函数空间和优化泛函计算 $L^{2}[T]$ 空间中的 Caputo 分数阶导数, 在理论和数值上均表现出良好的稳定性.

以上分数阶数值微分问题主要研究的分数阶为 0$\sim$1 阶, 基于现实驱动, 0$\sim$1 阶 Caputo 分数阶导数不足以描述一些复杂的物理现象. 因此本文考虑扰动数据 1$\sim$2 阶分数阶导数的正则化计算方法及其收敛性分析. 本文结构如下: 第 2 节为问题描述, 并阐述了本文构造的正则化方法的基本思想; 第 3 节给出了一些定理, 为后续的研究做准备; 第 4 节基于正则化参数的先验选取规则, 给出了先验误差估计的收敛阶; 第 5 节给出了关于正则化参数的后验选取规则的主要结果; 第 6 节为数值算例; 第 7 节为结论.

为了更好地描述本文所研究的内容, 首先回顾下 Caputo 分数阶导数的基本定义.

定义 2.1^[13] 设 $\alpha>0$, 且 $n-1<\alpha<n$, 对于定义在区间 $[T]$ 上的 $n$ 次可微函数 $q(t)$, 则其 $\alpha$ 阶 Caputo 分数阶导数为

其中 $\Gamma$ 为 Gamma 函数, $q^{(n)}(t)$ 表示对 $q(t)$ 求 $n$ 阶导数.

基于现实考虑, 精确数据 $q(t)$ 是未知的, 往往只能获得有噪声的观测数据, 设噪声数据 $q^{\delta}\in L^{\infty}[T]$ 满足

其中 $\delta>0$ 表示噪声水平.

本文考虑的 1$\sim$2 阶 Caputo 分数阶数值微分问题归结为: 利用满足误差估计 (2.2) 的观测数据 $q^{\delta}(t)$ 近似计算 $q(t)$ 的分数阶导数 $\left(D^{(\alpha)}q\right)(t)$.

本文主要在文献[11] 的基础上, 基于 Lavrentiev 正则化思想, 构造计算 1 至 2 阶 Caputo 分数阶导数的正则化方法, 该方法构造的基本思路归结如下:

第 1 步 基于 Lavrentiev 正则化方法, 通过计算正则化方程 (2.3) 求解观测数据 $q^{\delta}(t)$ 的一阶导数的正则化解 $u_{\beta}(t)$,

其中 $\beta>0$ 是正则化参数.

第 2 步 继续基于 Lavrentiev 正则化方法, 通过计算正则化方程 (2.4) 求解观测数据 $q^{\delta}(t)$ 的二阶导数的正则化解 $\omega_{\beta}(t)$,

第 3 步 由分数阶导数定义 (2.1) 得到观测数据 $q^{\delta}(t)$ 的 $\alpha(1<\alpha<2)$ 阶分数阶导数的正则化解 $\omega_{\alpha,\beta}^{\delta}(t)$,

在讨论正则化解的收敛性之前, 先证明以下结果.

定理 3.1 设 $q(t)$ 是 $[T]$ 上三阶连续可微函数, 满足 $q(0)=q'(0)=q''(0)=0$, 且存在常数 $E>0$ 使得 $\|q^{(3)}\|_{\infty}\leq E$, 那么以下误差估计成立

证 对 (2.3) 和 (2.4}) 式使用求导与积分法则, 可知 $u_{\beta}(t)$、$\omega_{\beta}(t)$ 的显示表达式为

构造辅助函数 $v_{\beta}(t)$、$g_{\beta}(t)$ 满足以下方程

定义函数 $g_{\alpha, \beta}(t)$ 为

显然 $v_{\beta}(t)$、$g_{\beta}(t)$ 的显式表达式为

根据 (3.2) 和 (3.5) 式, 有

根据文献[11] 可知 $\left\|u_{\beta}-v_{\beta}\right\|_{\infty} \leq \dfrac{2-{\rm e}^{-T / \beta}}{\beta} \delta$, 则有

对 (3.4) 式进行分部积分, 可得

由于 $q(0)=q'(0)=q''(0)=0$, 那么

由于 $\|q^{(3)}\|_{\infty}\leq E$, 那么

所以

因此

由此证得定理 3.1.

从定理 3.1 可知, 正则化方法的稳定性与正则化参数 $\beta$ 的选取有关. 接下来考虑正则化参数先验与后验选取规则及其在两类选取规则下分数阶导数的收敛率.

本节主要讨论先验参数选取规则及在先验选取规则下正则化分数阶导数的收敛性.

定理 4.1 若 $q(t)$ 是 $[T]$ 上三阶连续可微函数, 满足 $q(0)=q'(0)=q''(0)=0$, 且存在常数 $E>0$ 使得$\left\|q^{(3)}\right\|_{\infty} \leq E$, 则以下估计成立:

特别地, 当 $\beta=\sqrt[3]{2 \delta / E}$ 时

证 根据定理 3.1 可知 $\left\|\omega_{\beta}-q''\right\|_{\infty} \leq \dfrac{\left(2-{\rm e}^{-T / \beta}\right)^{2}}{\beta^{2}} \delta+2 \beta E\left(1-{\rm e}^{-T / \beta}\right)$, 因此

即有

令 $\beta=\sqrt[3]{2 \delta / E}$, 代入即得

即得证.

定理 4.1 给出了 1$\sim$2 阶正则化导数的先验误差估计, 该定理可以进行一般推广, 推得任意阶分数阶导数正则化解的先验收敛率. 在推广之前, 先利用数学归纳法证明如下引理:

引理 4.1 若 $q(t)$ 是 $[T]$ 上 $n+1$ 阶连续可微函数, 满足$q(0)=q^{\prime}(0)=...=q^{(n)}(0)=0$, 且存在常数 $E>0$ 使得 $\left\|q^{(n+1)}\right\|_{\infty} \leq E$, 则有

其中函数 $g^{i}_{\beta}(t)$ 满足方程

且 $g^{0}_{\beta}(t)=q(t)$.

证 根据求导与积分法则可知 $g^{i}_{\beta}(t)$ 的显式表达式为

显然当 $i=1$ 时, 由文献[11] 可知

当 $i=2$ 时, 根据 (3.7) 式有

设 $i=n-1$ 时, 猜想 $g_{\beta}^{n-1}(t)$ 有如下表达式

当 $i=n$ 时, 将上式代入 (4.3}) 式, 即可得

因此, 根据数学归纳法可知猜想成立, 即有

由于 $i=n$ 时, $q(0)=q^{\prime}(0)=...=q^{(n)}(0)=0$, 且$\left\|q^{(n+1)}\right\|_{\infty} \leq E$, 故有

引理得证.

根据 $2.2$ 节正则化思路, 定义噪声数据 $q^{\delta}(t)$ 的 $n$ 阶导数的正则化导数 $\omega^{n}_{\beta}(t)$ 为下列方程的解:

其中 $\omega^{0}_{\beta}(t)=q^{\delta}(t)\text{.}$

定义 $q^{\delta}(t)$ 的 $\alpha(n-1<\alpha<n)$ 阶正则化导数为:

则有如下一般化的正则化分数阶导数的收敛结论.

定理 4.2 若 $q(t)$ 是 $[T]$ 上 $n+1$ 阶连续可微函数, 满足$q(0)=q^{\prime}(0)=...=q^{(n)}(0)=0$, 且存在常数 $E>0$ 使得$\left\|q^{(n+1)}\right\|_{\infty} \leq E$, 则以下估计成立:

证 根据求导与积分法则可知 $\omega^{i}_{\beta}(t)$ 的显式表达式为

故有

则有

由于 $\omega^{0}_{\beta}(t)=q^{\delta}(t),g^{0}_{\beta}(t)=q(t)$, 且 $\|q-q^{\delta}\|_{\infty}\leq\delta$, 故可知

则根据引理 4.1, 由三角不等式可知

由于

因此, 可知

当 $\beta=\sqrt[n+1]{(2^n\delta)/(nE)}$ 时, 则有

定理得证.

本节讨论正则化参数后验选取准则及正则化分数阶导数的后验收敛率. 在推导后验误差估计前, 先证明以下稳定性结果.

定理 5.1 (稳定性估计) 假设 $q_{j}(t), j=1,2$ 是 $L^{\infty}[T]$ 上三阶连续可微函数, 且满足

则以下估计成立

证 令 $h(t)=\dfrac{q_{1}(t)+q_{2}(t)}{2}$, 则

假设函数 $f_{\beta}(t), F_{\beta}(t)$ 满足

定义函数 $F_{\alpha, \beta}(t)$ 为

则由定理 4.1 可得

当 $\beta=\sqrt[3]{\delta / E}$ 时,

根据三角不等式可知

$\begin{equation} \left\|D^{(\alpha)} q_{1}-D^{(\alpha)} q_{2}\right\| \leq \dfrac{8 T^{2-\alpha}}{(2-\alpha) \Gamma(2-\alpha)} \sqrt[3]{\delta E^{2}}.\end{equation}$ (5.1)

由此定理 5.1 得证.

由于参数 $\beta$ 的选取与测量数据有关, 令 $\tau>8+\delta^{-2 / 3}$, 假定观测数据满足

取后验正则化参数 $\beta_{\delta}>0$ 满足下列偏差准则

满足偏差准则 (5.3) 的后验正则化参数 $\beta$ 的存在性由下面引理 5.1 给出.

引理 5.1 ($\beta$ 的存在性) 令 $\varphi(\beta)=\left\| 2 \beta u_{\beta}-\beta^{2} \omega_{\beta} \right\|_{\infty}, \beta>0$, 则

i) $\lim\limits_{\beta \rightarrow \infty}\left\|2 \beta u_{\beta}-\beta^{2} \omega_{\beta}\right\|_{\infty}=\left\|q^{\delta}\right\|_{\infty}$;

ii) 当 $\beta$、$\delta$ 足够小时, $\varphi(\beta)<\tau \delta$;

iii) $\varphi(\beta)$ 是 $(0, \infty)$ 上的连续函数.

证 i) 根据 (3.1) 和 (3.2) 式可得

则根据三角不等式可知

而当 $\beta \rightarrow \infty$ 时, 有

因此

ii) 由于

其中根据分部积分法则, 有

由于 $q(0)=q'(0)=q''(0)=0,\left\|q^{(3)}\right\|_{\infty} \leq E$, 且

所以

故对于 $\forall \beta \in\left(0, \dfrac{\delta}{6 E T^{2}}\right)$, 以下不等关系成立:

由于 $\tau>8+\delta^{-2/3}$, 则当 $\beta, \delta$ 足够小时, 有

iii) 由于 $2 \beta u_{\beta}-\beta^{2} \omega_{\beta}=-\dfrac{1}{\beta^{2}} {\rm e}^{-t / \beta} \int_{0}^{t} \int_{0}^{s} {\rm e}^{\tau / \beta} q^{\delta}(\tau) {\rm d} \tau {\rm d} s+q^{\delta}(t)$, 若令

则 $\phi(\beta)$ 在 $(0, \infty)$ 上是连续函数, 那么关于 $\beta$ 的函数 $\left\|2 \beta u_{\beta}-\beta^{2} \omega_{\beta}\right\|_{\infty}$ 也是连续函数, 故 $\varphi(\beta)$ 在 $(0, \infty)$ 上是连续函数.

综上, 引理证毕.

由引理 5.1 可知, 满足偏差准则 (5.3) 的后验正则化参数是存在的. 下面证明根据偏差准则选择的后验参数对应正则化分数阶导数的后验误差估计.

定理 5.2 若 $q(t)$ 满足定理 4.1 的条件, 且观测数据满足 (2.2), 正则化参数 $\beta_{\delta}$ 由偏差准则 (5.3) 决定, 则正则化分数阶导数 $\omega_{\alpha,\beta_{\delta}}$ 满足以下误差估计:

证 令 $G(t)=\int_{0}^{t} \int_{0}^{s} g_{\beta_{\delta}}(\tau) {\rm d} \tau {\rm d} s, t \in[T]$, 则

显然, $G(0)=G'(0)=G''(0)=0$, 且 $G^{(3)}(t)=g_{\beta_{\delta}}'(t)$, 结合 (3.7) 式则有

由于 $\left\|q^{(3)}\right\|_{\infty} \leq E$, 所以

根据 $v_{\beta}(t)$、$g_{\beta}(t)$ 的定义可知

故根据定理 3.1 可得

因此由定理 5.1 可知

此外, 由于

故当 $\tau>8+\delta^{-2 / 3}$ 时, 有

根据 (5.5) 式可知

且由拉格朗日中值定理可知, 存在 $\xi \in(0, t)$ 使得

根据 (3.6) 式对 $t$ 求导, 可得$v_{\beta_{\delta}}'(t)=q''(t)-{\rm e}^{-t / \beta_{\delta}} \int_{0}^{t} {\rm e}^{s / \beta_{\delta}} q^{(3)}(s){\rm d}s\text{,}$

因此, 当 $t=0$ 时, $v_{\beta_{\delta}}(0)=v_{\beta_{\delta}}'(0)=0$. 则有

故

根据 (5.7) 和 (5.9) 式, 可得

因此 $\beta_{\delta}>\dfrac{1}{2 E T^{2}} \delta^{\frac{1}{3}}$, 故可知

代入 Caputo 分数阶导数的定义, 可得

因此

结合 (5.6) 式, 由三角不等式可得

由此, 定理得证.

在实际计算中, 由于部分函数的 Caputo 分数阶导数不存在解析形式, 因此本文采用高斯—克朗罗德 (Gauss-Kronrod) 数值积分进行数值求解. 在分析由数值积分得到的正则化分数阶导数的收敛性之前, 给出以下引理.

引理 5.2^[14] 若$f(t)\in C\left[a,b\right],F(x)=\int_{a}^{x}f(t){\rm d}t$, 则对于 Gauss 型数值积分 $G_{n}(x_{i})$, 有

记 $w_{\alpha,\beta_{\delta}}^{i}$ 为 $q(t)$ 在 $t_{i}$ 处的正则化数值微分, 由引理 5.2 可知, 当求积节点 $n$ 足够多时, 通过数值积分求解的正则化数值微分 $w_{\alpha,\beta_{\delta}}^{i}$ 满足

因此, 利用高斯型积分得到的正则化数值微分具有以下收敛性结果.

定理 5.3 若正则化参数 $\beta_{\delta}$ 满足偏差准则 (5.3), 且 $q^\delta(t)\in C[T]$, 取 Gauss-Kronrod 数值积分的误差控制为 $O(\delta^{1/3})$, 则有

证 由定理 5.2 可知, 对于 $\forall t_{i}\in[T]\text{, }$ 有

故由三角不等式可得

本节通过两个数值算例来验证方法的有效性. 在本文的两个算例中, 精确数据与噪声数据的关系如下:

其中 $\delta$ 为噪声水平, $\operatorname{rand}(i)$ 表示 $(0,1)$ 上服从均匀分布的随机数. 在每个数值算例中, $T=3$, 正则化参数根据后验参数选取规则选择, 并采用 $l_{2}$ 范数计算相对误差.

例 1: 精确数据 $q(t)=-2 t+\sin 2 t,~ t \in[T]$. 不同噪声水平下的相对误差见表1.

例 2: 精确数据 $q(t)= \begin{cases} \quad t^4, &0\leq t\leq\frac{1}{2} \\ \\2t^3-\frac{3}{2}t^2+\frac{1}{2}t-\frac{1}{16},&t>\frac{1}{2} \end{cases}, t \in[T]$. 不同噪声水平下的相对误差见表2.

根据上述两个算例的反演情况, 可以发现误差基本小于噪声水平, 且当噪声水平下降时, 反演的正则化解更加接近于精确解. 由此说明本文提出的正则化方法对于求解扰动函数的 Caputo 数值微分是可行且有效的.

本文提出了一种正则化方法求解带有噪声函数的 1 至 2 阶 Caputo 分数阶导数. 理论上给出了分数阶导数的先验误差估计和后验误差估, 数值上通过数值积分规则实现了 Caputo 分数阶导数的正则化计算, 为分数阶数值微分问题提供理论与计算方面的参考.