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

第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 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 在人工智能算法设计与自动控制的相关研究 中，无人机路径规划[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 年提出一种基于投影 (projec￾tion) 的算法[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 differ￾ence 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·

·200· 智能系 统学 报 第14卷 0.08 。B样条曲线拟合法 0.035 B样条曲线拟合法 一◆一多项式曲线拟合法 ◆一 多项式曲线拟合法 ◆ 投影法 0.030 投影法 单侧差分法 单侧差分法 0.06 ·离散几何法 ·一离散几何法 。一微中心差分法 0.025 。一微中心差分法 0.04 0.015 0.010 0.02 0.005 200 400 600 0 200 400 600 点数n 点数n (a)曲率MRE对比 (b)曲率RMSE对比 0.25 ◆一B样条曲线拟合法 0.025 ◆B样条曲线拟合法 多项式曲线拟合法 ◆ 多项式曲线拟合法 投影法 投影法 0.20 ◆ 单侧差分法 0.020 单侧差分法 一离散几何法 一离散几何法 微中心差分法 微中心差分法 0.15 0.015 0.10 0.05 0.005 200 400 600 200 400 600 点数n 点数n (C)挠率MRE对比 (d)挠率RMSE对比 图3 多变式内接弹簧线误差对比实验 Fig.3 Error comparison experiment of changeable inner spring curve 1.0 ◆B样条曲线拟合法 0.7 ◆一B样条曲线拟合法 ◆一多项式曲线拟合法 ◆一多项式曲线拟合法 ·一投影法 0.6 。一投影法 0.8 一。单侧差分法 一◆一单侧差分法 一离散几何法 0.5 一离散几何法 。一微中心差分法 。一微中心差分法 0.6 0.3 0.4 0.2 0.1 200 400 600 200 400 600 点数n 点数n (a)曲率MRE对比 (b)曲率RMSE对比

(b) 曲率RMSE对比 (c) 挠率MRE对比 (d) 挠率RMSE对比 0.08 0.06 0.04 0.02 0 曲率最大相对误差值 200 400 600 (a) 曲率MRE对比 点数n 点数n 200 400 600 点数n 0.035 0.030 0.025 0.020 0.015 0.010 0.005 0 曲率均方根误差值 0.25 0.20 0.15 0.10 0.05 0 挠率最大相对误差值 200 400 600 点数n 200 400 600 0.025 0.020 0.015 0.010 0.005 0 挠率均方根误差值 B样条曲线拟合法 多项式曲线拟合法 投影法 单侧差分法 离散几何法 微中心差分法 B样条曲线拟合法 多项式曲线拟合法 投影法 单侧差分法 离散几何法 微中心差分法 B样条曲线拟合法 多项式曲线拟合法 投影法 单侧差分法 离散几何法 微中心差分法 B样条曲线拟合法 多项式曲线拟合法 投影法 单侧差分法 离散几何法 微中心差分法 图 3 多变式内接弹簧线误差对比实验 Fig. 3 Error comparison experiment of changeable inner spring curve B样条曲线拟合法 多项式曲线拟合法 投影法 单侧差分法 离散几何法 微中心差分法 B样条曲线拟合法 多项式曲线拟合法 投影法 单侧差分法 离散几何法 微中心差分法 (b) 曲率RMSE对比 1.0 0.8 0.6 0.4 0.2 0 曲率最大相对误差值 200 400 600 点数n 200 400 600 点数n 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 曲率均方根误差值 (a) 曲率MRE对比 ·200· 智 能 系 统 学 报 第 14 卷

第1期 慕生鹏，等：三维离散曲线曲率挠率的微中心差分算法 ·201· 0.30 ◆B样条曲线拟合法 ◆一 多项式曲线拟合法 6 0.25 投影法 。一B样条曲线拟合法 单侧差分法 ◆一多项式曲线拟合法 一离散几何法 。一投影法 。一单侧差分法 写0.20 。一微中心差分法 4 一离散几何法 。一微中心差分法 0.15 0.10 0.05 200 400 600 200 400 600 点数n 点数n (c)挠率MRE对比 (d)挠率RMSE对比 图4循环曲线误差对比实验 Fig.4 Error comparison experiment of circular curve 0.10 ◆B样条曲线拟合法 0.10r ◆B样条曲线拟合法 ◆多项式曲线拟合法 ◆一 多项式曲线拟合法 ◆投影法 ◆ 投影法 0.08 单侧差分法 0.08 单侧差分法 一 离散几何法 一 离散几何法 微中心差分法 。微中心差分法 0.06 0.04 0.02 0.02 200 400 600 200 400 600 点数n 点数n (a)曲率MRE对比 (b)曲率RMSE对比 0 0.14 。B样条曲线拟合法 ◆一 多项式曲线拟合法 0.12 投影法 。B样条曲线拟合法 单侧差分法 ◆一多项式曲线拟合法 0.10 。一离散几何法 10 投影法 。微中心差分法 单侧差分法 昭0.08 10r 一 离散几何法 。一微中心差分法 罩06 0.04 10 0.02 10 100 200 300 400 500 200 400 600 点数n 点数n (c)挠率MRE对比 (d挠率RMSE对比 图5 Viviani曲线误差对比实验 Fig.5 Error comparison experiment of Viviani curve

B样条曲线拟合法 多项式曲线拟合法 投影法 单侧差分法 离散几何法 微中心差分法 B样条曲线拟合法 多项式曲线拟合法 投影法 单侧差分法 离散几何法 微中心差分法 (c) 挠率MRE对比 (d) 挠率RMSE对比 点数n 7 6 5 4 3 2 1 0 挠率最大相对误差值 200 400 600 点数n 200 400 600 0.30 0.25 0.20 0.15 0.10 0.05 0 挠率均方根误差值 图 4 循环曲线误差对比实验 Fig. 4 Error comparison experiment of circular curve (a) 曲率MRE对比 0.10 0.08 0.06 0.04 0.02 0 曲率最大相对误差值 200 400 600 点数n B样条曲线拟合法 多项式曲线拟合法 投影法 单侧差分法 离散几何法 微中心差分法 (c) 挠率MRE对比 200 400 500 点数n 101 100 10−1 10−2 10−3 10−4 0 挠率最大相对误差值 (d) 挠率RMSE对比 200 400 600 点数n 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0 挠率均方根误差值 B样条曲线拟合法 多项式曲线拟合法 投影法 单侧差分法 离散几何法 微中心差分法 (b) 曲率RMSE对比 0.10 0.08 0.06 0.04 0.02 0 200 400 600 点数n 曲率均方根误差值 B样条曲线拟合法 多项式曲线拟合法 投影法 单侧差分法 离散几何法 微中心差分法 10−5 100 300 B样条曲线拟合法 多项式曲线拟合法 投影法 单侧差分法 离散几何法 微中心差分法 图 5 Viviani 曲线误差对比实验 Fig. 5 Error comparison experiment of Viviani curve 第 1 期 慕生鹏，等：三维离散曲线曲率挠率的微中心差分算法 ·201·

·202· 令 能系 统学 报 第14卷 0.08 ◆B样条曲线拟合法 0.07 ◆B样条曲线拟合法 ◆多项式曲线拟合法 ◆一多项式曲线拟合法 一 投影法 0.06 投影法 ◆单侧差分法 单侧差分法 0.06 。一离散几何法 0.05 一 离散几何法 。微中心差分法 。微中心差分法 0.04 0.04 0.03 百 0.02 0.02 0.01 200 400 600 200 400 600 点数n 点数n (a)曲率MRE对比 (b)曲率RMSE对比 102 0.25 ◆B样条曲线拟合法 ◆ 多项式曲线拟合法 投影法 0.20 单侧差分法 离散几何法 10° 微中心差分法 ◆一B样条曲线拟合法 0.15 10 ◆一 多项式曲线拟合法 投影法 单侧差分法 0.10 10- 离散几何法 微中心差分法 10 0.05 10 100 200 300 400 500 100 200 300 400500 点数n 点数n (C)挠率MRE对比 (d)挠率RMSE对比 图6 Clelia曲线误差对比实验 Fig.6 Error comparison experiment of Clelia curve 0.10 0.035 ◆一B样条曲线拟合法 B样条曲线拟合法 ◆多项式曲线拟合法 ◆一 多项式曲线拟合法 。一投影法 0.030 一投影法 0.08 。一单侧差分法 。一单侧差分法 一离散几何法 一离散几何法 。微中心差分法 0.025 。一微中心差分法 0.06 0.020 0.015 0.04 0.010 0.02 0.005 200 400 600 200 400 600 点数n 点数n (a)曲率MRE对比 (b)曲率RMSE对比

B样条曲线拟合法 多项式曲线拟合法 投影法 单侧差分法 离散几何法 微中心差分法 B样条曲线拟合法 多项式曲线拟合法 投影法 单侧差分法 离散几何法 微中心差分法 B样条曲线拟合法 多项式曲线拟合法 投影法 单侧差分法 离散几何法 微中心差分法 (d) 挠率RMSE对比 (b) 曲率RMSE对比 0.08 0.06 0.04 0.02 0 曲率最大相对误差值 200 400 600 点数n 点数n 200 400 600 点数n (a) 曲率MRE对比 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0 曲率均方根误差值 挠率均方根误差值 100 200 300 400 500 0.25 0.20 0.15 0.10 0.05 0 (c) 挠率MRE对比 200 400 500 点数n 102 101 100 10−1 10−2 10−3 0 挠率最大相对误差值 10−4 100 300 B样条曲线拟合法 多项式曲线拟合法 投影法 单侧差分法 离散几何法 微中心差分法 图 6 Clelia 曲线误差对比实验 Fig. 6 Error comparison experiment of Clelia curve 0.10 0.08 0.06 0.04 0.02 0 曲率最大相对误差值 200 400 600 点数n (a) 曲率MRE对比 B样条曲线拟合法 多项式曲线拟合法 投影法 单侧差分法 离散几何法 微中心差分法 200 400 600 点数n (b) 曲率RMSE对比 0.035 0.030 0.025 0.020 0.015 0.010 0.005 0 曲率均方根误差值 B样条曲线拟合法 多项式曲线拟合法 投影法 单侧差分法 离散几何法 微中心差分法 ·202· 智 能 系 统 学 报 第 14 卷

第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·