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

慕生鵬，李紅軍，李世林

（北京林業大學 理學院，北京 100083）

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

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

近些年提出的關于三維空間離散曲線的曲率和撓率計算方法大致可以分為2類。1) 基于曲線擬合的方法，這類算法主要是對給定點及其近鄰點擬合一段光滑曲線，然后根據這段解析曲線的微分幾何公式計算給定點的曲率和撓率。其本質就是用擬合曲線的曲率和撓率來估算離散曲線的曲率和撓率。這類算法的典型算法有多項式擬合法(polynomial curve fitting)[16]和B樣條曲線擬合法(B-spline curve fitting)[17-19]。2) 基于離散結構的方法，這類算法無需對離散曲線進行光滑曲線的擬合，直接利用離散點之間的距離和內接多邊形等離散結構，運用差分形式或者離散逼近等策略估算出離散曲線上各個點的曲率和撓率?；陔x散結構的方法中比較典型的有投影法(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)的計算公式[16]為

式(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 多項式曲線擬合法

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

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

1.2 B樣條曲線擬合法

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

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

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

1.3 投影法

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

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

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

1.4 單側差分法

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

則一階導數、二階導數分別為

同理，三階導數為

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

1.5 離散幾何法

An等[22]在2011年提出一種稱為離散幾何法(discrete geometry)的方法，該方法基于對稱鄰域的最小二乘的導數計算過程，可以歸結為：利用離散導數的定義，采用約束最優化問題的求解方法來計算指定點在其對稱鄰域內的離散導數。

對于平面內的離散曲線，點 Pi=(xi,yi)及其鄰域{Pk|i-n1≤k≤i+n1}(n1為點 Pi一側選取的鄰近點個數)，則點 Pi處的離散導數為

上述離散導數也稱為點 Pi的一階離散導數。類似的在點 Pi處的二階離散導數為

同理，在點 Pi處 的 n階 離散導數為

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

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

式中： 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))分別為

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

1.6 微中心差分法

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

首先對微中心差分法名稱中的“微”和“中心”作簡單的說明。“微”是指在局部范圍的逐階差商計算。比如，先用改進的離散導數定義計算一階差商，如式(3)；再利用一階差商值計算離散曲線的二階差商，如式(4)；最后利用二階差商值計算三階差商，如式(5)?！爸行摹敝冈谀繕它c鄰域內，進行差分計算時，以離散曲線的目標點為中心進行差分計算。另外，方法中涉及求目標點的各階差商，我們指定，每階差商目標點單側需要上一階差商值的個數，記為 ki(i=1,2,3)，稱為鄰域半徑。根據上述說明，微中心差分算法的計算方法具體如下。定義其一階前向差商為

從而其一階差商為

式中 k1為鄰域半徑。

其二階前向差商定義為

二階后向差商定義為

從而其二階差商定義為

式中 k2為鄰域半徑。

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

三階后向差商定義為

從而其三階差商為

式中 k3為鄰域半徑。

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

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

2 實驗比較

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

2.1 實驗數據

實驗用如圖1所示的6條連續可微參數曲線進行分析。由于有曲線方程，可以用經典微分幾何的方法(式(1)、式(2))計算出曲線上任意一個點的準確曲率值和撓率值，有了準確值便于進行誤差分析。實驗分析用到的6條曲線參數方程如下：

1)三次撓曲線：

2) 多變式內接彈簧線：

3) 循環曲線：

4) Viviani曲線：

5) Clelia曲線：

6) 環面螺線：

需要說明的是，解析曲線方程用于生成離散曲線，并用于計算離散曲線上每個采樣點的準確的曲率值和撓率值。實驗中對每條空間參數曲線通過均勻采樣獲得相應的離散曲線，采樣點數n在定義域內從50個點變化到500個點，用于測試算法的有效性。第1章介紹的5種常用的離散曲線曲率和撓率計算方法以及本文提出的算法都采用離散曲線上的點的坐標進行計算，都用不到曲線方程。

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

2.2 誤差分析及時間效率比較

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

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

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

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

圖 2 三次撓曲線誤差對比實驗Fig. 2 Error comparison experiment of cubic deflection curve

圖 3 多變式內接彈簧線誤差對比實驗Fig. 3 Error comparison experiment of changeable inner spring curve

圖 4 循環曲線誤差對比實驗Fig. 4 Error comparison experiment of circular curve

圖 5 Viviani曲線誤差對比實驗Fig. 5 Error comparison experiment of Viviani curve

圖 6 Clelia曲線誤差對比實驗Fig. 6 Error comparison experiment of Clelia curve

圖 7 環面螺線誤差對比實驗Fig. 7 Error comparison experiment of toroidal spiral curve

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

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

從表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 6種算法的平均運行時間比較Table 2 Comparison of average running time of thesix algorithms s

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種算法相比，本文的方法誤差變化浮動相對較小。

圖 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)

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

3 結束語

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

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

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

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