二阶导数的解算有着十分广泛的应用, 在光谱分析、大地测量、 航空航天、机械制造、图像处理等众多领域都有很高的应用价值. 由于实际应用中,测量数据往往是以离散点的形式给出的, 要获得准确的计算结果具有相当的难度,因此, 离散数据的二阶导数解算方法在工程界倍受关注. 目前, 比较常用的离散数据二阶导数解算方法是采用近似函数对测量数据进行逼近, 然后再进行微分求导, 比如数值微分公式的外推算法[1]、中心差商方法[2, 3]、余弦变换法[4] 等等. 但是,当测量数据变化非常剧烈或变化规律比较复杂时, 近似函数往往难以对其进行比较精确的拟合, 这就很容易导致二阶导数的解算结果引入较大的截断误差. 由于近似函数不准确,多年来, 截断误差问题一直是困扰数据处理工作的一个棘手的难题. 2004年, 差值定理及其推论[5]的发现,为解决截断误差问题提供了良好的途径. 采用差值定理及其推论作为二阶导数解算的有效点判据[5], 具有很大的优越性,不仅能有效降低数据剧烈变化段微分求导的截断误差, 而且受近似函数的形式及拟合区间的大小等因素的影响很小,解算精度很高. 但是,由于测量误差的影响,可能会导致一些测量点被错误地判定为有效点, 给计算结果带来误差,有时甚至不能得到解算结果. 要解决这个问题, 至少有两种途径: 一是使测量数据满足一定的精度要求, 二是从算法上采取相应的措施. 本文立足于此,在等间隔采样的前提下, 从算法着手,研究了应用差值定理及其推论时, 在何种情况下测量误差不会干扰二阶导数的解算结果, 并对应用差值定理的一些相关问题做了初步的研究.
定理 2.1 设函数$F(t)$和$f(t)$在某定义域 $\Omega$内的任一点均存在$n \in N^{+}$阶导数,且 $G(t)= F^{(n-1)}(t)- f^{ (n-1)}(t)$,则$F^{(n)}(t_0 )= f^{ (n)}(t_0 ) \Leftrightarrow {G}'(t_0 )=0 (t_0 \in \Omega ).$
证
$$\left. {\begin{array}{*{20}{c}} {G(t) = {F^{(n - 1)}}(t) - {f^{(n - 1)}}(t)}&{{\rm{ }}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt}{\kern 1pt} }\\ {F(t){\rm{和}}f(t){\rm{在定义域 }}\Omega {\rm{ }}内的任一点均存在n阶导数}{} \end{array}} \right\} $$
$$ \left. {{\begin{array}{*{20}c} {\Rightarrow {G}'(t)=F^{(n)}(t)-f^{(n)}(t)} \hfill \\ {t_0 \in \Omega } \hfill \\ \end{array} }} \right\} $$
$$ \Rightarrow {G}'(t_0 )=F^{(n)}(t_0 )-f^{(n)}(t_0 ), $$
故
$$ F^{(n)}(t_0 )=f^{(n)}(t_0 )\Leftrightarrow F^{(n)}(t_0 )-f^{(n)}(t_0 )=0\Leftrightarrow {G}'(t_0 )=0. $$
证毕.
推论 2.1 设函数 $G(t)=F(t)-f(t)$的定义域为$\Omega $,且函数$F(t)$ 和$f(t)$的一阶导数、二阶导数均存在,则${F}'(t_1 )={f}'(t_1 )\Leftrightarrow {G}'(t_1 )=0(t_1 \in \Omega )$ 且${F}''(t_2 )={f}''(t_2 )\Leftrightarrow {G}''(t_2 )=0(t_2 \in \Omega )$.
定理2.1即差值定理,推论2.1是差值定理的推论.
差值定理及其推论有明确的物理含义. 差值定理的物理意义是``两个函数的$n-1$ 阶导数的差值曲线的驻点是它们的$n$阶导数相等的点", 其推论的物理意义是``两个函数的差值曲线的驻点是它们的一阶导数相等的点, 差值曲线的拐点是它们的二阶导数相等的点".
应用差值定理及其推论求解函数$f(t)$的二阶导数时, 用$F(t)$来表示拟合的近似函数,则根据差值定理的推论, 只需找到函数$F(t)$与$f(t)$的差值曲线的拐点, 即可把$F(t)$在该点的二阶导数近似作为$f(t)$的二阶导数. 因此, 离散测量数据的二阶导数解算问题就转化为如何克服测量误差的影响, 比较准确地判断差值曲线的拐点.
3.1 判断拐点的充要条件的理论分析
对于离散数据, 受采样率的限制,难以准确找出差值曲线的拐点, 故把拐点的一个邻近采样点近似作为拐点,称为近似拐点. 由于曲线的拐点是其斜率的极值点,因此,在判断差值曲线的拐点时, 不需知道差值曲线斜率的确切值,而只要知道相邻点的斜率孰大孰小即可.
设$k_0 $为近似拐点的斜率,$k_{-1} $为近似拐点前一点的斜率,$k_{+1} $为近似拐点后一点的斜率,则离散数据近似拐点的判别条件为
$$\left\{ {\begin{array}{l} k_{-1} >k_0 ,\\ k_{+1} >k_0 \\ \end{array}} \right. $$
或
$$ \left\{ {\begin{array}{*{20}{l}} {{k_{ - 1}} < {k_0},}\\ {{k_{ + 1}} < {k_0}.} \end{array}} \right. $$
可采用中心平滑[6]直线拟合法求差值曲线上某点的斜率,故上述判别条件可用下式表示
式中 $N$表示求斜率时中心平滑直线拟合的滤波半径,$t_i $ 表示差值曲线上下标为$i$的点的横坐标,$y_i $表示差值曲线上下标为$i$的点的纵坐标.
设采样间隔为$h$,则式(3.1)可以化为
$$ \left\{ {\begin{array}{l} \frac{(2N+1)h\sum\limits_{i=-N-1}^{N-1} {iy_i } -h(\sum\limits_{i=-N-1}^{N-1} i )(\sum\limits_{i=-N-1}^{N-1} {y_i } )}{(2N+1)h^2\sum\limits_{i=-N-1}^{N-1} {i^2} -h^2(\sum\limits_{i=-N-1}^{N-1} i )^2}>\frac{(2N+1)h\sum\limits_{i=-N}^N {iy_i } -h(\sum\limits_{i=-N}^N i )(\sum\limits_{i=-N}^N {y_i } )}{(2N+1)h^2\sum\limits_{i=-N}^N {i^2} -h^2(\sum\limits_{i=-N}^N i )^2} ,\\ \frac{(2N+1)h\sum\limits_{i=-N+1}^{N+1} {iy_i } -h(\sum\limits_{i=-N+1}^{N+1} i )(\sum\limits_{i=-N+1}^{N+1} {y_i } )}{(2N+1)h^2\sum\limits_{i=-N+1}^{N+1} {i^2} -h^2(\sum\limits_{i=-N+1}^{N+1} i )^2}>\frac{(2N+1)h\sum\limits_{i=-N}^N {iy_i } -h(\sum\limits_{i=-N}^N i )(\sum\limits_{i=-N}^N {y_i } )}{(2N+1)h^2\sum\limits_{i=-N}^N {i^2} -h^2(\sum\limits_{i=-N}^N i )^2} , \end{array}} \right. $$
即
$$ \left\{ {\begin{array}{l} \sum\limits_{i=-N}^N {y_i } -Ny_{-N-1} -(N+1)y_N >0,\\ \sum\limits_{i=-N}^N {y_i } -Ny_{N+1} -(N+1)y_{-N} <0 . \end{array}} \right. $$
同理,式(3.2)可化为
$$ \left\{ {\begin{array}{l} \sum\limits_{i=-N}^N {y_i } -Ny_{-N-1} -(N+1)y_N <0,\\[4mm] \sum\limits_{i=-N}^N {y_i } -Ny_{N+1} -(N+1)y_{-N} >0 . \end{array}} \right. $$
故差值曲线上下标为$0$的点为近似拐点的必要条件是
将上面的推导过程反推可知,必要条件也是充分条件.
3.2 实测数据拐点判断的充要条件
用实测数据计算差值曲线时, 往往会引入测量误差. 用$\bar {y}_i $ 表示差值曲线上下标为$i$的点含有误差时的纵坐标,用$s_i $表示下标为$i$的点的测量误差,则
将式(3.5)代入式(3.3),得
$$\left\{ {\begin{array}{l} \sum\limits_{i=-N}^N {(\bar {y}_i -s_i )} -N(\bar {y}_{-N-1} -s_{-N-1} )-(N+1)(\bar {y}_N -s_N )>0 ,\\[4mm] \sum\limits_{i=-N}^N {(\bar {y}_i -s_i )} -N(\bar {y}_{N+1} -s_{N+1} )-(N+1)(\bar {y}_{-N} -s_{-N} )<0 , \end{array}} \right. $$
同理,将式(3.5)代入式(3.4),得
将上面的推导过程反推可得式(3.3)和式(3.4),因此, 对于实测数据, 判断差值曲线上近似拐点的充要条件是式(3.6)和式(3.7).
3.3 差值定理的适用条件
设测量误差限为$S>0$,即$\vert s_i \vert\le S$,则由式(3.6)得
$$\left\{ {{\begin{array}{ll} \max \bigg(\sum\limits_{i=-N}^{N-1} {s_i } -Ns_{-N-1} -Ns_N\bigg )=4NS ,\\[3mm] \min \bigg(\sum\limits_{i=-N+1}^N {s_i } -Ns_{N+1} -Ns_{-N} \bigg)=-4NS, \end{array} }} \right. $$
且
$$ \left\{ {{\begin{array}{ll} \min \bigg(\sum\limits_{i=-N}^{N-1} {s_i } -Ns_{-N-1} -Ns_N\bigg )=-4NS,\\[3mm] \max \bigg(\sum\limits_{i=-N+1}^N {s_i } -Ns_{N+1} -Ns_{-N} \bigg)=4NS . \end{array} }} \right. $$
故当
成立时,必使式(3.6)或式(3.7)成立.
设 $f(t_i)$的测量值为$\bar {f}(t_i )$,则
将式(3.10)代入式(3.8)和式(3.9),得
式(3.11)和式(3.12)是判断差值曲线上下标为0的点为近似拐点的充分条件, 即当测量数据满足式(3.11)或式(3.12)时,测量误差不会干扰近似拐点的判断, 且此时下标为0 的点是差值曲线的近似拐点.
4.1 基本原理的验证
应用差值定理及其推论作为有效点判据, 对测试函数$f(t)=(t/10)^{-0.01t+10.1}/1000+100000$ (采样间隔为1) 进行二阶导数解算. 对$f(t)$进行最小二乘3阶拟合, 拟合区间为$30 \sim 50$,其差值曲线绘入图 1, 由三阶多项式微分得到的二阶导数及其拟合误差列入表 1.
对差值曲线采用中心平滑直线拟合,滤波半径$N=1$,$S=0.00001$, 以式(3.11)和式(3.12)作为近似拐点的判断条件, 求得近似拐点为$t=36$ 和$t=44$. 从图 1可看出, $t=36$和$t=44$的确是差值曲线的近似拐点. 从表 1中的数据可看出, $t=36$和$t=44$也刚好是二阶导数拟合误差最接近于0的点.
保持拟合区间为21个点,将拟合区间向前或向后移动, 用上述方法找到每一区间差值曲线的近似拐点所对应的三阶拟合函数的二阶导数作为测 试函数$f(t)=(t/10)^{-0.01t+10.1}/1000+100000$的二阶导数解算结果, 列入表 2.
由表 2可知,在采用了并非最佳拟合阶次的情况下, 二阶导数的解算误差比较小,计算结果良好.
4.2 与其它算法的比较
以11个点为拟合区间, 采用3阶多项式进行最小二乘拟合,取$N=1$,$S=0.00001$, 以式(3.11)和式(3.12)作为近似拐点的判断条件, 对测试函数$f(t)=(t/10)^{-0.01t+10.1}/1000+100000$ (采样间隔为1) 进行2 阶导数解算; 将拟合区间向后滑动, 用同样的方法获得每个区间的近似拐点处的2 阶导数解算结果,列入表 3.
以11个点为拟合区间,采用3阶多项式拟合,以中心平滑方法计算$f(t)$的2阶导数; 将拟合区间向后滑动,用同样的方法获得每个区间的中点的2阶导数解算结果,列入表 4.
以25个点为拟合区间,5个点为1个样条的长度,拟合区间端点处的1阶导数为已知条件, 采用3阶样条算法计算$f(t)$节点处的2阶导数; 将拟合区间向后滑动, 用同样的方法获得每个区间的节点处的2阶导数解算结果,列入表 5.
由表 3、表 4、表 5可知,在拟合模型相同的情况下, 采用差值定理进行2阶导数解算, 结果甚至比中心平滑算法和3阶样条算法还要好. 而且, 将表 3与表 2中的数据比较可知,以差值定理作为二阶导数解算的有效点判据, 计算结果受拟合区间的长度影响较小.
4.3 不同拟合模型的比较
以11个点为拟合区间,取$N=1$,$S=0.00001$, 分别采用3阶多项式、2阶多项式、指数函数为拟合模型进行最小二乘拟合, 以式(3.11)和式(3.12) 作为近似拐点的判断条件, 对测试函数$f(t)=(t/10)^{-0.01t+10.1}/1000+100000$ (采样间隔为1) 进行2阶导数解算,将结果分别列入表 6、表 7、表 8.
由表 6、表 7、表 8可知,当采用差值定理进行2阶导数解算时,计算结果受拟合模型的影响较小, 解算精度仍然较好.
上述分析表明,采用差值定理进行2阶导数解算是一种性能比较稳定、准确的方法, 式(3.11)和式(3.12)算法正确,可以作为判断近似拐点的充分条件使用.
5 实测数据解算与分析
用差值定理及其推论作为有效点判据,对一组实测数据进行二阶导数解算, 改变多项式拟合阶次、对差值曲线进行中心平滑直线拟合的滤波半径和测量误差限$S$的取值, 将解算的近似拐点数列入表 9.
由表 9可看出,在同一拟合区间中,当拟合阶次相同时, 调节差值曲线直线拟合的滤波半径后进行解算,可以弥补部分未解算出来的二阶导数; 在同一拟合区间中,当差值曲线直线拟合的滤波半径相同时,调整拟合阶次, 即改变拟合模型,也可弥补部分未解算出来的二阶导数; 在同一拟合区间中, 当拟合阶次和差值曲线直线拟合的滤波半径均相同时,随着$S$的增大, 解算结果会越来越少,因此,$S$的值应当由测量误差限决定, 在计算时不宜随意改变,也不宜将$S$不同时的解算结果合并在一起.
将上述$S=0.2$的各解算结果合并后,绘入图 2.
由图 2可看出,将各解算结果合并后所得的结果没有出现明显的跳点, 计算结果良好. 这说明: 对于实测数据,调节差值曲线直线拟合的滤波半径或改变测量数据的拟合模型, 均有可能弥补未解算出来的二阶导数; 而且,对于等间隔采样, 采样间隔$h$对近似拐点的判断没有影响,且近似拐点的判断与测量值、 拟合值、差值曲线中心平滑直线拟合的滤波半径$N$及测量误差有关. 因此,对于含有测量误差的实测数据,调节滤波半径$N$或改变拟合模型, 有可能使判别近似拐点的充分条件得到满足,从而使一些无法解算的测量数 据获得二阶导数的解算结果.
综上所述,可得应用差值定理及其推论解算测量数据二阶导数的具体步骤,即 先用近似曲线$F(t)$在最小二乘原则下对测量数据进行区间拟合,然后按式 (3.11)或式(3.12)判断是否满足差值定理的适用条件,即是否满 足差值曲线上近似拐点判别的充分条件; 若满足,则表明下标为0的点是近似拐点, 此时可将拟合曲线在该点的二阶导数作为测量数据二阶导数的解算结果; 然后将拟合 区间向后移动,用上述方法找到每一区间差值曲线的近似拐点所对应的二阶导数作为该 点的二阶导数解算结果; 若有一些点没有获得解算结果,则在一定范围内改变差值曲线 的直线拟合滤波半径$N$或改变测量数据的拟合函数,即改变拟合值,用前述方法重新计算, 尽可能获得解算结果; 最后,将各次二阶导数的解算结果合并,获得最终的解算结果.
滤波半径$N$的选择越小越好,一则可缩短计算时间,二可尽量避免差值曲线在滤波半 径$N$的范围内有多个近似拐点时给计算带来的不利影响. 对于近似曲线$F(t)$, 只要求其在差值曲线的近似拐点处的二阶导数能够逼近测量数据的二阶导数即可, 对其它部位则没有要求,故一般情况下,$F(t)$可采用各种适宜的拟合模型, 如适宜阶次的多项式.
当差值曲线在滤波半径$N$的范围内可能有多个近似拐点出现时,应当采用强化条件
对差值曲线的近似拐点进行判断, 以消除各近似拐点之间的相互干扰. 在此不详细讨论.
式(6.1)和式(6.2)中: $k_{-N} $表示近似拐点前第$N$点的斜率; $k_N $表示近似拐点后第$N$点的斜率; $N$表示求斜率时中心平滑直线拟合的滤波半径.
当测量数据等间隔采样时, 若采用中心平滑直线拟合法求解测量数据与近似拟合函数 的差值曲线的斜率, 则
(1) 式(3.6)和式(3.7)是判断差值曲线上近似拐点的充要条件.
(2) 当满足式(3.11)或式(3.12)的条件时, 测量误差不干扰近似拐点的判断, 且式(3.11)和式(3.12)是判断近似拐点的充分条件.
(3) 近似拐点的判断条件与采样间隔的大小无关, 与测量值、近似函数拟合值、 差值曲线直线拟合的滤波半径及测量误差有关.
2015, Vol. 35 Issue (6): 1207-1219
PDF (446 KB)
2015, Vol. 35 Issue (6): 1207-1219
PDF (446 KB)