【人工智能基础】三维离散曲线曲率挠率的微中心差分算法.pdf_大学文库

第14卷第1期智能系统学报 Vol.14 No.I 2019年1月 CAAI Transactions on Intelligent Systems Jan.2019 D0:10.11992/tis.201802008 网络出版地址：http:/kns.cnki.net/kcms/detail/23.1538.TP.20180927.1124.006html 三维离散曲线曲率挠率的微中心差分算法慕生鹏，李红军，李世林 (北京林业大学理学院，北京100083) 摘要：曲率和挠率是描述三维空间离散曲线的弯曲和扭曲程度的两个微分量。为了准确计算这两个微分量，从连续曲线的导数定义出发，提出微中心差分算法进行三维空间离散曲线的曲率和挠率计算。该算法基于差商平滑策略实现对单侧差分算法的一个有效扩展。与单侧差分算法相比，微中心差分算法不增加算法执行时间，但在计算精度方面有显著提升。实验分析是通过6条曲线的均匀采样获取离散曲线数据，与5种常用的曲率和挠率计算算法相比较，对这6种算法从采样密度对算法精度的影响、计算效率和抗噪声性能这3个方面进行了对比分析。实验结果表明，微中心差分算法总体效果最好。关键词：曲率；挠率；算法比较：离散曲线；微中心差分法；离散几何法；三维空间；差商；均匀采样中图分类号：TP311文献标志码：A文章编号：1673-4785(2019)01-0194-13 中文引用格式：慕生鹏，李红军，李世林.三维离散曲线曲率挠率的微中心差分算法J.智能系统学报，2019,14(1)： 194-206. 英文引用格式：MU Shengpeng,LI Hongjun,LI Shilin.An algorithm for estimating curvature and torsion of discrete curve in three-dimensional space based on microcentral differencelJ].CAAI transactions on intelligent systems,2019,14(1):194-206. An algorithm for estimating curvature and torsion of discrete curve in three-dimensional space based on microcentral difference MU Shengpeng,LI Hongjun,LI Shilin (College of Science,Beijing Forestry University,Beijing 100083,China) Abstract:The curvature and torsion of a 3D discrete curve reflect the degrees of its bending and distortion.To calculate these quantities accurately,following the definition of the derivative of the continuous curve,a microcentral difference algorithm,which is an extension to the one-side difference algorithm,is proposed based on the smoothing of the differ- ence quotient.Compared with the one-side difference algorithm,the microcentral difference algorithm fails to prolong the running time but it remarkably improves the calculation accuracy.Several experiments are conducted by uniform sampling from six continuous curves,which are then compared with the five traditional algorithms of curvature and tor- sion.The experimental results are analyzed from three aspects:the influence of the sampling density on the accuracy of the algorithm,the efficiency of calculation,and the anti-noise performance.The experimental results show the good per- formance of the proposed microcentral difference algorithm. Keywords:curvature;torsion;algorithm comparison;discrete curve;microcentral difference algorithm;discrete geo- metry method;three-dimensional space;difference quotient;uniformly sampling 在人工智能算法设计与自动控制的相关研究视觉研究，以及植物生长模拟四、结构工程分析图中，无人机路径规划四，弹道分析P,公路线性设计等问题的解决都可以基于离散曲线的曲率和挠率路径约束下的车辆行为研究，几何处理和机器的分析。曲率和挠率是空间曲线在固有运动下的收稿日期：2018-02-05.网络出版日期：2018-09-29. 不变量，直接描述了曲线在一点邻近的形状。基金项目：国家自然科学基金项目(61372190)：中央高校基本科研业务费专项资金项目(2015 ZCQ-LY-01). 其中，曲率揭示曲线在所在平面的弯曲程度，挠通信作者：李红军.E-mail:lihongjun69@bjfu.edu.cn. 率则刻画曲线离开既定平面的扭曲程度。也就是

DOI: 10.11992/tis.201802008 网络出版地址: http://kns.cnki.net/kcms/detail/23.1538.TP.20180927.1124.006.html 三维离散曲线曲率挠率的微中心差分算法慕生鹏，李红军，李世林（北京林业大学理学院，北京 100083）摘要：曲率和挠率是描述三维空间离散曲线的弯曲和扭曲程度的两个微分量。为了准确计算这两个微分量，从连续曲线的导数定义出发，提出微中心差分算法进行三维空间离散曲线的曲率和挠率计算。该算法基于差商平滑策略实现对单侧差分算法的一个有效扩展。与单侧差分算法相比，微中心差分算法不增加算法执行时间，但在计算精度方面有显著提升。实验分析是通过 6 条曲线的均匀采样获取离散曲线数据，与 5 种常用的曲率和挠率计算算法相比较，对这 6 种算法从采样密度对算法精度的影响、计算效率和抗噪声性能这 3 个方面进行了对比分析。实验结果表明，微中心差分算法总体效果最好。关键词：曲率；挠率；算法比较；离散曲线；微中心差分法；离散几何法；三维空间；差商；均匀采样中图分类号：TP311 文献标志码：A 文章编号：1673−4785(2019)01−0194−13 中文引用格式：慕生鹏, 李红军, 李世林. 三维离散曲线曲率挠率的微中心差分算法[J]. 智能系统学报, 2019, 14(1): 194–206. 英文引用格式：MU Shengpeng, LI Hongjun, LI Shilin. An algorithm for estimating curvature and torsion of discrete curve in three-dimensional space based on microcentral difference[J]. CAAI transactions on intelligent systems, 2019, 14(1): 194–206. An algorithm for estimating curvature and torsion of discrete curve in three-dimensional space based on microcentral difference MU Shengpeng，LI Hongjun，LI Shilin (College of Science, Beijing Forestry University, Beijing 100083, China) Abstract: The curvature and torsion of a 3D discrete curve reflect the degrees of its bending and distortion. To calculate these quantities accurately, following the definition of the derivative of the continuous curve, a microcentral difference algorithm, which is an extension to the one-side difference algorithm, is proposed based on the smoothing of the difference quotient. Compared with the one-side difference algorithm, the microcentral difference algorithm fails to prolong the running time but it remarkably improves the calculation accuracy. Several experiments are conducted by uniform sampling from six continuous curves, which are then compared with the five traditional algorithms of curvature and torsion. The experimental results are analyzed from three aspects: the influence of the sampling density on the accuracy of the algorithm, the efficiency of calculation, and the anti-noise performance. The experimental results show the good performance of the proposed microcentral difference algorithm. Keywords: curvature; torsion; algorithm comparison; discrete curve; microcentral difference algorithm; discrete geometry method; three-dimensional space; difference quotient; uniformly sampling 在人工智能算法设计与自动控制的相关研究中，无人机路径规划[1] ，弹道分析[2] ，公路线性设计[3] ，路径约束下的车辆行为研究[4] ，几何处理[5]和机器视觉研究[6] ，以及植物生长模拟[7] 、结构工程分析[8] 等问题的解决都可以基于离散曲线的曲率和挠率的分析。曲率和挠率是空间曲线在固有运动下的不变量[9] ，直接描述了曲线在一点邻近的形状。其中，曲率揭示曲线在所在平面的弯曲程度，挠率则刻画曲线离开既定平面的扭曲程度。也就是收稿日期：2018−02−05. 网络出版日期：2018−09−29. 基金项目：国家自然科学基金项目 (61372190)；中央高校基本科研业务费专项资金项目 (2015ZCQ-LY-01). 通信作者：李红军. E-mail：lihongjun69@bjfu.edu.cn. 第 14 卷第 1 期智能系统学报 Vol.14 No.1 2019 年 1 月 CAAI Transactions on Intelligent Systems Jan. 2019

第1期慕生鹏，等：三维离散曲线曲率挠率的微中心差分算法 ·195· 说，三维空间中的曲线是可以通过一根直线弯曲 ),()中要求3个函数x()、y()、()都是连续 (曲率)和扭曲（挠率）得到的。显然，曲率和挠率函数，且都3阶可导。能完全刻画曲线的局部行为。因此，有研究者但是对于没有方程的三维离散曲线的曲率挠根据曲线的曲率挠率特性，对空间中的曲线进行率计算，研究人员提出的各种方法中有代表性的分析1)、分类4或判断轨迹线的形状1。对于方法有：多项式曲线拟合法、B样条曲线拟合法、有解析或者参数方程的曲线，经典微分几何中有投影法、单侧差分法、离散几何法5种。本文在成熟方法计算曲线上任意一点的曲率和挠率。分析这5种算法的基础上，提出了微中心差分法。因为近些年应用日益广泛的离散空间曲线没 11多项式曲线拟合法有解析方程，经典微分几何中的求解方法不能直 2005年，Cazals等16提出了一种基于二次多接应用于离散曲线的情形，所以三维空间中离散项式曲线拟合(polynomial curve fitting,PCF)的方曲线的曲率和挠率计算成为近些年的热点研究问法来计算离散平面曲线在给定点P,的曲率、挠题。此外，在实际应用中，受人工、仪器、环境等率。该方法对点P,及其邻域(Pi-n≤k≤i+n(n为诸多因素影响，得到的离散曲线的数据（一个序点P,一侧选取的邻近点个数)，在以P:为原点的局列的三维坐标点)常常存在一定的偏差或者噪部坐标系中拟合一个二次多项式曲线。二次多项声，因此如何使算法更加准确和鲁棒，成为三维式曲线的数学模型为空间离散曲线的曲率和挠率计算的难点。 f(x)=Ao+Ax+Azx2 近些年提出的关于三维空间离散曲线的曲率通过最小二乘拟合，可求解该二次多项式曲和挠率计算方法大致可以分为2类。1)基于曲线线及其系数A,(i=0,1,2),利用多项式系数A可计算拟合的方法，这类算法主要是对给定点及其近邻出点P的曲率k(P)。类似的可推广拟合三次或三点拟合一段光滑曲线，然后根据这段解析曲线的次以上的多项式曲线，从而利用式(1)、式(2)计微分几何公式计算给定点的曲率和挠率。其本质算出三维空间离散曲线的曲率、挠率。就是用拟合曲线的曲率和挠率来估算离散曲线的 1.2B样条曲线拟合法曲率和挠率。这类算法的典型算法有多项式拟合 2017年，赵夫群等1用4次B样条曲线拟合法(polynomial curve fitting)Io和B样条曲线拟合 (B-spline curve fitting,BSCF)方法将空间离散曲线法(B-spline curve fitting).19.2)基于离散结构的拟合成光滑的空间曲线，从而以此来计算离散曲方法，这类算法无需对离散曲线进行光滑曲线的线的曲率和挠率。B样条曲线的定义为：给定拟合，直接利用离散点之间的距离和内接多边形 n=m+k+1个顶点，则可以定义m+1段k次参数曲等离散结构，运用差分形式或者离散逼近等策略线第段B样条曲线函数可表示为估算出离散曲线上各个点的曲率和挠率。基于离散结构的方法中比较典型的有投影法(projection)P、 e0=∑PNad 单侧差分法(one-sided difference)2和离散几何法 =0 (discrete geometry)2。本文根据导数定义以及单式中：i=0,1,…,n;{P40,P41,…,P4t}为定义第i段侧差分法的设计思路，融合差商平滑策略，提出曲线特征多边形的k+1个顶点；Nu(①为B样条基了微中心差分(microcentral difference)算法。底函数，表示为 t-l 1离散曲线的曲率挠率算法 N(0= 1.(-1yct+k-l- j=0 经典微分几何中，设曲线C的参数方程为对于离散曲线上的每一个离散点P,选取点 r=r(0=(x(),y(0,z()(a≤t≤)，若一阶导函数 P,及其邻域{Pli-n≤k≤i+(n为点P,一侧选取的 (0在区间[α，B上恒不为零，即参数化曲线r()是邻近点个数)内的点进行样条曲线拟合。类似的正则的，则曲率k(0)的计算公式2到为方法有：喻德生等8采用的3次均匀B样条进行曲线拟合，Kehtarnavaz等I采用的5次B样条进 K)= r(t)xr"(t) (1) r()l' 行曲线拟合，从而利用经典微分几何的式(1)、式挠率x()的计算公式为 (2)求解方法对离散曲线进行曲率、挠率的计算。 )=.r).r() 1.3投影法 Ir(t)xr"(t) (2) Mardia在1999年提出一种基于投影(projec 式(1)、式(2)中的符号“×”表示叉积；式(2)中 tion)的算法20，对于三维空间中的离散曲线，将的分子表示3个微分向量的混合积。()=(x(), 点投影到单位切线方向上。对于由n个点构成的

说，三维空间中的曲线是可以通过一根直线弯曲 (曲率) 和扭曲 (挠率) 得到的。显然，曲率和挠率能完全刻画曲线的局部行为[10]。因此，有研究者根据曲线的曲率挠率特性，对空间中的曲线进行分析[11-13] 、分类[14]或判断轨迹线的形状[15]。对于有解析或者参数方程的曲线，经典微分几何中有成熟方法计算曲线上任意一点的曲率和挠率。因为近些年应用日益广泛的离散空间曲线没有解析方程，经典微分几何中的求解方法不能直接应用于离散曲线的情形，所以三维空间中离散曲线的曲率和挠率计算成为近些年的热点研究问题。此外，在实际应用中，受人工、仪器、环境等诸多因素影响，得到的离散曲线的数据 (一个序列的三维坐标点) 常常存在一定的偏差或者噪声，因此如何使算法更加准确和鲁棒，成为三维空间离散曲线的曲率和挠率计算的难点。近些年提出的关于三维空间离散曲线的曲率和挠率计算方法大致可以分为 2 类。1) 基于曲线拟合的方法，这类算法主要是对给定点及其近邻点拟合一段光滑曲线，然后根据这段解析曲线的微分几何公式计算给定点的曲率和挠率。其本质就是用拟合曲线的曲率和挠率来估算离散曲线的曲率和挠率。这类算法的典型算法有多项式拟合法 (polynomial curve fitting)[16]和 B 样条曲线拟合法 (B-spline curve fitting)[17-19]。2) 基于离散结构的方法，这类算法无需对离散曲线进行光滑曲线的拟合，直接利用离散点之间的距离和内接多边形等离散结构，运用差分形式或者离散逼近等策略估算出离散曲线上各个点的曲率和挠率。基于离散结构的方法中比较典型的有投影法 (projection)[20] 、单侧差分法 (one-sided difference)[21]和离散几何法 (discrete geometry)[22]。本文根据导数定义以及单侧差分法的设计思路，融合差商平滑策略，提出了微中心差分 (microcentral difference) 算法。 1 离散曲线的曲率挠率算法 C r = r(t) = (x (t), y (t),z(t)) (α ⩽ t ⩽ β) r ′ (t) [ α, β] r(t) κ(t) 经典微分几何中，设曲线的参数方程为，若一阶导函数在区间上恒不为零，即参数化曲线是正则的，则曲率的计算公式[23]为 κ(t) = |r ′ (t)× r ′′(t)| |r ′ (t)| 3 (1) 挠率τ(t) 的计算公式[16]为 τ(t) = (r ′ (t),r ′′(t),r ′′′(t)) |r ′ (t)× r ′′(t)| 2 (2) 式 (1)、式 (2) 中的符号“×”表示叉积；式 (2) 中的分子表示 3 个微分向量的混合积。r(t)=(x(t)， y(t)，z(t)) 中要求 3 个函数 x(t)、y(t)、z(t) 都是连续函数，且都 3 阶可导。但是对于没有方程的三维离散曲线的曲率挠率计算，研究人员提出的各种方法中有代表性的方法有：多项式曲线拟合法、B 样条曲线拟合法、投影法、单侧差分法、离散几何法 5 种。本文在分析这 5 种算法的基础上，提出了微中心差分法。 1.1 多项式曲线拟合法 Pi Pi {Pk |i−n ⩽ k ⩽ i+n} n Pi Pi 2005 年，Cazals 等 [16]提出了一种基于二次多项式曲线拟合 (polynomial curve fitting，PCF) 的方法来计算离散平面曲线在给定点的曲率、挠率。该方法对点及其邻域 ( 为点一侧选取的邻近点个数)，在以为原点的局部坐标系中拟合一个二次多项式曲线。二次多项式曲线的数学模型为 f(x) = A0 + A1 x+ A2 x 2 Ai(i = 0,1,2) Ai Pi κ(Pi) 通过最小二乘拟合，可求解该二次多项式曲线及其系数，利用多项式系数可计算出点的曲率。类似的可推广拟合三次或三次以上的多项式曲线，从而利用式 (1)、式 (2) 计算出三维空间离散曲线的曲率、挠率。 1.2 B 样条曲线拟合法 n = m+k+1 m+1 k i 2017 年，赵夫群等[17]用 4 次 B 样条曲线拟合 (B-spline curve fitting，BSCF) 方法将空间离散曲线拟合成光滑的空间曲线，从而以此来计算离散曲线的曲率和挠率。B 样条曲线的定义为：给定个顶点，则可以定义段次参数曲线第段 B 样条曲线函数可表示为 Q(t) = ∑k l=0 Pi+lNl,k(t) i = 0,1,··· ,n {Pi+0,Pi+1,··· ,Pi+k} i k+1 Nl,k(t) 式中：；为定义第段曲线特征多边形的个顶点；为 B 样条基底函数，表示为 Nl,k(t) = 1 k! ∑k−l j=0 (−1)jC j k+1 (t+k−l− j) k Pi Pi {Pk |i−n ⩽ k ⩽ i+n} n Pi 对于离散曲线上的每一个离散点，选取点及其邻域 ( 为点一侧选取的邻近点个数) 内的点进行样条曲线拟合。类似的方法有：喻德生等[18]采用的 3 次均匀 B 样条进行曲线拟合，Kehtarnavaz 等 [19]采用的 5 次 B 样条进行曲线拟合，从而利用经典微分几何的式 (1)、式 (2) 求解方法对离散曲线进行曲率、挠率的计算。 1.3 投影法 n Mardia 在 1999 年提出一种基于投影 (projection) 的算法[20] ，对于三维空间中的离散曲线，将点投影到单位切线方向上。对于由个点构成的第 1 期慕生鹏，等：三维离散曲线曲率挠率的微中心差分算法 ·195·

·196· 智能系统学报第14卷离散曲线，其点P处的单位切向量定义为上述离散导数也称为点P的一阶离散导数。 T-PF i i=3,4,…,n-2 类似的在点P处的二阶离散导数为 t挡容易计算点P处的单位法向量N,利用T和 -xy-y) N,做叉积得到副法向量B,从而点P处的局部投影 =i- y"= 平面离散曲线定义为 (T&-N:T&.Ba ￥=7TT k=i-2.…,i+2 同理，在点P处的n阶离散导数为然后以点Y,为中心，拟合平面二次曲线 ar2+bx,从而三维空间离散曲线点P,的曲率k(P)、（一0- 挠率x(P)分别为 y=E K(P)= lY;-Y IIP:-P-il 2a 上述各阶离散导数是对连续可微函数在点 t(P)= Γ(b2+)P) P处的各阶导数的逼近。对于三维空间中的离散 1.4单侧差分法曲线，为了计算其离散导数，首先要对离散曲线 2016年Blankenburg等2通过利用一般参数进行参数化。可采用累积弦长参数化方法2来确方程的曲线导数的有限差商逼近，直接从离散点定离散点的参数值。在该方法中，离散曲线山在各中采用单侧差分(one-sided difference,.OSD)逼近点P,=(化，y,乙1≤i≤)处的累积弦长参数s,可确定为计算出空间离散曲线在给定点P的一阶差商、二阶差商和三阶差商，然后代入曲率和挠率的解析 ∑Ip1-pl 计算公式进行求解。在该方法中，定义差分为 S:= 2≤i≤n △f)=ft+k△)-ft+(k-1)△)，k=1,2,3 Pu-p 则一阶导数、二阶导数分别为 k=1 并且s1=0,从而离散曲线a为 f(t)=lim f(t+△)-f( △t 2=0 a={P}=(x(s(s,z(s)} f0=mf0+2a）-2f+a+f0_ 式中：，∈{s1,2,,Sn:坐标函数x=x(s)、y=(s)和 42 =(s)就是弦长s的离散函数。因此，平面上的离 ,△f()△fi() mAe 散导数的计算方法可以转化为用于计算三维空间 42 同理，三阶导数为中的离散曲线上各点的导数。从而，采用上述方广0=四60-20+9 法，对三维空间中的离散曲线山上的点P,其邻域 △3 △3 为{Pi-m≤k≤i+n1(n1为点P一侧选取的邻近点从而根据曲线曲率、挠率的微分定义，将上个数)，则点P处的离散导数(x(s),y(s),x(s)分别为述各阶导直接代入曲率、挠率的计算式(1)、式 t四 (2)即可求得各点处的曲率、挠率。 (si-si)(xi-xi) 1.5离散几何法 x(S)= An等2在2011年提出一种称为离散几何法 (discrete geometry)的方法，该方法基于对称邻域 - 的最小二乘的导数计算过程，可以归结为：利用 ∑5-sw-w 离散导数的定义，采用约束最优化问题的求解方 y(s)= 法来计算指定点在其对称邻域内的离散导数。对于平面内的离散曲线，点P:=(xy)及其邻域{Pi-n,≤k≤i+n1(m1为点P,一侧选取的邻近点个数)，则点P处的离散导数为 2s-s6-a z(S)= ,(x-x)0y-) k=i-n y= u- (-x)2 上述离散导数也称为点P的一阶离散导数。类似的可计算点P,的二阶、三阶离散导数。从而

离散曲线，其点 Pi处的单位切向量定义为 Ti= PiPi+1 |PiPi+1| , i = 3,4,··· ,n−2 Pi Ni Ti Ni Bi Pi 容易计算点处的单位法向量，利用和做叉积得到副法向量，从而点处的局部投影平面离散曲线定义为 Yk = ( Tk · Ni Tk ·Ti , Tk · Bi Tk ·Ti ) , k = i−2,· · ·,i+2 Yi ax2 +bx Pi κ(Pi) τ(Pi) 然后以点为中心，拟合平面二次曲线，从而三维空间离散曲线点的曲率、挠率分别为 κ(Pi) = ∥Yi −Yi−1∥ ∥Pi − Pi−1∥ τ(Pi) = 2a (b 2 +1)3/2 κ(Pi) 1.4 单侧差分法 Pi 2016 年 Blankenburg 等 [21]通过利用一般参数方程的曲线导数的有限差商逼近，直接从离散点中采用单侧差分 (one-sided difference，OSD) 逼近计算出空间离散曲线在给定点的一阶差商、二阶差商和三阶差商，然后代入曲率和挠率的解析计算公式进行求解。在该方法中，定义差分为 ∆fk(t) = f(t+k∆t)− f(t+(k−1)∆t) , k = 1,2,3 则一阶导数、二阶导数分别为 f ′ (t) = lim ∆t→0 f(t+∆t)− f(t) ∆t = lim ∆t→0 ∆f1(t) ∆t f ′′(t) = lim ∆t→0 f(t+2∆t)−2 f(t+∆t)+f(t) ∆t 2 = lim ∆t→0 ( ∆f2(t) ∆t 2 − ∆f1(t) ∆t 2 ) 同理，三阶导数为 f ′′′(t) = lim ∆t→0 ( ∆f3(t) ∆t 3 −2 ∆f2(t) ∆t 3 + ∆f1(t) ∆t 3 ) 从而根据曲线曲率、挠率的微分定义，将上述各阶导直接代入曲率、挠率的计算式 (1)、式 (2) 即可求得各点处的曲率、挠率。 1.5 离散几何法 An 等 [22]在 2011 年提出一种称为离散几何法 (discrete geometry) 的方法，该方法基于对称邻域的最小二乘的导数计算过程，可以归结为：利用离散导数的定义，采用约束最优化问题的求解方法来计算指定点在其对称邻域内的离散导数。 Pi = (xi , yi) {Pk |i−n1 ⩽ k ⩽ i+n1} n1 Pi Pi 对于平面内的离散曲线，点及其邻域 ( 为点一侧选取的邻近点个数)，则点处的离散导数为 y ′ i = ∑i+n1 k=i−n1 (xk − xi)(yk −yi) ∑i+n1 k=i−n1 (xk − xi) 2 Pi Pi 上述离散导数也称为点的一阶离散导数。类似的在点处的二阶离散导数为 yi ′′ = ∑i+n2 k=i−n2 (xk − xi)(y ′ k −y ′ i ) ∑i+n2 k=i−n2 (xk − xi) 2 同理，在点 Pi处的n阶离散导数为 y (n) i = ∑i+n k=i−n (xk − xi)(y (n−1) k −y (n−1) i ) ∑i+n k=i−n (xk − xi) 2 Pi ψd Pi= (xi , yi ,zi) 1 ⩽ i ⩽ n si 上述各阶离散导数是对连续可微函数在点处的各阶导数的逼近。对于三维空间中的离散曲线，为了计算其离散导数，首先要对离散曲线进行参数化。可采用累积弦长参数化方法[24]来确定离散点的参数值。在该方法中，离散曲线在各点 ( ) 处的累积弦长参数可确定为 si = ∑i−1 k=1 |pk+1 − pk | ∑n−1 k=1 |pk+1 − pk | , 2 ⩽ i ⩽ n 并且 s1 = 0 ，从而离散曲线 ψd为 ψd = {Pi} = {(x(si), y(si),z(si))} si ∈ {s1,s2,· · ·,sn} xi = x(si) yi=y(si) zi=z(si) si ψd Pi {Pk |i−n1 ⩽ k ⩽ i+n1} n1 Pi Pi (x ′ (si), y ′ (si),z ′ (si)) 式中：；坐标函数、和就是弦长的离散函数。因此，平面上的离散导数的计算方法可以转化为用于计算三维空间中的离散曲线上各点的导数。从而，采用上述方法，对三维空间中的离散曲线上的点，其邻域为 ( 为点一侧选取的邻近点个数)，则点处的离散导数分别为    x ′ (si) = ∑i+n1 j=i−n1 (sj − si)(xj − xi) ∑i+n1 j=i−n1 (sj − si) 2 y ′ (si) = ∑i+n1 j=i−n1 (sj − si)(yj −yi) ∑i+n1 j=i−n1 (sj − si) 2 z ′ (si) = ∑i+n1 j=i−n1 (sj − si)(zj −zi) ∑i+n1 j=i−n1 (sj − si) 2 Pi Pi 上述离散导数也称为点的一阶离散导数。类似的可计算点的二阶、三阶离散导数。从而 ·196· 智能系统学报第 14 卷

第1期慕生鹏，等：三维离散曲线曲率挠率的微中心差分算法 ·197· 利用式(1)、式(2)即可求得三维空间离散曲线的式中k2为邻域半径。曲率、挠率。同理，其三阶前向差商定义为 1.6微中心差分法单侧差分法2在进行离散曲线的曲率挠率计 g0=1f"0-f- △t 算时，会对给定的曲率挠率产生一定的偏移或者三阶后向差商定义为缩放，为了克服由单侧差分法引起的这种计算偏差，结合文献[22]中的思想和导数定义，把其单侧 f"0=员 f"(t+iA)-f"() At 差分法修改为对目标离散点的双侧差商取平均的从而其三阶差商为方法来得到目标离散点的各阶导数，并参考文献[25]中的性能度量算法思想和差分方式，将改 "0=2(g"0*"02k,A"+kA)-f"-kA》进的方法命名为微中心差分法(microcentral differ- (5) 式中k3为邻域半径。 ence algorithm). 因此对于微中心差分法，为了得到离散曲线首先对微中心差分法名称中的“微”和“中心” 上给定的目标离散点的曲率值，根据本文的计算作简单的说明。“微”是指在局部范围的逐阶差商方法，则需要知道目标离散点邻域内的离散点数计算。比如，先用改进的离散导数定义计算一阶 N,=2(k1+k2)+1,即目标离散点两侧各需k1+k2)个差商，如式(3)：再利用一阶差商值计算离散曲线离散点；而对于目标离散点的挠率计算，则至少的二阶差商，如式(4)：最后利用二阶差商值计算需要知道目标离散点邻域内的离散点数N, 三阶差商，如式（⑤）。“中心”指在目标点邻域内，进行差分计算时，以离散曲线的目标点为中心进 2k1+2+k3)+1,其中目标离散点两侧各(k1+k2+k) 个离散点，k1、k、k均为大于等于1的整数。根行差分计算。另外，方法中涉及求目标点的各阶差商，我们指定，每阶差商目标点单侧需要上一据式(3)（⑤）得到其各阶差商后，则可根据曲线曲阶差商值的个数，记为k(i=1,2,3),称为邻域半率、挠率的微分定义，将上述各阶导数式(3)~ 径。根据上述说明，微中心差分算法的计算方法 (⑤)代入曲率、挠率的计算式(1)、式(2)，即可求得具体如下。定义其一阶前向差商为离散曲线上各点处的曲率、挠率。这里提出的微 1点f0-ft-i) 中心差分计算方法是对单侧差分方法的一个依据 △t 导数定义进行的一个扩展，理论上是符合连续函数条件下可导的要求。一阶后向差商为总体来看，上述6种方法分为2种思路：1)一 1月f+i)-f@ f()= 对近邻点进行曲线拟合，它们之间的区别在于近 At 邻点数的确定以及拟合曲线模型的假设；2)对近从而其一阶差商为邻点进行离散求导，各方法之间的区别在于取点 f0=0+o)2台 1ft+i)-ft-i△ 的策略以及如何利用近邻点的信息得到各阶导函 △ 数的可靠逼近，从而得到不同的曲率和挠率的计 (3) 算方法。式中k为邻域半径。其二阶前向差商定义为 2实验比较 f0= 1f-fa-i) 本节通过实验对上述6种算法计算的曲率和挠率的计算精度进行评估。为了便于比较，本文二阶后向差商定义为中实验方法的k1、k2、k均取值为1。算法评估指 1台f+i)-f@ 标参考文献[26]的评价指标，选取最大相对误差 0台 △t (max relative error,MRE)和均方根误差(root mean square error,.RMSE),即对于获得的离散曲线，在从而其二阶差商定义为每一个离散点P,处的偏差定义为。=9，-，其中 )-A@kA-F(-kA 9为实验方法的计算值，为对应的准确值。因 ) 此，对于给定的离散曲线，其最大相对误差为EE=

利用式 (1)、式 (2) 即可求得三维空间离散曲线的曲率、挠率。 1.6 微中心差分法单侧差分法[22]在进行离散曲线的曲率挠率计算时，会对给定的曲率挠率产生一定的偏移或者缩放，为了克服由单侧差分法引起的这种计算偏差，结合文献[22]中的思想和导数定义，把其单侧差分法修改为对目标离散点的双侧差商取平均的方法来得到目标离散点的各阶导数，并参考文献[25]中的性能度量算法思想和差分方式，将改进的方法命名为微中心差分法 (microcentral difference algorithm)。 ki(i = 1,2,3) 首先对微中心差分法名称中的“微”和“中心” 作简单的说明。“微”是指在局部范围的逐阶差商计算。比如，先用改进的离散导数定义计算一阶差商，如式 (3)；再利用一阶差商值计算离散曲线的二阶差商，如式 (4)；最后利用二阶差商值计算三阶差商，如式 (5)。“中心”指在目标点邻域内，进行差分计算时，以离散曲线的目标点为中心进行差分计算。另外，方法中涉及求目标点的各阶差商，我们指定，每阶差商目标点单侧需要上一阶差商值的个数，记为，称为邻域半径。根据上述说明，微中心差分算法的计算方法具体如下。定义其一阶前向差商为 f ′ q (t) = 1 k1 ∑k1 i=1 f(t)− f(t−i∆t) ∆t 一阶后向差商为 f ′ h (t) = 1 k1 ∑k1 i=1 f(t+i∆t)− f(t) ∆t 从而其一阶差商为 f ′ (t) = 1 2 ( f ′ q (t)+ f ′ h (t) ) = 1 2k1 ∑k1 i=1 f(t+i∆t)− f(t−i∆t) ∆t (3) 式中 k1为邻域半径。其二阶前向差商定义为 f ′′ q (t) = 1 k2 ∑k2 i=1 f ′ (t)− f ′ (t−i∆t) ∆t 二阶后向差商定义为 f ′′ h (t) = 1 k2 ∑k2 i=1 f ′ (t+i∆t)− f ′ (t) ∆t 从而其二阶差商定义为 f ′′(t)= 1 2 ( f ′′ q (t)+ f ′′ h (t) ) = 1 2k2∆t ( f ′ (t+k2∆t)− f ′ (t−k2∆t) ) (4) 式中 k2为邻域半径。同理，其三阶前向差商定义为 f ′′′ q (t) = 1 k3 ∑k3 i=1 f ′′(t)− f ′′(t−i∆t) ∆t 三阶后向差商定义为 f ′′′ h (t) = 1 k3 ∑k3 i=1 f ′′(t+i∆t)− f ′′(t) ∆t 从而其三阶差商为 f ′′′(t)= 1 2 ( f ′′′ q (t)+ f ′′′ h (t) ) = 1 2k3∆t (f ′′(t+k2∆t)− f ′′(t−k2∆t)) (5) 式中 k3为邻域半径。 N1 = 2(k1 +k2)+1 (k1 +k2) (k1 +k2 +k3)+1 (k1 +k2 +k3) 因此对于微中心差分法，为了得到离散曲线上给定的目标离散点的曲率值，根据本文的计算方法，则需要知道目标离散点邻域内的离散点数，即目标离散点两侧各需个离散点；而对于目标离散点的挠率计算，则至少需要知道目标离散点邻域内的离散点数 N2= 2 ，其中目标离散点两侧各个离散点，k1、k2、k3 均为大于等于 1 的整数。根据式 (3)~(5) 得到其各阶差商后，则可根据曲线曲率、挠率的微分定义，将上述各阶导数式 (3)~ (5) 代入曲率、挠率的计算式 (1)、式 (2)，即可求得离散曲线上各点处的曲率、挠率。这里提出的微中心差分计算方法是对单侧差分方法的一个依据导数定义进行的一个扩展，理论上是符合连续函数条件下可导的要求。总体来看，上述 6 种方法分为 2 种思路：1) 一对近邻点进行曲线拟合，它们之间的区别在于近邻点数的确定以及拟合曲线模型的假设；2) 对近邻点进行离散求导，各方法之间的区别在于取点的策略以及如何利用近邻点的信息得到各阶导函数的可靠逼近，从而得到不同的曲率和挠率的计算方法。 2 实验比较 Pi e i φ=φi−φi φi φi EMRE = 本节通过实验对上述 6 种算法计算的曲率和挠率的计算精度进行评估。为了便于比较，本文中实验方法的 k1、k2、k3 均取值为 1。算法评估指标参考文献[26]的评价指标，选取最大相对误差 (max relative error，MRE) 和均方根误差 (root mean square error，RMSE)，即对于获得的离散曲线，在每一个离散点处的偏差定义为，其中为实验方法的计算值，为对应的准确值。因此，对于给定的离散曲线，其最大相对误差为第 1 期慕生鹏，等：三维离散曲线曲率挠率的微中心差分算法 ·197·

·198· 智能系统学报第14卷 max(p,-)/pl)均方根误差为ERsE= 3)循环曲线：最大相对误差体现算法的稳定性，均方根误差反 r(t)=(cos(t),sin(t),2sin(t),t∈[0,2π] 映算法的整体有效性。实验平台为个人笔记本电 4)Viviani曲线： r(t)=(cos2(t),cos(t)sin(t),sin(t)),tE[0,2] 脑，计算机配置是Inter(R)Core(TM)i5-7300HQ 5)Clelia曲线： CPU@2.50GHz处理器和8GB内存，程序运行环 r(t)=(cos(2t)cos(t),cos(2r)sin(t),sin(2r)),t[0,2] 境是MATLAB R2014b。 6)环面螺线： 2.1实验数据 r(0=(2+cos(t)cos(2),(2+cos(t)sin(2r),),t∈[0,2 实验用如图1所示的6条连续可微参数曲线需要说明的是，解析曲线方程用于生成离散进行分析。由于有曲线方程，可以用经典微分几何曲线，并用于计算离散曲线上每个采样点的准确的方法（式(1）、式(2)计算出曲线上任意一个点的曲率值和挠率值。实验中对每条空间参数曲线的准确曲率值和挠率值，有了准确值便于进行误通过均匀采样获得相应的离散曲线，采样点数差分析。实验分析用到的6条曲线参数方程如下： n在定义域内从50个点变化到500个点，用于测 1)三次挠曲线：试算法的有效性。第1章介绍的5种常用的离散 r(t0=(2t2,f),te[-2,2] 曲线曲率和挠率计算方法以及本文提出的算法都 2)多变式内接弹簧线：采用离散曲线上的点的坐标进行计算，都用不到 r(t)=(2cos(3t)+cos(2t),2sin(3t)+sin(21),t),tE [0,2n] 曲线方程。 10 14 1.2 10r 0.18 0 1.0 0.16 0.8 15 0.14 0.6 0.12 -10 4 0.4 0.10 0-4202 0.2 0.08 54-20 0.06 (a)三次挠曲线 (b)多变式内接弹簧线 2 1.0 0.7 0.8 0.6 0 0.5 8 0 1.0 0.4 -1 05 0.5 -0.5 0.2 1.0 0.2 -0.5 0.5 0.1 -10 0 0.50 (©)循环曲线 (d)Viviani曲线 ■0.7 0.6 0.5 1 0.20 0 0 0.4 -1 0.15 0.3 0.10 0.5 1.0 0.2 .0 0.1 -2 0.05 -1.0 0.50 (①环面螺线 (e)Clelia曲线图1用于实验的6条空间曲线（线颜色与挠率值对应） Fig.1 Six 3D curves used in the experiment (The color is corresponding to torsion value) 2.2误差分析及时间效率比较线拟合法(PCF)、投影法(Proj)、单侧差分法通过上述三维空间6条离散曲线来定量地评 (OSD)、离散几何法(DG)这5种算法的曲率和挠价本文提出的微中心差方法(MCD),并将其与第率的计算结果进行比较。图2~7分别给出了这 1章介绍的B样条曲线拟合法(BSCF)、多项式曲 6条离散曲线随着采样密度增加，曲率和挠率的

max(|(φi−φi)/φi |) ERMSE= √ 1 n ∑n i=1 (φi−φi) 2 , 均方根误差为。最大相对误差体现算法的稳定性，均方根误差反映算法的整体有效性。实验平台为个人笔记本电脑，计算机配置是 Inter(R) Core(TM) i5-7300HQ CPU @2.50 GHz 处理器和 8 GB 内存，程序运行环境是 MATLAB R2014b。 2.1 实验数据实验用如图 1 所示的 6 条连续可微参数曲线进行分析。由于有曲线方程，可以用经典微分几何的方法 (式 (1)、式 (2)) 计算出曲线上任意一个点的准确曲率值和挠率值，有了准确值便于进行误差分析。实验分析用到的 6 条曲线参数方程如下： 1) 三次挠曲线： r(t) = (2t,t 2 ,t 3 ), t ∈ [−2,2] 2) 多变式内接弹簧线： r(t) = (2 cos(3t)+cos(2t),2 sin(3t)+sin(2t),t), t ∈ [0,2π] 3) 循环曲线： r(t) = (cos(t),sin(t),2 sin(t)), t ∈ [0,2π] 4) Viviani 曲线： r(t) = (cos2 (t), cos(t) sin(t),sin(t)), t ∈ [0,2π] 5) Clelia 曲线： r(t) = (cos(2t) cos(t), cos(2t) sin(t),sin(2t)), t ∈ [0,2π] 6) 环面螺线： r(t) = ((2+cos(t)) cos(2t),(2+cos(t)) sin(2t),t), t ∈ [0,2π] 需要说明的是，解析曲线方程用于生成离散曲线，并用于计算离散曲线上每个采样点的准确的曲率值和挠率值。实验中对每条空间参数曲线通过均匀采样获得相应的离散曲线，采样点数 n 在定义域内从 50 个点变化到 500 个点，用于测试算法的有效性。第 1 章介绍的 5 种常用的离散曲线曲率和挠率计算方法以及本文提出的算法都采用离散曲线上的点的坐标进行计算，都用不到曲线方程。 2.2 误差分析及时间效率比较通过上述三维空间 6 条离散曲线来定量地评价本文提出的微中心差方法 (MCD)，并将其与第 1 章介绍的 B 样条曲线拟合法 (BSCF)、多项式曲线拟合法 (PCF)、投影法 (Proj)、单侧差分法 (OSD)、离散几何法 (DG) 这 5 种算法的曲率和挠率的计算结果进行比较。图 2～7 分别给出了这 6 条离散曲线随着采样密度增加，曲率和挠率的 (a) 三次挠曲线 1.4 1.2 1.0 0.8 0.6 0.4 0.2 10 −10 −2 0 2 4 −4 z 0 4 2 0 y x (b) 多变式内接弹簧线 0.18 0.16 0.14 0.12 0.10 0.08 0.06 10 5 0 z y x 5 0 −5 −4 −2 0 2 4 (c) 循环曲线 1.0 0.8 0.6 0.4 0.2 0 z 2 1 0 1.0 1.0 0.5 0.5 0 −0.5 0 −0.5 −1.0 y x (d) Viviani 曲线 0.7 0.6 0.5 0.4 0.3 0.2 0.1 z 1 0 −1 0.5 0 −0.5 y 0 0.5 1.0 x (e) Clelia 曲线 0.7 0.6 0.5 0.4 0.3 0.2 0.1 z 1 −1 0 1 0 −1.0 −0.5 0 0.5 1.0 y x y x (f) 环面螺线 0.20 0.15 0.10 0.05 z 1 −1 0 4 4 2 2 0 −2 0 −2 −4 图 1 用于实验的 6 条空间曲线 (线颜色与挠率值对应) Fig. 1 Six 3D curves used in the experiment (The color is corresponding to torsion value) ·198· 智能系统学报第 14 卷

第1期慕生鹏，等：三维离散曲线曲率挠率的微中心差分算法 ·199· 最大相对误差(MRE)与均方根误差(RMSE)的变看，本文提出的算法在所有离散曲线的计算结果化。从图2~7的误差比较图可以看出，随着采样中几乎都能取得最佳效果。密度的不断增加，这6种方法计算的曲率和挠率图2~7分别给出了6条离散曲线随着采样密的最大相对误差和均方根误差都逐渐地减小，并度增加的挠率最大相对误差(MRE)。从图中可以且采样密度越高，最大相对误差值和均方根误差看出，随着采样密度的不断增加，本文方法的挠值越接近于0，这表明这6种算法都能有效地计率最大相对误差逐渐地减小，并且采样密度越算出三维空间离散曲线的曲率和挠率。大，最大相对误差值和均方根误差值就越接近于当然，这6种方法在收敛速度和稳定性方面 0,这表明挠率计算越准确，也表明本文方法在计存在一些差异。从图2~7的曲率MRE对比子图算三维空间离散曲线的挠率时收敛性和稳定性较可以看出：投影法和单侧差分法的最大误差最好。而单侧差分法、投影法和B样条曲线拟合法大，且随着采样密度增加，误差一直比较大；离散的最大相对误差有时会有较大波动。几何法与本文提出的微中心差分法的误差较小，图2~7中(c)图呈现的挠率计算的均方根误特别是当采样较稀疏时，曲率估算的最大相对误差变化结果表明，随着采样密度增大，本文方法差最小。从图2~7的这6个图中的每个子图与离散几何方法和多项式拟合方法的均方根误差 (©)图的结果，也说明从均方根指标评价，也能得能单调递减收敛于0，而单侧差分法、投影法的误出同样结论。需要指出的是，从均方根误差指标差收敛速度较差。 0.20 一。B样条曲线拟合法 0.030 ◆B样条曲线拟合法 ◆多项式曲线拟合法 ◆ 多项式曲线拟合法 ◆一投影法投影法单侧差分法 0.025 ◆一单侧差分法 0.15 离散几何法离散几何法。一微中心差分法 0.020 。一微中心差分法 0.10 015 里0.010 0.05 0.005 200 400 600 200 400 600 点数n 点数n (a)曲率MRE对比 (b)曲率RMSE对比 0.30r ,一B样条曲线拟合法 0.12 ·一B样条曲线拟合法 ◆一多项式曲线拟合法 ◆一多项式曲线拟合法 0.25 投影法投影法单侧差分法 0.10 单侧差分法 ·一离散几何法离散几何法 0.20 。一微中心差分法烂0.08 微中心差分法 0.15 0.10 器0.04 0.05 0.02 200 400 600 200 400 600 点数n 点数n (c)挠率MRE对比 (d)挠率RMSE对比图2三次挠曲线误差对比实验 Fig.2 Error comparison experiment of cubic deflection curve

n 最大相对误差 (MRE) 与均方根误差 (RMSE) 的变化。从图 2~7 的误差比较图可以看出，随着采样密度的不断增加，这 6 种方法计算的曲率和挠率的最大相对误差和均方根误差都逐渐地减小，并且采样密度越高，最大相对误差值和均方根误差值越接近于 0，这表明这 6 种算法都能有效地计算出三维空间离散曲线的曲率和挠率。当然，这 6 种方法在收敛速度和稳定性方面存在一些差异。从图 2~7 的曲率 MRE 对比子图可以看出：投影法和单侧差分法的最大误差最大，且随着采样密度增加，误差一直比较大；离散几何法与本文提出的微中心差分法的误差较小，特别是当采样较稀疏时，曲率估算的最大相对误差最小。从图 2~7 的这 6 个图中的每个子图 (c) 图的结果，也说明从均方根指标评价，也能得出同样结论。需要指出的是，从均方根误差指标看，本文提出的算法在所有离散曲线的计算结果中几乎都能取得最佳效果。 n 图 2~7 分别给出了 6 条离散曲线随着采样密度增加的挠率最大相对误差 (MRE)。从图中可以看出，随着采样密度的不断增加，本文方法的挠率最大相对误差逐渐地减小，并且采样密度越大，最大相对误差值和均方根误差值就越接近于 0，这表明挠率计算越准确，也表明本文方法在计算三维空间离散曲线的挠率时收敛性和稳定性较好。而单侧差分法、投影法和 B 样条曲线拟合法的最大相对误差有时会有较大波动。图 2~7 中 (c) 图呈现的挠率计算的均方根误差变化结果表明，随着采样密度增大，本文方法与离散几何方法和多项式拟合方法的均方根误差能单调递减收敛于 0，而单侧差分法、投影法的误差收敛速度较差。 (a) 曲率MRE对比 (b) 曲率RMSE对比 (c) 挠率MRE对比 (d) 挠率RMSE对比 0.20 0.15 0.10 0.05 曲率最大相对误差值 0 200 400 600 点数n 点数n 0.030 0.025 0.020 0.015 0.010 0.005 0 曲率均方根误差值 200 400 600 0.30 0.25 0.20 0.15 0.10 0.05 0 挠率最大相对误差值 200 400 600 点数n 200 400 600 点数n 0.12 0.10 0.08 0.06 0.04 0.02 0 挠率均方根误差值 B样条曲线拟合法多项式曲线拟合法投影法单侧差分法离散几何法微中心差分法 B样条曲线拟合法多项式曲线拟合法投影法单侧差分法离散几何法微中心差分法 B样条曲线拟合法多项式曲线拟合法投影法单侧差分法离散几何法微中心差分法 B样条曲线拟合法多项式曲线拟合法投影法单侧差分法离散几何法微中心差分法图 2 三次挠曲线误差对比实验 Fig. 2 Error comparison experiment of cubic deflection curve 第 1 期慕生鹏，等：三维离散曲线曲率挠率的微中心差分算法 ·199·

第1期慕生鹏，等：三维离散曲线曲率挠率的微中心差分算法 ·203· 。一B样条曲线拟合法 0.05 ·B样条曲线拟合法 ◆一多项式曲线拟合法 ◆一多项式曲线拟合法 ◆一投影法 ◆一投影法 ◆一单侧差分法 0.04 单侧差分法 ·一离散几何法 ·高散几何法。一微中心差分法过。微中心差分法 0.01 一+ 100 200300 400 500 200 400 600 点数n 点数n (c)挠率MRE对比 (d)挠率RMSE对比图7环面螺线误差对比实验 Fig.7 Error comparison experiment of toroidal spiral curve 图2~7反映了随着点的密度增加，曲率和挠速度最快，本文提出的微中心差分法处于中间率的误差的收敛性和稳定性。下面针对固定的点位置。密度，比较一下各个方法的精度。实验中我们表26种算法的平均运行时间比较根据上述的密度实验，取采样点数量适中误差 Table 2 Comparison of average running time of the 比较稳定的点密度，即n=200的时候进行6条 six algorithms 曲线的曲率和挠率估计精度的比较，结果如表1 算法 BSCF PCF Proj OSD DG MCD 所示。时间56.655963.69350.05380.00090.00700.0120 表16条曲线的平均误差比较(F200) Table 1 Comparison of average error(n=200) 23噪声数据实验算法BSCF 为了测试每个算法的抗噪声鲁棒性，在采样 PCF Proj OSD DG MCD 点数n为200时，设置噪声水平为0.00010.1，噪曲率0.01290.00590.03930.02530.00140.0004 声为均值为零、标准差为离散曲线上邻近的两离挠率0.00530.00120.03300.03020.00080.0003 散点之间的平均距离与噪声水平的乘积的高斯噪从表1可以看出，本文提出的MCD方法计算声。实验用Clelia曲线生成了不同噪声水平的离的曲率和挠率平均误差分别为0.0004和0.0003，散曲线。比另5种算法中效果最佳的DG算法的曲率和挠 6种算法计算的曲率最大相对误差和曲率均率平均误差（分别为0.0014和0.0008），分别降低方根误差分别如图8(a)、(b)所示，从图中可以看了69.74%和55.36%。可见，总体来看，本文提出出：噪声水平越高，算法的误差变化越大；各算法的算法无论是曲率估计还是挠率估计，精度都是在噪声水平不大于0.01时，均能较准确地计算离最高的。散曲线的曲率值；噪声水平大于0.01时，其误差关于运行效率，从算法原理可以看出，所有算明显变大，尤其是单侧差分法和投影法，而本文的方法误差变化浮动相对较小。法时间复杂度都是线性级的，即O(),因此给出 6种算法计算的挠率最大相对误差和均方根平均实验时间便可以说明计算的时间效率。用误差分别如图8(c)、(d)所示，从图中可以看出：含 6种算法分别计算6条曲线的曲率和挠率，对每噪声的离散曲线在计算挠率时变得更加敏感，当条离散曲线的不同采样密度进行50次实验，总噪声水平大于0.001时，6种算法得到的挠率值与的平均运行时间如表2所示。从表2可以看真实值就已有了明显的偏差，并且随着噪声水平出，离散结构法(Proj、OSD、DG、MCD)相比曲变大，挠率的最大相对误差和均方根误差相比曲线拟合法(BSCF、PCF),计算效率大幅提高；在率也会更快地增大。与其他5种算法相比，本文离散结构法中，单侧差分算法(OSD)平均运行的方法误差变化浮动相对较小

图 2～7 反映了随着点的密度增加，曲率和挠率的误差的收敛性和稳定性。下面针对固定的点密度，比较一下各个方法的精度。实验中我们根据上述的密度实验，取采样点数量适中误差比较稳定的点密度，即 n=200 的时候进行 6 条曲线的曲率和挠率估计精度的比较，结果如表 1 所示。从表 1 可以看出，本文提出的 MCD 方法计算的曲率和挠率平均误差分别为 0.000 4 和 0.000 3，比另 5 种算法中效果最佳的 DG 算法的曲率和挠率平均误差 (分别为 0.001 4 和 0.000 8)，分别降低了 69.74% 和 55.36%。可见，总体来看，本文提出的算法无论是曲率估计还是挠率估计，精度都是最高的。关于运行效率，从算法原理可以看出，所有算法时间复杂度都是线性级的，即 O(n)，因此给出平均实验时间便可以说明计算的时间效率。用 6 种算法分别计算 6 条曲线的曲率和挠率，对每条离散曲线的不同采样密度进行 50 次实验，总的平均运行时间如表 2 所示。从表 2 可以看出，离散结构法 (Proj、OSD、DG、MCD) 相比曲线拟合法 (BSCF、PCF)，计算效率大幅提高；在离散结构法中，单侧差分算法 (OSD) 平均运行速度最快，本文提出的微中心差分法处于中间位置。 2.3 噪声数据实验 n 为了测试每个算法的抗噪声鲁棒性，在采样点数为 200 时，设置噪声水平为 0.000 1~0.1，噪声为均值为零、标准差为离散曲线上邻近的两离散点之间的平均距离与噪声水平的乘积的高斯噪声。实验用 Clelia 曲线生成了不同噪声水平的离散曲线。 6 种算法计算的曲率最大相对误差和曲率均方根误差分别如图 8(a)、(b) 所示，从图中可以看出：噪声水平越高，算法的误差变化越大；各算法在噪声水平不大于 0.01 时，均能较准确地计算离散曲线的曲率值；噪声水平大于 0.01 时，其误差明显变大，尤其是单侧差分法和投影法，而本文的方法误差变化浮动相对较小。 6 种算法计算的挠率最大相对误差和均方根误差分别如图 8(c)、(d) 所示，从图中可以看出：含噪声的离散曲线在计算挠率时变得更加敏感，当噪声水平大于 0.001 时，6 种算法得到的挠率值与真实值就已有了明显的偏差，并且随着噪声水平变大，挠率的最大相对误差和均方根误差相比曲率也会更快地增大。与其他 5 种算法相比，本文的方法误差变化浮动相对较小。 (d) 挠率RMSE对比挠率均方根误差值点数n 200 400 600 0.05 0.04 0.03 0.02 0.01 0 B样条曲线拟合法多项式曲线拟合法投影法单侧差分法离散几何法微中心差分法 (c) 挠率MRE对比 200 400 500 点数n 102 101 0 挠率最大相对误差值 100 100 300 B样条曲线拟合法多项式曲线拟合法投影法单侧差分法离散几何法微中心差分法图 7 环面螺线误差对比实验 Fig. 7 Error comparison experiment of toroidal spiral curve 表 1 6 条曲线的平均误差比较 (n=200) Table 1 Comparison of average error (n = 200) 算法 BSCF PCF Proj OSD DG MCD 曲率 0.012 9 0.005 9 0.039 3 0.025 3 0.001 4 0.000 4 挠率 0.005 3 0.001 2 0.033 0 0.030 2 0.000 8 0.000 3 表 2 6 种算法的平均运行时间比较 Table 2 Comparison of average running time of the six algorithms s 算法 BSCF PCF Proj OSD DG MCD 时间 56.655 9 63.693 5 0.053 8 0.000 9 0.007 0 0.012 0 第 1 期慕生鹏，等：三维离散曲线曲率挠率的微中心差分算法 ·203·