数值微分, 即由给定函数的测量数据近似求其导数, 这类问题在科学研究及工程实践中有广泛的应用价值.如在数字图像处理中, 微分算子方法是常用的边缘检测算法^[1].图像的边缘, 即图像强度变化急剧的部分, 一般可通过寻找图像强度的一阶导数的极大值或二阶导数的交叉零点得到.又如磁共振电阻抗成像技术中, 介质的导电率是利用生物组织内部电流产生的磁场信息重建的.基于主磁场方向磁感应强度的调和 $B_z$算法是当前主流的重建算法, 而如何由磁感应强度的测量数据计算其二阶导数是该算法的关键步骤^[2].数值微分的具体应用还体现在天气预报、期权定价、放射性废弃物处理等方面.

数值微分是一类典型的不适定问题, 其求解的本质困难在于解的不稳定性, 即当测量数据有较小扰动时, 对其直接求导产生的误差可能是任意大的.为此, 如何构造稳定的微分算法一直是数值微分的研究重点.总的来说, 已有算法可分为有限差分法、正则化方法、平滑化方法等几类.有限差分法是直接利用差分格式构造待求导函数的近似导数, 该方法格式简单且应用广泛.当测量数据伴有误差时, 差分格式的稳定性一般依赖于差分步长的选取^[3-5].正则化方法是指借助正则化的思想, 将不适定的数值微分问题转化为适定的近似问题进行求解, 进而得到稳定近似导数的微分算法^[6-10].平滑化方法是将待求导函数进行平滑处理后得到一个光滑的近似函数, 对该函数求导可得稳定的近似导数^{[11, 12]}.该方法在数字图像处理及计算机视觉的诸多领域, 如图像去噪、复原、特征检测与匹配等, 中应用十分广泛.

除稳定性要求外, 数值微分算法还要保证近似导数逼近精确导数的精度.然而, 上述微分算法在保证稳定性的同时, 也减弱了近似导数逼近精确导数的效果.保证近似导数的稳定性和精度已成为数值微分算法相互矛盾的两个性能要求, 解决这一冲突的常规做法是引入一个折衷参数.适当选取折衷参数可在一定程度上保证两方面性能的最优化, 即近似导数能在有效抑制噪声的同时较好地逼近精确导数.在已有的算法中, 有限差分法中的差分步长、正则化方法中的正则化参数、平滑化方法中的尺度参数等都起到了这样的作用.然而, 受制于方法自身的缺陷, 折衷难免会导致某一方面性能的降低(这将在下一节中继续讨论).

为此, 在折衷思想不变的框架下, 找寻一种更适当的方法显得很有必要.在本文中, 我们将引入一种能有效利用函数局部信息反演其局部特性的方法---局部正则化方法^[13-15], 并将其用于构造稳定的数值微分算法.该微分算法能保证近似导数在有效抑制噪声的同时, 较好地反应精确导数有间断或急剧变化的部分.这一点我们将通过数值实验来说明.已有的局部正则化理论是在连续函数空间中建立的, 考虑到实际应用中待求导函数一般不具有足够的光滑性, 我们将其理论推广到更一般的可积函数空间中, 并给出正则化参数的先验选取策略及相应近似导数的误差估计.

众所周知, 微分和积分是两种互逆的基本运算, 求微分可等价为求积分的逆, 即解第一类的积分方程.引入正则化方法求解该积分方程, 可构造稳定的数值微分算法.考虑函数 $f\in H^1[0,1]$, $u=f'\in L^2[0,1]$为其一阶导数.不失一般性, 我们设 $f(0)=0$, 则 $f$和 $u$满足第一类的Volterra型积分方程

$ \begin{eqnarray} A[u](x):=\int_0^xu(s)ds=f(x), \quad x\in[0,1]. \end{eqnarray} $

(2.1)

易知, 算子 $A$为空间 $L^2[0,1]$上的有界线性紧算子.考虑有噪声的测量数据 $f^\delta\in L^2[0,1]$, 算子方程

$ \begin{eqnarray} A[u](x)=f^\delta(x) \end{eqnarray} $

(2.2)

的求解需引入正则化方法.经典的Tikhonov正则化方法可用于求解该方程, 相应的正则化解 $u^{\alpha, \delta}$满足

$ \begin{eqnarray} \alpha u^{\alpha, \delta}(x)+A^*A[u^{\alpha, \delta}](x)=A^*[f^\delta](x), \end{eqnarray} $

(2.3)

其中 $A^*[v](x):=\int_x^1v(s)ds$为算子 $A$的伴随算子.对任意的 $v\in L^2[0,1]$, 有

$ \begin{eqnarray} (A[v], v)=\int_0^1v(x)\int_0^xv(s)dsdx=\int_0^1\int_0^xv(s)dsd(\int_0^xv(s)ds)=\frac{1}{2}(\int_0^xv(s)ds)^2\geq0. \nonumber \end{eqnarray} $

这说明算子 $A$为非负算子.因此, 还可引入Lavrentiev正则化方法求解算子方程(2.1), 相应的正则化解 $v^{\alpha, \delta}$满足

$ \begin{eqnarray} \alpha v^{\alpha, \delta}(x)+A[v^{\alpha, \delta}](x)=f^\delta(x). \end{eqnarray} $

(2.4)

此时, 正则化解具有解析表达式

$ \begin{eqnarray} v^{\alpha, \delta}(x)=\frac{-1}{\alpha^2}\int_0^xe^{\frac{s-x}{\alpha}}f^\delta(s)ds+\frac{f^\delta(x)}{\alpha}, \quad x\in[0,1]. \end{eqnarray} $

(2.5)

众所周知, Tikhonov正则化方法求得的正则化解是过于光滑的^[16], 这也反应在方程(2.3) 右端的积分项上.区间 $[x, 1]$上函数 $f^\delta$的积分运算实际上是对其进行的一种平滑化处理.与之相反, Lavrentiev正则化方法求得的正则化解对噪声是过于敏感的.这一点可通过数值实验来说明, 如取 $ f(x)=|x-\frac{1}{2}|-\frac{1}{2}, x\in[0,1], $对其添加 $1\%$的高斯白噪声得 $f^\delta$.在图 1中, 我们给出了利用两种经典的正则化方法求得的正则化解 $u^{\alpha, \delta}$, $v^{\alpha, \delta}$以及精确导数 $u$的图像, 其中Tikhonov正则化的参数取为 $\alpha=0.005$, Lavrentiev正则化的参数取为 $\alpha=0.02$.值得注意的是, 正则化解 $u^{\alpha, \delta}$在右端点附近及 $v^{\alpha, \delta}$在左端点附近逼近精确导数的效果均不理想, 这是由正则化方法自身的缺陷所决定的^[16].从图 1中可以看出, $u^{\alpha, \delta}$足够光滑, 但不能反应 $x=0.5$处精确导数的间断情况; $v^{\alpha, \delta}$虽能反应精确导数的间断情况, 但对噪声过于敏感.为此, 我们将在下一节中引入求解算子方程(2.1) 的局部正则化方法.该方法能有效克服两类经典正则化方法的不足, 在有效抑制噪声的同时较好地反应精确导数有间断或急剧变化的部分.

设 $r\in(0, R]$, 其中 $R>0$是一个很小的正数, 为引入局部正则化方法, 我们将函数 $f$和 $u$的定义区间延伸至 $[0, 1+R]$上, 此时有

$ \int_0^{x+\rho}u(s)ds=f(x+\rho), \quad x\in[0,1], \quad \rho\in[0, r]. $

(3.1)

分解等式左边的积分项可得

$ \int_0^xu(s)ds+\int_x^{x+\rho}u(s)ds=f(x+\rho), \quad x\in[0,1], \quad \rho\in[0, r].\nonumber $

等式两边关于 $\rho$在区间 $[0, r]$上积分可得

$ \int_0^r\int_0^xu(s)dsd\rho+\int_0^r\int_x^{x+\rho}u(s)dsd\rho=\int_0^rf(x+\rho)d\rho, \quad x\in[0,1], \nonumber $

即

$ r\int_0^xu(s)ds+\int_0^r\int_0^{\rho}u(x+s)dsd\rho=\int_0^rf(x+\rho)d\rho, \quad x\in[0,1]. $

(3.2)

注意到(3.2) 式仍然是精确导数 $u$满足的方程.当我们将给定函数 $f$换成伴有噪声的测量数据 $f^\delta\in L^2[0,1]$时, 为保证近似导数的稳定性, 方程(3.2) 的求解需引入正则化方法.当 $R>0$足够小时, $u$在小区间 $[x, x+r]$上的取值可近似为与 $r$无关的 $u(x)$, 即上式左端第二项中的被积函数 $u(x+s)$可近似替换为 $u(x)$, 进而有正则化方程

$ r\int_0^xu(s)ds+\frac{r^2}{2}u(x)=\int_0^rf^\delta(x+\rho)d\rho, \quad x\in[0,1]. $

(3.3)

此时, 参数 $r\in(0, R]$起到了正则化参数的作用^[13].(3.3) 式两端同时除以 $r$可得

$ \begin{eqnarray} \int_0^xu(s)ds+\frac{r}{2}u(x)=f_r^\delta(x), \quad x\in[0,1], \end{eqnarray} $

(3.4)

其中 $f_r^\delta(x)=\frac{1}{r}\int_0^rf^\delta(x+\rho)d\rho$.由于 $r>0$, 算子 $A$为有界线性算子, 方程(3.4) 存在唯一解 $u^{r, \delta}\in L^2[0,1]$, 且该解连续依赖于右端数据 $f_r^\delta$^[14].不难说明, 此时的正则化解 $u^{r, \delta}$具有解析表达式

$ \begin{eqnarray} u^{r, \delta}(x)=\frac{-4}{r^3}\int_0^xe^{\frac{2(s-x)}{r}}\int_0^rf^\delta(s+\rho)d\rho ds+\frac{2}{r^2}\int_0^rf^\delta(x+\rho)d\rho, \quad x\in[0,1]. \end{eqnarray} $

(3.5)

在本节中, 我们将在可积函数空间中给出正则化参数的先验选取策略及相应近似导数的误差估计.为估计近似导数的误差并给出其收敛速度, 我们设精确导数 $u\in H^1[0, 1+R]$, 即 $u'\in L^2[0, 1+R]$.相应的, 有如下结论(如不特别说明, 文中出现的范数均指空间 $L^2[0,1]$上的 $L^2$范数).

定理1  设精确导数 $u\in H^1[0, 1+R]$, 测量数据 $f^\delta$满足 $\|f^\delta-f\|_{L^2[0, 1+R]}\leq\delta$, 则正则化方程(3.4) 的解 $u^{r, \delta}$满足

$ \begin{eqnarray} \|u^{r, \delta}-u\|\leq\frac{2\delta}{r}+Cr, \end{eqnarray} $

(4.1)

其中 $C$的取值与 $\delta$无关.进一步当取 $r(\delta)=\sqrt{2\delta/C}$时, 有 $\|u^{r, \delta}-u\|=O(\sqrt{\delta})$.

证  已知正则化解 $u^{r, \delta}$满足方程

$ r\int_0^xu^{r, \delta}(s)ds+\frac{r^2}{2}u^{r, \delta}(x)=\int_0^rf^\delta(x+\rho)d\rho, \quad x\in[0,1]. $

(4.2)

记 $e^{r, \delta}=u^{r, \delta}-u$, $(4.2)$式减去 $(3.2)$式可得

$ r\int_0^xe^{r, \delta}(s)ds+\frac{r^2}{2}e^{r, \delta}(x)=\int_0^r[f^\delta(x+\rho)-f(x+\rho)]d\rho+\int_0^r\int_0^\rho[u(x+s)-u(x)]dsd\rho, \nonumber $

其中 $x\in[0,1]$.故有

$ \begin{eqnarray} e^{r, \delta}(x)+\frac{2}{r}\int_0^xe^{r, \delta}(s)ds=F_r^\delta(x)+U_r(x), \quad x\in[0,1], \end{eqnarray} $

(4.3)

其中

$ F_r^\delta(x)=\frac{2}{r^2}\int_0^r[f^\delta(x+\rho)-f(x+\rho)]d\rho, \quad U_r(x)=\frac{2}{r^2}\int_0^r\int_0^\rho[u(x+s)-u(x)]dsd\rho.\nonumber $

方程(4.3) 两边关于误差函数 $e^{r, \delta}$做内积, 由算子 $A$的非负性可知

$ (e^{r, \delta}, e^{r, \delta})\leq(e^{r, \delta}, F_r^\delta)+(e^{r, \delta}, U_r), $

进而有 $\|e^{r, \delta}\|\leq\|F_r^\delta\|+\|U_r\|.$借助Hölder不等式可估计 $F_r^\delta$和 $U_r$的上界, 其中

$ \begin{eqnarray} \|F_r^\delta\|^2 &=& \frac{4}{r^4}\int_0^1\{\int_0^r[f^\delta(x+\rho)-f(x+\rho)]d\rho\}^2dx \nonumber \\ &\leq& \frac{4}{r^4}\int_0^1r\int_0^r[f^\delta(x+\rho)-f(x+\rho)]^2d\rho dx \nonumber \\ &=& \frac{4}{r^3}\int_0^r\int_\rho^{1+\rho}[f^\delta(x)-f(x)]^2dxd\rho \leq \frac{4\delta^2}{r^2}, \nonumber \end{eqnarray} $

因而有 $\|F_r^\delta\|\leq\frac{2\delta}{r}$.关于第二项 $U_r$, 有

$ \begin{eqnarray} \|U_r\|^2 &=& \frac{4}{r^4}\int_0^1\{\int_0^r\int_0^\rho[u(x+s)-u(x)]dsd\rho\}^2dx \nonumber \\ &\leq& \frac{4}{r^4}\int_0^1r\int_0^r\rho\int_0^\rho[u(x+s)-u(x)]^2dsd\rho dx \nonumber \\ &=& \frac{4}{r^3}\int_0^r\rho\int_0^\rho\int_0^1[u(x+s)-u(x)]^2dxdsd\rho. \nonumber \end{eqnarray} $

记 $u_s(x)=u(x+s)-u(x), x\in[0,1], s\in(0, R]$, 则有 $u_s\in H^1[0,1]$, 存在常数 $C_1$以及线性映射 $E: H^1[0,1]\rightarrow H^1(R)$, 使得 ${\left\| {{u_s}} \right\|_{{H^1}[0,1]}} \le {\left\| {E{u_s}} \right\|_{{H^1}(R)}} \le {C_1}{\left\| {{u_s}} \right\|_{{H^1}[0,1]}}$成立^[17].记 $Eu_s$的Fourier变换为 $\widehat{Eu}_s$, 则

$ \begin{eqnarray} \|u_s\|^2 &\leq& \|Eu_s\|_{L^2(R)}^2=\|\widehat{Eu}_s\|_{L^2(R)}^2 = \int_R|\widehat{Eu}(\xi)\cdot e^{is\xi}-\widehat{Eu}(\xi)|^2d\xi \nonumber \\ &=& \int_R|\widehat{Eu}(\xi)|^2\cdot|1-e^{is\xi}|^2d\xi = \int_R|\widehat{Eu}(\xi)|^2\cdot[1-\cos(s\xi)]d\xi \nonumber \\ &\leq& s^2\int_R|\widehat{Eu}(\xi)|^2\xi^2d\xi \leq s^2\|Eu\|_{H^1(R)}\leq C_1s^2\|u\|_{H^1[0,1]}\leq C_1Ms^2. \nonumber \end{eqnarray} $

因此 $ \|U_r\|^2\leq\frac{4MC_1}{r^3}\int_0^r\rho\int_0^\rho s^2dsd\rho=\frac{MC_1}{3}r^2.$进而有 $ \|u^{r, \delta}-u\|\leq\frac{2\delta}{r}+Cr, $进一步当取 $r(\delta)=\sqrt{2\delta/C}$时, 有 $\|u^{r, \delta}-u\|=O(\sqrt\delta)$.定理证毕.

在本节中, 我们将通过两个算例说明局部正则化方法求解一阶微分问题的数值有效性.首先考虑第二节中提到的算例.

算例1  考虑函数 $f(x)=|x-\frac{1}{2}|-\frac{1}{2}, x\in[0,1]$, 对其添加不同水平的高斯白噪声可得模拟的测量数据 $f^\delta$, 分别利用不同的正则化方法求其一阶导数.

将局部正则化方法求得的正则化解以及第二节中利用Tikhonov正则化和Lavrentiev正则化方法求得的正则化解分别记为 $u_{{\rm Loc}}$, $u_{{\rm Tikh}}$和 $u_{{\rm Lavr}}$.当噪声水平为 $1\%$时, 我们在图 2中给出了 $u_{{\rm Tikh}}$和 $u_{{\rm Loc}}$的图像, 其中局部正则化和Lavrentiev正则化方法的参数均先验取为 $\sqrt\delta/5$, Tikhonov正则化方法的参数取为该值的四分之一.从图 1和图 2中可以看出, 与两类经典的正则化方法相比, 局部正则化方法能在有效抑制噪声的同时, 很好地反应精确导数有间断的部分.为进一步说明局部正则化方法的有效性, 我们在表 1中比较了不同噪声水平下, 不同正则化解的相对误差.以 $u_{{\rm Loc}}$为例, 其相对误差定义为 $\|u_{{\rm Loc}}-u\|/\|u\|$, 该误差的大小能同时反映近似导数的精度和稳定性.从表 1中可以看出, 在 $f$的不同噪声水平下, 近似导数 $u_{{\rm Loc}}$的相对误差是最小的.随着噪声水平 $\delta$的提高, $u_{{\rm Lavr}}$的相对误差的增加速度是最快的, 这说明该正则化方法抑制噪声的能力较弱, 对噪声过于敏感.虽然Tikhonov在抑制噪声方面的能力是最强的, 但当噪声水平较小时, $u_{{\rm Tikh}}$的相对误差最大, 这与该正则化方法的过度光滑化不无关系.事实上, 正则化解 $u_{{\rm Loc}}$在抑制噪声方面的能力与 $u_{{\rm Tikh}}$相仿.综上所述, 与 $u_{{\rm Lavr}}$和 $u_{{\rm Tikh}}$相比, $u_{{\rm Loc}}$在抑制噪声和逼近精确导数两方面具有明显的优越性.

算例2  考虑

$ \begin{equation} f(x)=\left\{\begin{aligned} -\frac{1}{2}x^2+\frac{x}{2}, \quad\quad x\in[0, \frac{1}{2}], \\ \frac{1}{2}x^2-\frac{x}{2}+\frac{1}{4}, \quad x\in(\frac{1}{2}, 1]. \end{aligned} \right.\nonumber \end{equation} $

对其添加不同水平的高斯白噪声可得模拟的测量数据 $f^\delta$, 分别利用不同的正则化方法求其一阶导数.

当噪声水平为 $1\%$时, 我们在图 3中给出了近似导数 $u_{{\rm Tikh}}$, $u_{{\rm Lavr}}$, $u_{{\rm Loc}}$以及精确导数 $u$的图像, 其中正则化参数的选取与算例1相同.从图中可以看出, 局部正则化方法能在有效抑制噪声的同时, 较好地反应精确导数的局部变化情况.为进一步说明局部正则化方法的有效性, 我们在表 2中比较了不同噪声水平下, 不同正则化解的相对误差.类似于算例1的分析可知 $u_{{\rm Loc}}$在抑制噪声和逼近精确导数两方面具有优越性.

在本文中, 我们给出了求解一阶数值微分问题的局部正则化方法.在可积函数空间中讨论了正则化参数的先验选取策略及相应近似导数的误差估计.数值实验表明相对于经典的Tikhonov正则化和Lavrentiev正则化方法, 局部正则化方法求得的近似导数能在有效抑制噪声的同时, 很好地反应精确导数有间断或急剧变化的部分.鉴于该微分算法在保持近似导数的精度和稳定性两方面的优越性, 我们将进一步考虑该算法在具体问题如图像边缘检测、医学成像反演问题等中的应用.此外, 正则化参数的后验选取策略及相应近似导数的误差分析也是我们将要考虑的问题.