吳易澤，張 旭

(上海工程技術大學 機械與汽車工程學院，上海 201620)

0 引言

隨著科技水平的不斷進步，作為計算機輔助幾何設計(Computer Aided Geometric Design, CAGD)中的一個重要分支，曲線曲面造型技術一直發(fā)揮著極其重要的作用，特別是在航空航天、船舶、汽車等對制造精度和外觀有著較高要求的制造業(yè)領域。此外，曲線曲面造型技術對產(chǎn)品質量、物理性能等都將會有重要的影響。因此，曲線曲面造型技術受到工業(yè)界和學術界的廣泛關注，并且一直是CAGD中的研究熱點。

曲線的光順算法種類很多?；诠忭樞Ч煞譃檎w光順和局部光順[1]；基于光順原理可分為節(jié)點刪除和節(jié)點插入；還有基于能量法的光順和基于小波技術的光順等。Farin等[2]提出一種B樣條曲線光順算法。Sapidis等[3]提出一種B樣條曲線的自動光順算法。龍小平[4]提出一種B樣條曲線的局部能量算法。高斌等[5]提出一種基于曲率變化約束的平面B樣條曲線光順方法，通過曲率單調性約束條件，考慮曲線段的離散能量及原始曲線的誤差控制，從而實現(xiàn)了曲線光順。王士瑋等[6]通過將曲線光順問題建模成基于稀疏模型的優(yōu)化問題，從而提出一種新的曲線光順算法。

近年來，基于小波分析的曲線曲面光順方法受到了國內外專家學者的廣泛關注。Amati[7]提出了一種基于小波分析的曲線光順算法。Wang等[8]基于小波分析，提出了非均勻有理B樣條(Non-Uniform Rational B-Splines, NURBS)曲線曲面的光順方法。趙罡等[9]針對一般的非均勻B樣條曲線，提出了一種基于非均勻B樣條小波的曲線光順算法。

事實上，由于數(shù)字曲線與離散信號之間的相似性，可通過信號分析的方法來處理數(shù)字曲線的光順性。顯然，數(shù)字曲線是一種非平穩(wěn)信號。處理非平穩(wěn)信號的方法包括傅里葉變換、小波分析、Gabor展開等。上述方法雖然都能在不同程度上處理非平穩(wěn)信號，但都有其局限性，無法從根本上解決傅里葉變換在信號分析中的問題。為解決該問題，Huang等[10]提出一種算法——經(jīng)驗模態(tài)分解(Empirical Mode Decomposition, EMD)，成功地解決了這個問題。為有效地區(qū)分出復雜的非平穩(wěn)信號中的信號和噪聲，該算法自適應地將其從高頻到低頻分解成有限個本征模態(tài)函數(shù)(Intrinsic Mode Function, IMF)及一個殘余項。

基于EMD方法，秦緒佳等[11]提出一種通過構造近似均值曲線并以此將原始曲線進行一維參數(shù)化的光順方法。譚小俊等[12]提出一種雙變量EMD分析的平面離散曲線光順方法，取得了不錯的光順效果?？墒沁@些方法僅去除了一次EMD分解后所得的IMF分量或多次EMD分解后所得的多個IMF分量，去除的干擾的IMF分量中可能含有有效的特征信息，保留下的IMF分量也可能含有噪聲信息，且分解的次數(shù)需要根據(jù)特定的曲線設定。此外，EMD分解還存在某些不足，如端點效應、模態(tài)混疊。

針對EMD方法存在的不足之處，Wu等[13]提出一種算法，即集合經(jīng)驗模態(tài)分解(Ensemble Empirical Mode Decomposition, EEMD)，成功地彌補了該算法的不足，也即EEMD方法能完美避免尺度混合，并能使最終分解的IMF分量保持物理上的唯一性。基于此，本文考慮將EEMD方法應用于曲線光順，若對EEMD分解后的若干個分量逐一進行降噪處理，在增加工作量的同時，還會降低光順效果，并且會增大誤差。因而本文提出采用游程檢測法重構EEMD分解后的分量，將其重構為高頻分量和低頻分量。此時噪聲主要集中在高頻分量，單獨去除高頻分量在產(chǎn)生較大誤差的同時，會去除曲線的有效特征。由于奇異譜分析(Singular Spectrum Analysis, SSA)方法在非線性、非平穩(wěn)信號的分析方面，屬于一種有效的分析方法，它可以顯著地降低信號中的噪聲，有效地提高信噪比[14-15]，并能保留噪聲信號頻帶相混疊的有用信號，突出信號的信息特征。因此，本文采用奇異譜分析對高頻分量進行降噪處理，在一定程度上保留了曲線的有效信息，避免了曲線光順后的擬合誤差過大。試驗表明，本文提出的EEMD分解、游程檢測法重構和SSA降噪三者相結合的方法對曲線具有很好的光順性。

1 算法原理

1.1 EMD與EEMD算法

EMD是由Huang等[10]提出的處理非平穩(wěn)和非線性信號的信號處理方法，該方法可將復雜的原始數(shù)據(jù)信號按不同尺度從高頻到低頻自適應地分解成若干個IMF分量和一個剩余分量Res。

EMD分解是一個“篩選”信號的過程，其詳細分解步驟參考文獻[16-19]。分解后的每個IMF分量需滿足兩個條件：

(1)對于任意分解得到的IMF分量中，過零點、極值點的個數(shù)之差小于等于1；

(2)由局部極大值點、局部極小值點分別確定的上、下包絡線均值為0。

EEMD作為一種信號處理算法，其在非線性信號、非平穩(wěn)信號的分析處理上具有舉足輕重的地位，能夠利用噪聲特性有效地克服由于EMD方法容易產(chǎn)生模態(tài)混疊現(xiàn)象的問題[13]。該方法本質上是對添加了高斯白噪聲的信號進行多次EMD分解，同時利用高斯白噪聲的隨機性，可以經(jīng)過反復試驗降低相應的IMF分量的模式混疊程度。該方法的具體分解步驟如下：

(1)將高斯白噪聲n(t)添加到目標信號y(t)，得到一個新的信號x(t)：

x(t)=y(t)+n(t)。

(1)

式中，高斯白噪聲n(t)服從下式：

(2)

式中，ε為噪音強度參數(shù)；N為總體個數(shù)；εn為信號的標準差。

(2)按照EMD分解方法對信號x(t)進行分解，得到1組IMF分量ci(t)和1個余量r(t)：

(3)

式中n為IMF分量個數(shù)。

(3)重復步驟(1)和步驟(2)，每次添加不同的白噪聲nj(t)，得到M組IMF分量和剩余分量。

(4)利用高斯白噪聲頻譜的零均值原理，將分解出的M組IMF分量和剩余分量的整體平均值作為最終結果，目標信號對應的IMF分量ci(t)可表示為：

(4)

1.2 游程檢測法

游程檢測亦稱“連續(xù)性檢驗”，是一種利用游程總數(shù)來判斷樣本隨機性的統(tǒng)計檢驗方法[20-21]，其中游程是指樣本序列中連續(xù)出現(xiàn)的變量值的次數(shù)。

(5)

式中，st由一組相互獨立統(tǒng)計的隨機分布的“0-1”序列組成，將每段連續(xù)相同符號(0或1)的序列定義為一個游程，將每一段游程中的數(shù)據(jù)個數(shù)定義為一個游程長度。

在游程檢測法中，可以根據(jù)游程總數(shù)的大小判斷出對應的EEMD分量的波動程度。進而，依據(jù)游程檢測法設定高頻、低頻分量的游程閾值和游程長度閾值。若游程總數(shù)大于閾值，同時游程長度較短，此時序列變化較大；反之，游程總數(shù)小，同時游程長度長，此時序列變化小。在每個分量的游程數(shù)和游程長度經(jīng)過綜合考量之后，EEMD分解后的分量被重構為高頻、低頻兩種分量。

1.3 奇異譜分析算法

1978年，Colebrook提出了奇異譜分析算法，它是一種主成分分析方法[22]。由于其時效性強、計算簡單，非常適用于非線性、非平穩(wěn)信號的處理[23]。

假設有某一一維信號x(i)(i=1,2,…n)，給定嵌入維數(shù)為m(m

(6)

令S為時滯矩陣的m×m維協(xié)方差矩陣，則

(7)

采用奇異譜分析對協(xié)方差矩陣S進行分解，從而得到m個奇異值λi(i=1,2,…,m)。為了將該信號的奇異譜圖構造出來，通過將得到的m個奇異值進行降序排列λ1>λ2>…>λm≥0。其中，奇異值大小表示信號和噪聲在奇異譜圖中能量大小的相對關系。將值較大的奇異值點看成信號點，而將值較小的奇異值點看成噪聲點。λk對應的特征向量Ek稱作經(jīng)驗正交投影函數(shù)，采樣信號x(i)在特征向量Ek上的正交投影系數(shù)就是第k個主分量：

(8)

如果已知各主分量與經(jīng)驗正交函數(shù)，反求原始信號序列的過程如下：

(9)

2 光順過程

在曲線的光順過程中，首先將空間離散數(shù)字曲線的x,y,z三個變量視為3個一維數(shù)字信號；其次對每個變量的數(shù)字信號序列分別進行EEMD分解，將每個變量分解成多個IMF分量和一個Res余量；接著分別對每個變量分解后的所有分量使用游程檢測法，將其重構為高頻、低頻分量；隨后對重構后的高頻分量使用SSA進行降噪；最終將降噪后的高頻分量與低頻分量進行重構，進而得到光順后的曲線。其曲線光順流程圖如圖1所示。

具體步驟如下：

(1)原始數(shù)字曲線信號的EEMD分解。將空間離散數(shù)字曲線的x,y,z三個變量視為3個一維數(shù)字信號，然后對每個變量的數(shù)字信號序列分別進行EEMD分解，將其分解成n個具有不同特征的IMF1～IMFn分量和一個Res余量。

(2)由于在EEMD分解中，得到了數(shù)量相對較多的IMF分量，若對每個分量逐一進行降噪處理，會去除曲線的有效信息，同時會加大光順工作量和增加曲線光順的誤差，這就需要對每個變量的序列分量進行重構。其重構步驟如下：

1)計算游程數(shù)、最大游程長度。計算每個變量所有的IMF分量、Res余量的游程數(shù)以及與之相對應的最大游程長度；

2)確定閾值。根據(jù)每個分量的游程數(shù)和最大游程長度，選擇游程數(shù)變化率最大的游程數(shù)作為閾值；

3)劃分高頻、低頻分量。比較閾值與每個分量的游程數(shù)，得到高頻、低頻分量；

4)根據(jù)判斷結果，對x,y,z三個變量的各個分量分別進行疊加，最終得到重構后的高頻、低頻分量。

(3)奇異譜分析降噪。對進行游程檢測法后得到的高頻分量進行奇異譜分析降噪。在奇異譜分析降噪的過程中，重構信號時階次p的選擇是一個重要環(huán)節(jié)。其降噪步驟如下：

1)使用Cao方法[24]確定最小嵌入維數(shù)m(本文選取m=40)，同時構建時滯矩陣X；

2)構建矩陣X的m維協(xié)方差矩陣S，對其進行奇異值分解，從而得到m個奇異值；

3)確定信號的奇異譜圖，并將奇異值按降序排列，從而確定信號的奇異譜圖；

4)選取主分量作為重構后的高頻分量，其中，主分量的重構主要選取前面幾個奇異值較大的點。

(4)光順結果重構。將降噪后的高頻分量與低頻分量進行重構，最后得到光順后的曲線。

(5)光順結果分析。采用3種不同的曲線光順算法進行光順性結果對比。采用曲率圖作為評價指標。曲率圖中曲率變化平穩(wěn)，說明曲線光順；曲率變化變化大，說明曲線不光順。其中，曲率計算公式為：

k=

(10)

3 試驗與結果分析

3.1 基于EEMD的空間離散曲線上的x,y,z三變量分解

根據(jù)1.1節(jié)所述的EEMD分解步驟，對空間離散曲線上的x,y,z三個變量分別進行EEMD分解。通過查閱文獻[10]，本文所選取的曲線離散數(shù)據(jù)點個數(shù)為100，則原始數(shù)據(jù)點分解為100組IMF分量，故設置M=100，α為顯著性水平，一般設α=0.25，ε為噪音強度參數(shù)，取值范圍為(0～1)，由于數(shù)據(jù)點個數(shù)較少，取值為1/4，故設ε=0.25。

EEMD分解步驟完成后，最終x變量的時間序列被分解成5個IMF分量和1個剩余分量Res；y變量的時間序列被分解成5個IMF分量和1個剩余分量Res；z變量的時間序列被分解成5個IMF分量和1個剩余分量Res。其原始曲線及其分解結果如圖2所示。

觀察圖2可以發(fā)現(xiàn)，空間離散數(shù)字曲線的x,y,z三個變量序列的非線性、非平穩(wěn)性特性相當明顯。其中，IMF1～IMF5分量分別反映了離散數(shù)字曲線的數(shù)據(jù)點序列在不同頻率上的變化信息，余項Res則展示出曲線在光順后的變化。

3.2 游程檢測法重構

由于各個變量經(jīng)過EEMD分解后，得到的分量過多，各個變量上的噪聲在每個IMF分量上分布不均，且在IMF1～IMF5以及Res上噪聲呈現(xiàn)由高到低的趨勢。因此，在對所有分量進行降噪處理之前，需要先對其進行游程檢測法重構，得到高頻分量和低頻分量。

分別對x,y,z三個變量的5個IMF分量和1個Res余量使用游程檢測法，將得到3個變量序列的游程數(shù)、最大游程長度，結果如表1所示。

表1 分量的游程數(shù)與最大游程長度

續(xù)表1

分析表1，綜合考慮游程數(shù)和最大游程長度，設游程數(shù)的閾值為10，設最大游程長度的閾值為50，將游程數(shù)大于10且最大游程長度小于50的分量合并為高頻分量，反之則合并為低頻分量。

最終得到x變量的高頻分量由IMF1～IMF5疊加而成，x變量的低頻分量由Res疊加而成；將y變量的IMF1～IMF4疊加為其高頻分量，IMF5和Res疊加為其低頻分量；z變量的IMF1～IMF4疊加為其高頻分量，IMF5和Res疊加為其低頻分量。最后重構得到x,y,z三個變量的高頻分量和低頻分量，如圖3所示。

3.3 奇異譜分析降噪

由3.2節(jié)所述可知，數(shù)據(jù)噪聲主要集中在高頻分量。因此，對高頻分量進行降噪處理是必不可少的。由于在信號降噪中SSA的效果較優(yōu)，若對含噪數(shù)據(jù)運用SSA降噪，如果信號重構階次過大則會殘留大量噪聲，反之會損失有效數(shù)據(jù)信息。因此，SSA信號重構階次的選擇是個關鍵問題，合適的階次可以很好地保留有效數(shù)據(jù)，反之不合適的階次則會導致數(shù)據(jù)信號失真或者噪聲過多。重構階次的選擇有多種，通過對比，本文最終使用奇異譜構造三角形求余弦值的方法確定重構階次，進而重構數(shù)字變量，具體步驟參考文獻[25]。

根據(jù)余弦定理公式，對其變形得：

式中：a為三角形的第一個邊長；b為三角形的第二個邊長；c為三角形的第三個邊長；C為邊長c所對的角。

眾所周知，在[0,π]上，余弦函數(shù)單調遞減，因此最大的余弦值對應著最小的角。因此，拐點位置即為最大的余弦值所處的位置，即重構階次p值。通過連接首尾奇異值、中間奇異值，構建三角形，如圖4所示。

分析圖4可知，圖中特征值按照降序排列，橫、縱坐標軸分別表示奇異值個數(shù)、與之相對應的奇異值，圖中奇異值按照降序排列。正常情況下，將圖4中奇異值趨于平坦的部分認定為噪聲奇異值，因此，由較大的奇異值組成了變量的主要成分，小的奇異值組成了噪聲。由此可知，x,y,z高頻分量的重構階次如表2所示。

表2 重構階次P值的選取

選取前p個奇異值進行信號重構，計算公式如式(12)所示：

(12)

得到降噪后的高頻分量，如圖5所示。

3.4 光順結果與對比分析

在驗證本文算法有效性的同時，選取2種典型算法與本文算法一起對同一條曲線進行對比試驗。3種優(yōu)化算法的光順結果如圖6所示。

本文所選取的2種典型的光順方法分別為EMD法和曲率法。其中，EMD法是通過EMD分解原始曲線數(shù)據(jù)，得到頻率不同的IMF分量，去除高頻分量，保留低頻分量，最終得到光順后的曲線，該算法是從數(shù)據(jù)實際出發(fā)來對數(shù)據(jù)進行處理；曲率法是通過曲率變化圖使目標曲線曲率單調均勻變化的曲線光順方法，該算法是從曲線的光順準則出發(fā)來對曲線進行處理。

觀察圖6可知，相比于EMD法、曲率法，從曲線光順效果上看，本文算法效果最優(yōu)。其中，在曲線拐點處的光順處理上，本文算法比EMD法、曲率法的光順效果更好。因此，初步觀察，直觀上可知本文提出的算法具有一定的優(yōu)勢。

為進一步驗證本文算法的有效性，算出原始曲線的曲率以及各個方法的曲率，其曲率圖如圖7所示。

分析圖7可知，本文提出的光順方法所得的曲線曲率變化最為平緩，且曲率值較低。本文方法與EMD法相比，光順效果略優(yōu)，與曲率法相比，優(yōu)勢較為明顯。

為了更加客觀地評價3種方法的優(yōu)劣性，現(xiàn)計算出各個曲線的曲率值，其曲率結果分析如表3所示。

表3 各方法曲率結果分析

分析表3可知：從最大曲率上看，本文算法所得最大曲率最小，比原始曲線小了2.426 1，比EMD方法小了0.242，比曲率法小了0.100 5；從最小曲率上看，本文算法所得最小曲率最大，比原始曲線大了0.010 4，比EMD法大了0.009，比曲率法大了0.006；從平均曲率上看，本文算法所得平均曲率最小，比原始曲線小了0.076 1，比EMD法小了0.002 6，比曲率法小了0.021 9。綜合分析最大曲率、最小曲率和平均曲率可知，本文算法在保留曲線有效特征的同時，也有效地去除了曲線的噪聲。相比之下，本文算法的光順效果最優(yōu)。

從運行時間上看，EMD法運行時間最短，本文算法次之，曲率法運行時間最長。本文算法運行時間雖慢于EMD法，但本文方法是多次EMD方法的疊加，因此其運行時間慢于EMD法。雖然運行效率略低，但曲率值變化更為平緩，效果更優(yōu)。

表3雖然對上述3種方法的曲線曲率進行了對比，但本文算法僅針對曲線進行光順，相反曲率法還可以針對曲面進行光順，因此，本文算法還存在一些局限，只能對曲線進行光順。

4 實例分析

為了進一步驗證本文算法，了解本文算法在工程實際中的應用，現(xiàn)給出以下3個實例分析。

實例1圖8是帶噪聲的離散數(shù)據(jù)光順，它是在圖8a的模擬數(shù)據(jù)中加入了高斯噪聲。圖8a為局部加入了0.02 mm的高斯白噪聲的初始截面數(shù)據(jù)；圖8b是經(jīng)由本文算法處理后得到的光順結果。觀察圖8可以發(fā)現(xiàn)，經(jīng)由本文算法處理后，加入高斯噪聲的模擬數(shù)據(jù)取得了不錯的光順效果。

實例2圖9為某汽車車身部分截面數(shù)據(jù)的光順，在驗證本文算法的重構及光順效果有效性的過程中，選擇圖9a中畫虛線內的汽車后蓋輪廓截面數(shù)據(jù)作為驗證的原始數(shù)據(jù)。圖9a是汽車車身的截面數(shù)據(jù)。圖9b是汽車后蓋輪廓截面數(shù)據(jù)的曲率梳。觀察圖9，通過對比汽車后蓋處理前后的曲率梳圖，顯然，經(jīng)由本文算法的處理，汽車后蓋得到了很好地光順。

實例3圖10為某渦輪葉片榫頭部分截面數(shù)據(jù)的光順，在驗證本文算法的重構及光順效果有效性的過程中，選擇圖10a中畫虛線內的部分榫頭截面數(shù)據(jù)作為驗證的原始數(shù)據(jù)。圖10a為渦輪葉片榫頭的截面數(shù)據(jù)，數(shù)據(jù)采集密度D≤0.02 mm。圖10b為渦輪葉片榫頭部分數(shù)據(jù)的曲率梳圖，圖10c為本文方法光順后的曲率梳圖。觀察圖10可以看出，經(jīng)過本文算法的處理，渦輪葉片榫頭得到了相應程度地光順。

5 結束語

本文利用集合經(jīng)驗模態(tài)分解、游程檢測法重構和奇異譜分析降噪三種方法相結合，提出一種將空間離散數(shù)字曲線看成離散數(shù)字信號的曲線光順算法，其核心部分在于：

(1)使用EEMD對空間離散數(shù)字曲線進行多尺度分解，一方面解決了EMD分解后出現(xiàn)的邊界效應和模態(tài)混疊問題；另一方面有效降低了原始空間離散數(shù)字曲線的x,y,z三個變量序列之間的相互干擾，以及有效降低了其非平穩(wěn)性；同時較為深入地挖掘了空間離散數(shù)字曲線的信息。

(2)由于EEMD分解時，所得到的分量過多導致光順工作量變大，采用了游程檢測法重構EEMD分解后的分量。在減少計算量、降低誤差的同時，可以更加直觀地反映原始空間離散數(shù)字曲線的特征。

(3)在對重構后的分量進行處理時，考慮到原始空間離散數(shù)字曲線的噪聲主要集中在高頻分量，直接刪除高頻分量會去除曲線的有效信息，采用了SSA對高頻分量進行降噪，從而避免了去除曲線的有效特征，使得曲線的光順盡可能符合原始曲線。

本文算法還存在一些不足：①在EEMD分解中，雖然在一定程度上克服了邊界效應和模態(tài)混疊的問題，但未能徹底解決；②在使用SSA降噪過程中，僅對高頻分量去噪。然而在實際情況中，還有一些噪聲存在于低頻分量中，本文未對其進行處理，致使光順后的曲線中還存留一些噪聲信息，采取何種方法對低頻分量進行處理，這是本文以后研究的方向。