三維離散曲線曲率撓率的微中心差分算法

慕生鵬，李紅軍，李世林

（北京林業(yè)大學(xué) 理學(xué)院，北京 100083）

在人工智能算法設(shè)計(jì)與自動(dòng)控制的相關(guān)研究中，無人機(jī)路徑規(guī)劃[1]，彈道分析[2]，公路線性設(shè)計(jì)[3]，路徑約束下的車輛行為研究[4]，幾何處理[5]和機(jī)器視覺研究[6]，以及植物生長(zhǎng)模擬[7]、結(jié)構(gòu)工程分析[8]等問題的解決都可以基于離散曲線的曲率和撓率的分析。曲率和撓率是空間曲線在固有運(yùn)動(dòng)下的不變量[9]，直接描述了曲線在一點(diǎn)鄰近的形狀。其中，曲率揭示曲線在所在平面的彎曲程度，撓率則刻畫曲線離開既定平面的扭曲程度。也就是說，三維空間中的曲線是可以通過一根直線彎曲(曲率)和扭曲(撓率)得到的。顯然，曲率和撓率能完全刻畫曲線的局部行為[10]。因此，有研究者根據(jù)曲線的曲率撓率特性，對(duì)空間中的曲線進(jìn)行分析[11-13]、分類[14]或判斷軌跡線的形狀[15]。對(duì)于有解析或者參數(shù)方程的曲線，經(jīng)典微分幾何中有成熟方法計(jì)算曲線上任意一點(diǎn)的曲率和撓率。

因?yàn)榻┠陸?yīng)用日益廣泛的離散空間曲線沒有解析方程，經(jīng)典微分幾何中的求解方法不能直接應(yīng)用于離散曲線的情形，所以三維空間中離散曲線的曲率和撓率計(jì)算成為近些年的熱點(diǎn)研究問題。此外，在實(shí)際應(yīng)用中，受人工、儀器、環(huán)境等諸多因素影響，得到的離散曲線的數(shù)據(jù)(一個(gè)序列的三維坐標(biāo)點(diǎn))常常存在一定的偏差或者噪聲，因此如何使算法更加準(zhǔn)確和魯棒，成為三維空間離散曲線的曲率和撓率計(jì)算的難點(diǎn)。

近些年提出的關(guān)于三維空間離散曲線的曲率和撓率計(jì)算方法大致可以分為2類。1) 基于曲線擬合的方法，這類算法主要是對(duì)給定點(diǎn)及其近鄰點(diǎn)擬合一段光滑曲線，然后根據(jù)這段解析曲線的微分幾何公式計(jì)算給定點(diǎn)的曲率和撓率。其本質(zhì)就是用擬合曲線的曲率和撓率來估算離散曲線的曲率和撓率。這類算法的典型算法有多項(xiàng)式擬合法(polynomial curve fitting)[16]和B樣條曲線擬合法(B-spline curve fitting)[17-19]。2) 基于離散結(jié)構(gòu)的方法，這類算法無需對(duì)離散曲線進(jìn)行光滑曲線的擬合，直接利用離散點(diǎn)之間的距離和內(nèi)接多邊形等離散結(jié)構(gòu)，運(yùn)用差分形式或者離散逼近等策略估算出離散曲線上各個(gè)點(diǎn)的曲率和撓率?；陔x散結(jié)構(gòu)的方法中比較典型的有投影法(projection)[20]、單側(cè)差分法(one-sided difference)[21]和離散幾何法(discrete geometry)[22]。本文根據(jù)導(dǎo)數(shù)定義以及單側(cè)差分法的設(shè)計(jì)思路，融合差商平滑策略，提出了微中心差分(microcentral difference)算法。

1 離散曲線的曲率撓率算法

經(jīng)典微分幾何中，設(shè)曲線 C的參數(shù)方程為r=r(t)=(x(t),y(t),z(t))(α≤t≤β)，若一階導(dǎo)函數(shù)r′(t) 在區(qū)間 [α,β] 上 恒不為零，即參數(shù)化曲線 r(t)是正則的，則曲率 κ(t)的計(jì)算公式[23]為

撓率 τ(t)的計(jì)算公式[16]為

式(1)、式(2)中的符號(hào)“×”表示叉積；式(2)中的分子表示3個(gè)微分向量的混合積。r(t)=(x(t)，y(t)，z(t))中要求3個(gè)函數(shù)x(t)、y(t)、z(t)都是連續(xù)函數(shù)，且都3階可導(dǎo)。

但是對(duì)于沒有方程的三維離散曲線的曲率撓率計(jì)算，研究人員提出的各種方法中有代表性的方法有：多項(xiàng)式曲線擬合法、B樣條曲線擬合法、投影法、單側(cè)差分法、離散幾何法5種。本文在分析這5種算法的基礎(chǔ)上，提出了微中心差分法。

1.1 多項(xiàng)式曲線擬合法

2005年，Cazals等[16]提出了一種基于二次多項(xiàng)式曲線擬合(polynomial curve fitting，PCF)的方法來計(jì)算離散平面曲線在給定點(diǎn) Pi的曲率、撓率。該方法對(duì)點(diǎn) Pi及 其鄰域 {Pk|i-n≤k≤i+n}(n為點(diǎn) Pi一 側(cè)選取的鄰近點(diǎn)個(gè)數(shù))，在以 Pi為原點(diǎn)的局部坐標(biāo)系中擬合一個(gè)二次多項(xiàng)式曲線。二次多項(xiàng)式曲線的數(shù)學(xué)模型為

通過最小二乘擬合，可求解該二次多項(xiàng)式曲線及其系數(shù) Ai(i=0,1,2)， 利用多項(xiàng)式系數(shù) Ai可計(jì)算出點(diǎn) Pi的 曲率 κ(Pi)。類似的可推廣擬合三次或三次以上的多項(xiàng)式曲線，從而利用式(1)、式(2)計(jì)算出三維空間離散曲線的曲率、撓率。

1.2 B樣條曲線擬合法

2017年，趙夫群等[17]用4次B樣條曲線擬合(B-spline curve fitting，BSCF)方法將空間離散曲線擬合成光滑的空間曲線，從而以此來計(jì)算離散曲線的曲率和撓率。B樣條曲線的定義為：給定n=m+k+1個(gè) 頂點(diǎn)，則可以定義 m+1段 k次參數(shù)曲線第i段B樣條曲線函數(shù)可表示為

式中： i=0,1,···,n ； { Pi+0,Pi+1,···,Pi+k} 為定義第i段曲線特征多邊形的 k+1個(gè) 頂點(diǎn)； Nl,k(t)為B樣條基底函數(shù)，表示為

對(duì)于離散曲線上的每一個(gè)離散點(diǎn) Pi，選取點(diǎn)Pi及其鄰域 {Pk|i-n≤k≤i+n}(n 為點(diǎn) Pi一側(cè)選取的鄰近點(diǎn)個(gè)數(shù))內(nèi)的點(diǎn)進(jìn)行樣條曲線擬合。類似的方法有：喻德生等[18]采用的3次均勻B樣條進(jìn)行曲線擬合，Kehtarnavaz等[19]采用的5次B樣條進(jìn)行曲線擬合，從而利用經(jīng)典微分幾何的式(1)、式(2)求解方法對(duì)離散曲線進(jìn)行曲率、撓率的計(jì)算。

1.3 投影法

Mardia在1999年提出一種基于投影(projection)的算法[20]，對(duì)于三維空間中的離散曲線，將點(diǎn)投影到單位切線方向上。對(duì)于由 n個(gè)點(diǎn)構(gòu)成的離散曲線，其點(diǎn) Pi處 的單位切向量定義為

容易計(jì)算點(diǎn) Pi處的單位法向量 Ni，利用 Ti和Ni做叉積得到副法向量 Bi， 從而點(diǎn) Pi處的局部投影平面離散曲線定義為

然后以點(diǎn) Yi為中心，擬合平面二次曲線ax2+bx ，從而三維空間離散曲線點(diǎn) Pi的 曲率 κ(Pi)、撓率 τ(Pi)分別為

1.4 單側(cè)差分法

2016年Blankenburg等[21]通過利用一般參數(shù)方程的曲線導(dǎo)數(shù)的有限差商逼近，直接從離散點(diǎn)中采用單側(cè)差分(one-sided difference，OSD)逼近計(jì)算出空間離散曲線在給定點(diǎn) Pi的一階差商、二階差商和三階差商，然后代入曲率和撓率的解析計(jì)算公式進(jìn)行求解。在該方法中，定義差分為

則一階導(dǎo)數(shù)、二階導(dǎo)數(shù)分別為

同理，三階導(dǎo)數(shù)為

從而根據(jù)曲線曲率、撓率的微分定義，將上述各階導(dǎo)直接代入曲率、撓率的計(jì)算式(1)、式(2)即可求得各點(diǎn)處的曲率、撓率。

1.5 離散幾何法

An等[22]在2011年提出一種稱為離散幾何法(discrete geometry)的方法，該方法基于對(duì)稱鄰域的最小二乘的導(dǎo)數(shù)計(jì)算過程，可以歸結(jié)為：利用離散導(dǎo)數(shù)的定義，采用約束最優(yōu)化問題的求解方法來計(jì)算指定點(diǎn)在其對(duì)稱鄰域內(nèi)的離散導(dǎo)數(shù)。

對(duì)于平面內(nèi)的離散曲線，點(diǎn) Pi=(xi,yi)及其鄰域{Pk|i-n1≤k≤i+n1}(n1為點(diǎn) Pi一側(cè)選取的鄰近點(diǎn)個(gè)數(shù))，則點(diǎn) Pi處的離散導(dǎo)數(shù)為

上述離散導(dǎo)數(shù)也稱為點(diǎn) Pi的一階離散導(dǎo)數(shù)。類似的在點(diǎn) Pi處的二階離散導(dǎo)數(shù)為

同理，在點(diǎn) Pi處 的 n階 離散導(dǎo)數(shù)為

上述各階離散導(dǎo)數(shù)是對(duì)連續(xù)可微函數(shù)在點(diǎn)Pi處的各階導(dǎo)數(shù)的逼近。對(duì)于三維空間中的離散曲線，為了計(jì)算其離散導(dǎo)數(shù)，首先要對(duì)離散曲線進(jìn)行參數(shù)化?？刹捎美鄯e弦長(zhǎng)參數(shù)化方法[24]來確定離散點(diǎn)的參數(shù)值。在該方法中，離散曲線 ψd在各點(diǎn)Pi=(xi,yi,zi)(1 ≤i≤n) 處的累積弦長(zhǎng)參數(shù) si可確定為

并且 s1=0 ， 從而離散曲線 ψd為

式中： si∈{s1,s2,···,sn}； 坐標(biāo)函數(shù) xi=x(si)、 yi=y(si)和zi=z(si)就 是弦長(zhǎng) si的離散函數(shù)。因此，平面上的離散導(dǎo)數(shù)的計(jì)算方法可以轉(zhuǎn)化為用于計(jì)算三維空間中的離散曲線上各點(diǎn)的導(dǎo)數(shù)。從而，采用上述方法，對(duì)三維空間中的離散曲線 ψd上的點(diǎn) Pi，其鄰域?yàn)閧Pk|i-n1≤k≤i+n1}(n1為 點(diǎn) Pi一側(cè)選取的鄰近點(diǎn)個(gè)數(shù))，則點(diǎn) Pi處 的離散導(dǎo)數(shù) (x′(si),y′(si),z′(si))分別為

上述離散導(dǎo)數(shù)也稱為點(diǎn) Pi的一階離散導(dǎo)數(shù)。類似的可計(jì)算點(diǎn) Pi的二階、三階離散導(dǎo)數(shù)。從而利用式(1)、式(2)即可求得三維空間離散曲線的曲率、撓率。

1.6 微中心差分法

單側(cè)差分法[22]在進(jìn)行離散曲線的曲率撓率計(jì)算時(shí)，會(huì)對(duì)給定的曲率撓率產(chǎn)生一定的偏移或者縮放，為了克服由單側(cè)差分法引起的這種計(jì)算偏差，結(jié)合文獻(xiàn)[22]中的思想和導(dǎo)數(shù)定義，把其單側(cè)差分法修改為對(duì)目標(biāo)離散點(diǎn)的雙側(cè)差商取平均的方法來得到目標(biāo)離散點(diǎn)的各階導(dǎo)數(shù)，并參考文獻(xiàn)[25]中的性能度量算法思想和差分方式，將改進(jìn)的方法命名為微中心差分法(microcentral difference algorithm)。

首先對(duì)微中心差分法名稱中的“微”和“中心”作簡(jiǎn)單的說明?！拔ⅰ笔侵冈诰植糠秶闹痣A差商計(jì)算。比如，先用改進(jìn)的離散導(dǎo)數(shù)定義計(jì)算一階差商，如式(3)；再利用一階差商值計(jì)算離散曲線的二階差商，如式(4)；最后利用二階差商值計(jì)算三階差商，如式(5)?！爸行摹敝冈谀繕?biāo)點(diǎn)鄰域內(nèi)，進(jìn)行差分計(jì)算時(shí)，以離散曲線的目標(biāo)點(diǎn)為中心進(jìn)行差分計(jì)算。另外，方法中涉及求目標(biāo)點(diǎn)的各階差商，我們指定，每階差商目標(biāo)點(diǎn)單側(cè)需要上一階差商值的個(gè)數(shù)，記為 ki(i=1,2,3)，稱為鄰域半徑。根據(jù)上述說明，微中心差分算法的計(jì)算方法具體如下。定義其一階前向差商為

從而其一階差商為

式中 k1為鄰域半徑。

其二階前向差商定義為

二階后向差商定義為

從而其二階差商定義為

式中 k2為鄰域半徑。

同理，其三階前向差商定義為

三階后向差商定義為

從而其三階差商為

式中 k3為鄰域半徑。

因此對(duì)于微中心差分法，為了得到離散曲線上給定的目標(biāo)離散點(diǎn)的曲率值，根據(jù)本文的計(jì)算方法，則需要知道目標(biāo)離散點(diǎn)鄰域內(nèi)的離散點(diǎn)數(shù)N1=2(k1+k2)+1， 即目標(biāo)離散點(diǎn)兩側(cè)各需 (k1+k2)個(gè)離散點(diǎn)；而對(duì)于目標(biāo)離散點(diǎn)的撓率計(jì)算，則至少需要知道目標(biāo)離散點(diǎn)鄰域內(nèi)的離散點(diǎn)數(shù)N2=2(k1+k2+k3)+1， 其中目標(biāo)離散點(diǎn)兩側(cè)各(k1+k2+k3)個(gè)離散點(diǎn)，k1、k2、k3均為大于等于1的整數(shù)。根據(jù)式(3)～(5)得到其各階差商后，則可根據(jù)曲線曲率、撓率的微分定義，將上述各階導(dǎo)數(shù)式(3)～(5)代入曲率、撓率的計(jì)算式(1)、式(2)，即可求得離散曲線上各點(diǎn)處的曲率、撓率。這里提出的微中心差分計(jì)算方法是對(duì)單側(cè)差分方法的一個(gè)依據(jù)導(dǎo)數(shù)定義進(jìn)行的一個(gè)擴(kuò)展，理論上是符合連續(xù)函數(shù)條件下可導(dǎo)的要求。

總體來看，上述6種方法分為2種思路：1)一對(duì)近鄰點(diǎn)進(jìn)行曲線擬合，它們之間的區(qū)別在于近鄰點(diǎn)數(shù)的確定以及擬合曲線模型的假設(shè)；2)對(duì)近鄰點(diǎn)進(jìn)行離散求導(dǎo)，各方法之間的區(qū)別在于取點(diǎn)的策略以及如何利用近鄰點(diǎn)的信息得到各階導(dǎo)函數(shù)的可靠逼近，從而得到不同的曲率和撓率的計(jì)算方法。

2 實(shí)驗(yàn)比較

本節(jié)通過實(shí)驗(yàn)對(duì)上述6種算法計(jì)算的曲率和撓率的計(jì)算精度進(jìn)行評(píng)估。為了便于比較，本文中實(shí)驗(yàn)方法的k1、k2、k3均取值為1。算法評(píng)估指標(biāo)參考文獻(xiàn)[26]的評(píng)價(jià)指標(biāo)，選取最大相對(duì)誤差(max relative error，MRE)和均方根誤差(root mean square error，RMSE)，即對(duì)于獲得的離散曲線，在每一個(gè)離散點(diǎn) Pi處的偏差定義為其中為實(shí)驗(yàn)方法的計(jì)算值， φi為對(duì)應(yīng)的準(zhǔn)確值。因此，對(duì)于給定的離散曲線，其最大相對(duì)誤差為EMRE=, 均方根誤差為 E最大相對(duì)誤差體現(xiàn)算法的穩(wěn)定性，均方根誤差反映算法的整體有效性。實(shí)驗(yàn)平臺(tái)為個(gè)人筆記本電腦，計(jì)算機(jī)配置是Inter(R) Core(TM) i5-7300HQ CPU @2.50 GHz處理器和8 GB內(nèi)存，程序運(yùn)行環(huán)境是MATLAB R2014b。

2.1 實(shí)驗(yàn)數(shù)據(jù)

實(shí)驗(yàn)用如圖1所示的6條連續(xù)可微參數(shù)曲線進(jìn)行分析。由于有曲線方程，可以用經(jīng)典微分幾何的方法(式(1)、式(2))計(jì)算出曲線上任意一個(gè)點(diǎn)的準(zhǔn)確曲率值和撓率值，有了準(zhǔn)確值便于進(jìn)行誤差分析。實(shí)驗(yàn)分析用到的6條曲線參數(shù)方程如下：

1)三次撓曲線：

2) 多變式內(nèi)接彈簧線：

3) 循環(huán)曲線：

4) Viviani曲線：

5) Clelia曲線：

6) 環(huán)面螺線：

需要說明的是，解析曲線方程用于生成離散曲線，并用于計(jì)算離散曲線上每個(gè)采樣點(diǎn)的準(zhǔn)確的曲率值和撓率值。實(shí)驗(yàn)中對(duì)每條空間參數(shù)曲線通過均勻采樣獲得相應(yīng)的離散曲線，采樣點(diǎn)數(shù)n在定義域內(nèi)從50個(gè)點(diǎn)變化到500個(gè)點(diǎn)，用于測(cè)試算法的有效性。第1章介紹的5種常用的離散曲線曲率和撓率計(jì)算方法以及本文提出的算法都采用離散曲線上的點(diǎn)的坐標(biāo)進(jìn)行計(jì)算，都用不到曲線方程。

圖 1 用于實(shí)驗(yàn)的6條空間曲線(線顏色與撓率值對(duì)應(yīng))Fig. 1 Six 3D curves used in the experiment (The color is corresponding to torsion value)

2.2 誤差分析及時(shí)間效率比較

通過上述三維空間6條離散曲線來定量地評(píng)價(jià)本文提出的微中心差方法(MCD)，并將其與第1章介紹的B樣條曲線擬合法(BSCF)、多項(xiàng)式曲線擬合法(PCF)、投影法(Proj)、單側(cè)差分法(OSD)、離散幾何法(DG)這5種算法的曲率和撓率的計(jì)算結(jié)果進(jìn)行比較。圖2～7分別給出了這6條離散曲線隨著采樣密度增加，曲率和撓率的最大相對(duì)誤差(MRE)與均方根誤差(RMSE)的變化。從圖2～7的誤差比較圖可以看出，隨著采樣密度 n的不斷增加，這6種方法計(jì)算的曲率和撓率的最大相對(duì)誤差和均方根誤差都逐漸地減小，并且采樣密度越高，最大相對(duì)誤差值和均方根誤差值越接近于0，這表明這6種算法都能有效地計(jì)算出三維空間離散曲線的曲率和撓率。

當(dāng)然，這6種方法在收斂速度和穩(wěn)定性方面存在一些差異。從圖2～7的曲率MRE對(duì)比子圖可以看出：投影法和單側(cè)差分法的最大誤差最大，且隨著采樣密度增加，誤差一直比較大；離散幾何法與本文提出的微中心差分法的誤差較小，特別是當(dāng)采樣較稀疏時(shí)，曲率估算的最大相對(duì)誤差最小。從圖2～7的這6個(gè)圖中的每個(gè)子圖(c)圖的結(jié)果，也說明從均方根指標(biāo)評(píng)價(jià)，也能得出同樣結(jié)論。需要指出的是，從均方根誤差指標(biāo)看，本文提出的算法在所有離散曲線的計(jì)算結(jié)果中幾乎都能取得最佳效果。

圖2～7分別給出了6條離散曲線隨著采樣密度增加的撓率最大相對(duì)誤差(MRE)。從圖中可以看出，隨著采樣密度 n的不斷增加，本文方法的撓率最大相對(duì)誤差逐漸地減小，并且采樣密度越大，最大相對(duì)誤差值和均方根誤差值就越接近于0，這表明撓率計(jì)算越準(zhǔn)確，也表明本文方法在計(jì)算三維空間離散曲線的撓率時(shí)收斂性和穩(wěn)定性較好。而單側(cè)差分法、投影法和B樣條曲線擬合法的最大相對(duì)誤差有時(shí)會(huì)有較大波動(dòng)。

圖2～7中(c)圖呈現(xiàn)的撓率計(jì)算的均方根誤差變化結(jié)果表明，隨著采樣密度增大，本文方法與離散幾何方法和多項(xiàng)式擬合方法的均方根誤差能單調(diào)遞減收斂于0，而單側(cè)差分法、投影法的誤差收斂速度較差。

圖 2 三次撓曲線誤差對(duì)比實(shí)驗(yàn)Fig. 2 Error comparison experiment of cubic deflection curve

圖 3 多變式內(nèi)接彈簧線誤差對(duì)比實(shí)驗(yàn)Fig. 3 Error comparison experiment of changeable inner spring curve

圖 4 循環(huán)曲線誤差對(duì)比實(shí)驗(yàn)Fig. 4 Error comparison experiment of circular curve

圖 5 Viviani曲線誤差對(duì)比實(shí)驗(yàn)Fig. 5 Error comparison experiment of Viviani curve

圖 6 Clelia曲線誤差對(duì)比實(shí)驗(yàn)Fig. 6 Error comparison experiment of Clelia curve

圖 7 環(huán)面螺線誤差對(duì)比實(shí)驗(yàn)Fig. 7 Error comparison experiment of toroidal spiral curve

圖2～7反映了隨著點(diǎn)的密度增加，曲率和撓率的誤差的收斂性和穩(wěn)定性。下面針對(duì)固定的點(diǎn)密度，比較一下各個(gè)方法的精度。實(shí)驗(yàn)中我們根據(jù)上述的密度實(shí)驗(yàn)，取采樣點(diǎn)數(shù)量適中誤差比較穩(wěn)定的點(diǎn)密度，即n=200的時(shí)候進(jìn)行6條曲線的曲率和撓率估計(jì)精度的比較，結(jié)果如表1所示。

表 1 6條曲線的平均誤差比較(n=200)Table 1 Comparison of average error (n = 200)

從表1可以看出，本文提出的MCD方法計(jì)算的曲率和撓率平均誤差分別為0.000 4和0.000 3，比另5種算法中效果最佳的DG算法的曲率和撓率平均誤差(分別為0.001 4和0.000 8)，分別降低了69.74%和55.36%?？梢姡傮w來看，本文提出的算法無論是曲率估計(jì)還是撓率估計(jì)，精度都是最高的。

關(guān)于運(yùn)行效率，從算法原理可以看出，所有算法時(shí)間復(fù)雜度都是線性級(jí)的，即O(n)，因此給出平均實(shí)驗(yàn)時(shí)間便可以說明計(jì)算的時(shí)間效率。用6種算法分別計(jì)算6條曲線的曲率和撓率，對(duì)每條離散曲線的不同采樣密度進(jìn)行50次實(shí)驗(yàn)，總的平均運(yùn)行時(shí)間如表2所示。從表2可以看出，離散結(jié)構(gòu)法(Proj、OSD、DG、MCD)相比曲線擬合法(BSCF、PCF)，計(jì)算效率大幅提高；在離散結(jié)構(gòu)法中，單側(cè)差分算法(OSD)平均運(yùn)行速度最快，本文提出的微中心差分法處于中間位置。

表 2 6種算法的平均運(yùn)行時(shí)間比較Table 2 Comparison of average running time of thesix algorithms s

2.3 噪聲數(shù)據(jù)實(shí)驗(yàn)

為了測(cè)試每個(gè)算法的抗噪聲魯棒性，在采樣點(diǎn)數(shù) n為200時(shí)，設(shè)置噪聲水平為0.000 1～0.1，噪聲為均值為零、標(biāo)準(zhǔn)差為離散曲線上鄰近的兩離散點(diǎn)之間的平均距離與噪聲水平的乘積的高斯噪聲。實(shí)驗(yàn)用Clelia曲線生成了不同噪聲水平的離散曲線。

6種算法計(jì)算的曲率最大相對(duì)誤差和曲率均方根誤差分別如圖8(a)、(b)所示，從圖中可以看出：噪聲水平越高，算法的誤差變化越大；各算法在噪聲水平不大于0.01時(shí)，均能較準(zhǔn)確地計(jì)算離散曲線的曲率值；噪聲水平大于0.01時(shí)，其誤差明顯變大，尤其是單側(cè)差分法和投影法，而本文的方法誤差變化浮動(dòng)相對(duì)較小。

6種算法計(jì)算的撓率最大相對(duì)誤差和均方根誤差分別如圖8(c)、(d)所示，從圖中可以看出：含噪聲的離散曲線在計(jì)算撓率時(shí)變得更加敏感，當(dāng)噪聲水平大于0.001時(shí)，6種算法得到的撓率值與真實(shí)值就已有了明顯的偏差，并且隨著噪聲水平變大，撓率的最大相對(duì)誤差和均方根誤差相比曲率也會(huì)更快地增大。與其他5種算法相比，本文的方法誤差變化浮動(dòng)相對(duì)較小。

圖 8 Clelia曲線的不同噪聲水平下的MSE和RMSE(n =200)Fig. 8 Max relative error and root mean square error of curvatures of Clelia curve under different noise levels (n =200)

上述實(shí)驗(yàn)表明，噪聲對(duì)6種算法的計(jì)算精度都有影響，尤其是當(dāng)噪聲水平較大時(shí)，撓率變化表現(xiàn)得更加明顯。另外，在同樣的噪聲水平下，撓率的最大相對(duì)誤差和均方根誤差波動(dòng)幅度相比曲率的最大相對(duì)誤差和均方根誤差波動(dòng)幅度更大，這也說明離散曲線的撓率計(jì)算比曲率計(jì)算對(duì)噪聲更加敏感。

3 結(jié)束語

針對(duì)三維空間離散曲線的曲率和撓率的計(jì)算方法設(shè)計(jì)問題，本文主要工作包括如下2個(gè)方面：

1)根據(jù)2016年Blankenburg等提出的單側(cè)差分方法和導(dǎo)數(shù)定義，結(jié)合差商平滑策略，設(shè)計(jì)了微中心差分法，該方法能夠有效計(jì)算離散曲線的曲率和撓率。從整個(gè)計(jì)算過程可以發(fā)現(xiàn)，本文的方法不依賴于任何曲線模型，從離散曲線中直接計(jì)算所得，從而對(duì)各類離散曲線均具有較好的普適性。

2)概述了多項(xiàng)式曲線擬合法、B樣條曲線擬合法、投影法、單側(cè)差分法和離散幾何法5種典型算法；詳細(xì)闡述了本文提出的微中心差分法的基本思想和計(jì)算公式。采用6條離散曲線數(shù)據(jù)進(jìn)行了不同算法的對(duì)比實(shí)驗(yàn)，實(shí)驗(yàn)表明，本文提出的微中心差分法的平均精度最高。

從實(shí)驗(yàn)結(jié)果可以看出，目前的這些算法的噪聲魯棒性并不是很好，特別是在有奇異點(diǎn)的時(shí)候，未來的研究可以進(jìn)一步研究魯棒性更佳的算法。此外，未來也可以研究利用曲率和撓率進(jìn)行特征線的分析、提取和類型識(shí)別等。