Loading [MathJax]/jax/output/HTML-CSS/jax.js
  数学物理学报  2015, Vol. 35 Issue (6): 1207-1219   PDF (446 KB)    
扩展功能
加入收藏夹
复制引文信息
加入引用管理器
Email Alert
RSS
本文作者相关文章
梁红
饶世钧
差值定理在离散数据二阶导数解算中的应用
梁红1 , 饶世钧2    
1 中国人民解放军 91550部队94分队 辽宁大连 116023;
2 海军大连舰艇学院信息作战系 辽宁大连 116018
收稿日期: 2015-04-24; 修订日期: 2015-10-29
作者简介: 梁红,springnote@126.com;饶世钧,rsj0207@163.com
摘要: 为解决离散条件下应用差值定理进行微分求导过程中受测量误差影响,导致计算结果不准确或某些测量数据无法解算的问题,在等间隔采样条件下对差值定理的适用条件进行了分析和推导,得出了判断差值曲线近似拐点的充要条件;分析得出了数据含有测量误差的情况下差值定理的适用条件,采用测试函数进行了验证,并与中心平滑算法和三阶样条算法以及不同的拟合模型的解算结果进行了比较;通过实测数据解算,分析并验证了影响差值定理使用效果的影响因素;提出了应用差值定理及其推论解算二阶导数的具体方法.
关键词: 差值定理     导数     数字滤波     截断误差     拐点    
Application of the Difference Theorem to Calculating Second Order Derivatives of Discrete Data
Liang Hong1, Rao Shijun2     
1 PLA Unit 91550, Liaoning Dalian 116023;
2 Department of Communication and Operation, Dalian Naval Academy, Liaoning Dalian 116018
Received: 2015-04-24; Revised:2015-10-29
Abstract: When the Difference Theorem is applied to calculating the derivative of discrete data, the measure-error can lead to inaccurate result or some of the data having no result. To resolve these problems, the necessary and sufficient conditions for identifying approximative inflexions of the difference-curve are obtained by analyzing and deducing applying conditions of the Difference Theorem on condition that the data are at equal intervals. In the condition of measure-error existing, the applying conditions of the Difference Theorem are obtained by analyzing, and test and verified by using test-function. The results are compared with those from center-smoothing algorithm, cubic spline algorithm and different fitting models. The factors which can affect application effects of the Difference Theorem are analyzed and validated by using measurement data. The detailed method of calculating second order derivative with the Difference Theorem and it's deduction is put forward.
Key words: The difference theorem     Derivative     Digital filtering     Truncation error     Inflexion    
1 引言

二阶导数的解算有着十分广泛的应用, 在光谱分析、大地测量、 航空航天、机械制造、图像处理等众多领域都有很高的应用价值. 由于实际应用中,测量数据往往是以离散点的形式给出的, 要获得准确的计算结果具有相当的难度,因此, 离散数据的二阶导数解算方法在工程界倍受关注. 目前, 比较常用的离散数据二阶导数解算方法是采用近似函数对测量数据进行逼近, 然后再进行微分求导, 比如数值微分公式的外推算法[1]、中心差商方法[2, 3]、余弦变换法[4] 等等. 但是,当测量数据变化非常剧烈或变化规律比较复杂时, 近似函数往往难以对其进行比较精确的拟合, 这就很容易导致二阶导数的解算结果引入较大的截断误差. 由于近似函数不准确,多年来, 截断误差问题一直是困扰数据处理工作的一个棘手的难题. 2004年, 差值定理及其推论[5]的发现,为解决截断误差问题提供了良好的途径. 采用差值定理及其推论作为二阶导数解算的有效点判据[5], 具有很大的优越性,不仅能有效降低数据剧烈变化段微分求导的截断误差, 而且受近似函数的形式及拟合区间的大小等因素的影响很小,解算精度很高. 但是,由于测量误差的影响,可能会导致一些测量点被错误地判定为有效点, 给计算结果带来误差,有时甚至不能得到解算结果. 要解决这个问题, 至少有两种途径: 一是使测量数据满足一定的精度要求, 二是从算法上采取相应的措施. 本文立足于此,在等间隔采样的前提下, 从算法着手,研究了应用差值定理及其推论时, 在何种情况下测量误差不会干扰二阶导数的解算结果, 并对应用差值定理的一些相关问题做了初步的研究.

2 差值定理及其推论

定理 2.1 设函数F(t)和f(t)在某定义域 Ω内的任一点均存在n∈N+阶导数,且 G(t)=F(n−1)(t)−f(n−1)(t),则F(n)(t0)=f(n)(t0)⇔G′(t0)=0(t0∈Ω).

证

G(t)=F(n−1)(t)−f(n−1)(t)F(t)和f(t)在定义域Ω内的任一点均存在n阶导数}

⇒G′(t)=F(n)(t)−f(n)(t)t0∈Ω}

⇒G′(t0)=F(n)(t0)−f(n)(t0),

故

F(n)(t0)=f(n)(t0)⇔F(n)(t0)−f(n)(t0)=0⇔G′(t0)=0.

证毕.

推论 2.1 设函数 G(t)=F(t)−f(t)的定义域为Ω,且函数F(t) 和f(t)的一阶导数、二阶导数均存在,则F′(t1)=f′(t1)⇔G′(t1)=0(t1∈Ω) 且F″(t2)=f″(t2)⇔G″(t2)=0(t2∈Ω).

定理2.1即差值定理,推论2.1是差值定理的推论.

差值定理及其推论有明确的物理含义. 差值定理的物理意义是``两个函数的n−1 阶导数的差值曲线的驻点是它们的n阶导数相等的点", 其推论的物理意义是``两个函数的差值曲线的驻点是它们的一阶导数相等的点, 差值曲线的拐点是它们的二阶导数相等的点".

应用差值定理及其推论求解函数f(t)的二阶导数时, 用F(t)来表示拟合的近似函数,则根据差值定理的推论, 只需找到函数F(t)与f(t)的差值曲线的拐点, 即可把F(t)在该点的二阶导数近似作为f(t)的二阶导数. 因此, 离散测量数据的二阶导数解算问题就转化为如何克服测量误差的影响, 比较准确地判断差值曲线的拐点.

3 差值定理适用条件的分析与推导

3.1 判断拐点的充要条件的理论分析

对于离散数据, 受采样率的限制,难以准确找出差值曲线的拐点, 故把拐点的一个邻近采样点近似作为拐点,称为近似拐点. 由于曲线的拐点是其斜率的极值点,因此,在判断差值曲线的拐点时, 不需知道差值曲线斜率的确切值,而只要知道相邻点的斜率孰大孰小即可.

设k0为近似拐点的斜率,k−1为近似拐点前一点的斜率,k+1为近似拐点后一点的斜率,则离散数据近似拐点的判别条件为

{k−1>k0,k+1>k0

或

{k−1<k0,k+1<k0.

可采用中心平滑[6]直线拟合法求差值曲线上某点的斜率,故上述判别条件可用下式表示

{(2N+1)N−1∑i=−N−1tiyi−(N−1∑i=−N−1ti)(N−1∑i=−N−1yi)(2N+1)N−1∑i=−N−1t2i−(N−1∑i=−N−1ti)2>(2N+1)N∑i=−Ntiyi−(N∑i=−Nti)(N∑i=−Nyi)(2N+1)N∑i=−Nt2i−(N∑i=−Nti)2,(2N+1)N+1∑i=−N+1tiyi−(N+1∑i=−N+1ti)(N+1∑i=−N+1yi)(2N+1)N+1∑i=−N+1t2i−(N+1∑i=−N+1ti)2>(2N+1)N∑i=−Ntiyi−(N∑i=−Nti)(N∑i=−Nyi)(2N+1)N∑i=−Nt2i−(N∑i=−Nti)2 (3.1)

或

{(2N+1)N−1∑i=−N−1tiyi−(N−1∑i=−N−1ti)(N−1∑i=−N−1yi)(2N+1)N−1∑i=−N−1t2i−(N−1∑i=−N−1ti)2<(2N+1)N∑i=−Ntiyi−(N∑i=−Nti)(N∑i=−Nyi)(2N+1)N∑i=−Nt2i−(N∑i=−Nti)2,(2N+1)N+1∑i=−N+1tiyi−(N+1∑i=−N+1ti)(N+1∑i=−N+1yi)(2N+1)N+1∑i=−N+1t2i−(N+1∑i=−N+1ti)2<(2N+1)N∑i=−Ntiyi−(N∑i=−Nti)(N∑i=−Nyi)(2N+1)N∑i=−Nt2i−(N∑i=−Nti)2, (3.2)

式中 N表示求斜率时中心平滑直线拟合的滤波半径,ti 表示差值曲线上下标为i的点的横坐标,yi表示差值曲线上下标为i的点的纵坐标.

设采样间隔为h,则式(3.1)可以化为

{(2N+1)hN−1∑i=−N−1iyi−h(N−1∑i=−N−1i)(N−1∑i=−N−1yi)(2N+1)h2N−1∑i=−N−1i2−h2(N−1∑i=−N−1i)2>(2N+1)hN∑i=−Niyi−h(N∑i=−Ni)(N∑i=−Nyi)(2N+1)h2N∑i=−Ni2−h2(N∑i=−Ni)2,(2N+1)hN+1∑i=−N+1iyi−h(N+1∑i=−N+1i)(N+1∑i=−N+1yi)(2N+1)h2N+1∑i=−N+1i2−h2(N+1∑i=−N+1i)2>(2N+1)hN∑i=−Niyi−h(N∑i=−Ni)(N∑i=−Nyi)(2N+1)h2N∑i=−Ni2−h2(N∑i=−Ni)2,

即

{N∑i=−Nyi−Ny−N−1−(N+1)yN>0,N∑i=−Nyi−NyN+1−(N+1)y−N<0.

同理,式(3.2)可化为

{N∑i=−Nyi−Ny−N−1−(N+1)yN<0,N∑i=−Nyi−NyN+1−(N+1)y−N>0.

故差值曲线上下标为0的点为近似拐点的必要条件是

{N∑i=−Nyi−Ny−N−1−(N+1)yN>0,N∑i=−Nyi−NyN+1−(N+1)y−N<0 (3.3)

或

{N∑i=−Nyi−Ny−N−1−(N+1)yN<0,N∑i=−Nyi−NyN+1−(N+1)y−N>0. (3.4)

将上面的推导过程反推可知,必要条件也是充分条件.

3.2 实测数据拐点判断的充要条件

用实测数据计算差值曲线时, 往往会引入测量误差. 用ˉyi 表示差值曲线上下标为i的点含有误差时的纵坐标,用si表示下标为i的点的测量误差,则

yi=ˉyi−si (3.5)

将式(3.5)代入式(3.3),得

{N∑i=−N(ˉyi−si)−N(ˉy−N−1−s−N−1)−(N+1)(ˉyN−sN)>0,N∑i=−N(ˉyi−si)−N(ˉyN+1−sN+1)−(N+1)(ˉy−N−s−N)<0,

即

{N∑i=−Nˉyi−Nˉy−N−1−(N+1)ˉyN>N−1∑i=−Nsi−Ns−N−1−NsN,N∑i=−Nˉyi−NˉyN+1−(N+1)ˉy−N<N∑i=−N+1si−NsN+1−Ns−N. (3.6)

同理,将式(3.5)代入式(3.4),得

{N∑i=−Nˉyi−Nˉy−N−1−(N+1)ˉyN<N−1∑i=−Nsi−Ns−N−1−NsN,N∑i=−Nˉyi−NˉyN+1−(N+1)ˉy−N>N∑i=−N+1si−NsN+1−Ns−N. (3.7)

将上面的推导过程反推可得式(3.3)和式(3.4),因此, 对于实测数据, 判断差值曲线上近似拐点的充要条件是式(3.6)和式(3.7).

3.3 差值定理的适用条件

设测量误差限为S>0,即|si|≤S,则由式(3.6)得

{max(N−1∑i=−Nsi−Ns−N−1−NsN)=4NS,min(N∑i=−N+1si−NsN+1−Ns−N)=−4NS,

且

{min(N−1∑i=−Nsi−Ns−N−1−NsN)=−4NS,max(N∑i=−N+1si−NsN+1−Ns−N)=4NS.

故当

{N∑i=−Nˉyi−Nˉy−N−1−(N+1)ˉyN>4NS,N∑i=−Nˉyi−NˉyN+1−(N+1)ˉy−N<−4NS (3.8)

或

{N∑i=−Nˉyi−Nˉy−N−1−(N+1)ˉyN<−4NS,N∑i=−Nˉyi−NˉyN+1−(N+1)ˉy−N>4NS (3.9)

成立时,必使式(3.6)或式(3.7)成立.

设 f(ti)的测量值为ˉf(ti),则

ˉyi=ˉf(ti)−F(ti). (3.10)

将式(3.10)代入式(3.8)和式(3.9),得

{N∑i=−N(ˉf(ti)−F(ti))−N(ˉf(t−N−1)−F(t−N−1))−(N+1)(ˉf(tN)−F(tN))>4NS,N∑i=−N(ˉf(ti)−F(ti))−N(ˉf(tN+1)−F(tN+1))−(N+1)(ˉf(t−N)−F(t−N))<−4NS (3.11)

或

{N∑i=−N(ˉf(ti)−F(ti))−N(ˉf(t−N−1)−F(t−N−1))−(N+1)(ˉf(tN)−F(tN))<−4NS,N∑i=−N(ˉf(ti)−F(ti))−N(ˉf(tN+1)−F(tN+1))−(N+1)(ˉf(t−N)−F(t−N))>4NS. (3.12)

式(3.11)和式(3.12)是判断差值曲线上下标为0的点为近似拐点的充分条件, 即当测量数据满足式(3.11)或式(3.12)时,测量误差不会干扰近似拐点的判断, 且此时下标为0 的点是差值曲线的近似拐点.

4 仿真数据验证

4.1 基本原理的验证

应用差值定理及其推论作为有效点判据, 对测试函数f(t)=(t/10)−0.01t+10.1/1000+100000 (采样间隔为1) 进行二阶导数解算. 对f(t)进行最小二乘3阶拟合, 拟合区间为30∼50,其差值曲线绘入图 1, 由三阶多项式微分得到的二阶导数及其拟合误差列入表 1.

图 1 函数f(t)=(t/10)−0.01t+10.1/1000+100000 的三阶拟合差值曲线

表1 测试函数f(t)=(t/10)−0.01t+10.1/1000+100000的二阶导数微分解算结果

对差值曲线采用中心平滑直线拟合,滤波半径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 算法测试结果

由表 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.

表3 应用差值定理的解算结果

以11个点为拟合区间,采用3阶多项式拟合,以中心平滑方法计算f(t)的2阶导数; 将拟合区间向后滑动,用同样的方法获得每个区间的中点的2阶导数解算结果,列入表 4.

表4 应用中心平滑算法的解算结果

以25个点为拟合区间,5个点为1个样条的长度,拟合区间端点处的1阶导数为已知条件, 采用3阶样条算法计算f(t)节点处的2阶导数; 将拟合区间向后滑动, 用同样的方法获得每个区间的节点处的2阶导数解算结果,列入表 5.

表5 应用3阶样条算法的解算结果

由表 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 采用3阶多项式拟合模型的解算结果

表7 采用拟合模型ae0.105t+100000.0的解算结果

表8 采用2阶多项式拟合模型的解算结果

由表 6、表 7、表 8可知,当采用差值定理进行2阶导数解算时,计算结果受拟合模型的影响较小, 解算精度仍然较好.

上述分析表明,采用差值定理进行2阶导数解算是一种性能比较稳定、准确的方法, 式(3.11)和式(3.12)算法正确,可以作为判断近似拐点的充分条件使用.

5 实测数据解算与分析

用差值定理及其推论作为有效点判据,对一组实测数据进行二阶导数解算, 改变多项式拟合阶次、对差值曲线进行中心平滑直线拟合的滤波半径和测量误差限S的取值, 将解算的近似拐点数列入表 9.

表9 某实测数据解算结果统计表

由表 9可看出,在同一拟合区间中,当拟合阶次相同时, 调节差值曲线直线拟合的滤波半径后进行解算,可以弥补部分未解算出来的二阶导数; 在同一拟合区间中,当差值曲线直线拟合的滤波半径相同时,调整拟合阶次, 即改变拟合模型,也可弥补部分未解算出来的二阶导数; 在同一拟合区间中, 当拟合阶次和差值曲线直线拟合的滤波半径均相同时,随着S的增大, 解算结果会越来越少,因此,S的值应当由测量误差限决定, 在计算时不宜随意改变,也不宜将S不同时的解算结果合并在一起.

将上述S=0.2的各解算结果合并后,绘入图 2.

图 2 二阶导数的解算结果

由图 2可看出,将各解算结果合并后所得的结果没有出现明显的跳点, 计算结果良好. 这说明: 对于实测数据,调节差值曲线直线拟合的滤波半径或改变测量数据的拟合模型, 均有可能弥补未解算出来的二阶导数; 而且,对于等间隔采样, 采样间隔h对近似拐点的判断没有影响,且近似拐点的判断与测量值、 拟合值、差值曲线中心平滑直线拟合的滤波半径N及测量误差有关. 因此,对于含有测量误差的实测数据,调节滤波半径N或改变拟合模型, 有可能使判别近似拐点的充分条件得到满足,从而使一些无法解算的测量数 据获得二阶导数的解算结果.

6 二阶导数解算方法

综上所述,可得应用差值定理及其推论解算测量数据二阶导数的具体步骤,即 先用近似曲线F(t)在最小二乘原则下对测量数据进行区间拟合,然后按式 (3.11)或式(3.12)判断是否满足差值定理的适用条件,即是否满 足差值曲线上近似拐点判别的充分条件; 若满足,则表明下标为0的点是近似拐点, 此时可将拟合曲线在该点的二阶导数作为测量数据二阶导数的解算结果; 然后将拟合 区间向后移动,用上述方法找到每一区间差值曲线的近似拐点所对应的二阶导数作为该 点的二阶导数解算结果; 若有一些点没有获得解算结果,则在一定范围内改变差值曲线 的直线拟合滤波半径N或改变测量数据的拟合函数,即改变拟合值,用前述方法重新计算, 尽可能获得解算结果; 最后,将各次二阶导数的解算结果合并,获得最终的解算结果.

滤波半径N的选择越小越好,一则可缩短计算时间,二可尽量避免差值曲线在滤波半 径N的范围内有多个近似拐点时给计算带来的不利影响. 对于近似曲线F(t), 只要求其在差值曲线的近似拐点处的二阶导数能够逼近测量数据的二阶导数即可, 对其它部位则没有要求,故一般情况下,F(t)可采用各种适宜的拟合模型, 如适宜阶次的多项式.

当差值曲线在滤波半径N的范围内可能有多个近似拐点出现时,应当采用强化条件

{k−N>k−N+1,⋮k−2>k−1;k−1>k0,k+1>k0; k+2>k+1,⋮kN>kN−1 (6.1)

或

{k−N<k−N+1,⋮k−2<k−1;k−1<k0,k+1<k0;k+2<k+1,⋮kN<kN−1 (6.2)

对差值曲线的近似拐点进行判断, 以消除各近似拐点之间的相互干扰. 在此不详细讨论.

式(6.1)和式(6.2)中: k−N表示近似拐点前第N点的斜率; kN表示近似拐点后第N点的斜率; N表示求斜率时中心平滑直线拟合的滤波半径.

7 结论

当测量数据等间隔采样时, 若采用中心平滑直线拟合法求解测量数据与近似拟合函数 的差值曲线的斜率, 则

(1) 式(3.6)和式(3.7)是判断差值曲线上近似拐点的充要条件.

(2) 当满足式(3.11)或式(3.12)的条件时, 测量误差不干扰近似拐点的判断, 且式(3.11)和式(3.12)是判断近似拐点的充分条件.

(3) 近似拐点的判断条件与采样间隔的大小无关, 与测量值、近似函数拟合值、 差值曲线直线拟合的滤波半径及测量误差有关.

参考文献
[1] 王燕. 二阶导数的五点数值微分公式及外推算法. 天津理工大学学报, 2009, 25(4):37-39
[2] 姜洁. 二阶导数的中心差商. 锦州师范学院学报, 2000, 21(1):59-60
[3] 王兴华, 王何宇. 数值差商公式研究. 中国科学 A辑:数学, 2005, 35(6):712-720
[4] 张卫星. 计算重力异常垂向二阶导数的DCT法-模型实验. 甘肃地质, 2009, 18(3):85-88
[5] 梁红. 利用差值定理降低飞行器速度和加速度拟合的截断误差. 飞行器测控学报, 2005, 24(3):51-54
[6] 孙中豪, 杜娟, 王子龙. 基于白噪声正交多项式滤波的GPS测速方法分析. 测绘地理信息, 2013, 38(5):21-24
差值定理在离散数据二阶导数解算中的应用
梁红 , 饶世钧