基于改進凸包算法的葉片型面特征參數提取

（華中科技大學數字制造裝備與技術國家重點實驗室，湖北 武 漢 430074）

葉片在航空航天、汽車、能源等領域，都有非常廣泛的應用。據有關專家預測，到2025年，新增民用飛機（客機和貨機）將達到2.72萬架，全球機隊規模將翻一番；到2020年，國內通用航空飛機將超過1萬架。巨大的國際、國內市場，使得葉片的需求量迅速增大。

葉片型面一般為復雜曲面，加工工序比較復雜，其型面品質對發動機的性能起著決定性的影響[1]。如何檢測和評估葉片的型面品質，已經成為航空制造領域的一個關鍵性技術問題之一。

目前，國內外主要采用人工卡板測量法和三坐標測量機（CMM）測量法，來對葉片型面進行檢測。這兩種方法都是離線式檢測方法，檢測速度比較慢。用激光掃描儀、結構光測量儀等非接觸式測量設備進行葉片型面的檢測，可以克服傳統人工卡板測量法和CMM測量法的不足，并且實現葉片型面的在機檢測，其還將成為葉片類復雜曲面型面精密檢測的重要手段。利用非接觸式測量得到葉片的點云數據，高效快速地提取葉片型面特征參數，具有重要的現實意義和實用價值。

如何根據測量得到的葉片數據，準確提取出葉片型截面凸包點云，是提取葉片型面特征參數的關鍵步驟之一。目前帶有葉片檢測模塊的軟件如Geomagic Qualify和Polyworks等均未提供葉片葉身軸線方向，葉片型面的截面選擇具有任意性，從而使得提取的結果并不能真正體現葉片型面的特征。經典的凸包算法包括卷包裹（Jarris）法、格雷厄姆（Graham）方法以及分治算法等。這些算法在處理海量點云數據時，存在效率低下、精度難以保證等諸多問題。

劉人午提出先將數據點集進行一次掃描，得到橫向和縱向排序點表并建立初始凸包，再運用增點法從外向內判別數據點是否加入凸包表的改進算法，成功解決了數據量超大時計算效率低下問題[2]。

但該算法存在兩處不足之處：一是直接將大量的原始點集進行排序非常耗時，二是其判別算法結束的條件不具有通用性，從而使得到的結果不準確。

王文軍提出以葉片截面在其最小外接圓上的各點之間距離最大者作為弦長[3]，由于沒有在弦線上投影變換，使得計算結果不準確。

俞學蘭在《基于MATLAB的葉片參數辨識》一文中詳細講解了求弦長的方法[4]，但該方法比較繁瑣，且不能用于前后緣是橢圓型葉片弦長的計算中。

本文提出了基于主成份分析法的葉片葉身軸線方向的提取方法，并以此軸線方向為葉片型截面的法矢，準確截取葉片型面，得到葉片型截面點云數據。

本文改進了文獻[2]的算法，用矩形區域腐蝕法確定凸包邊界區域，刪除矩形區域內的大量冗余點，大大提高了算法效率。

在葉片型面檢測中，本文用改進的凸包算法來對測得的葉片型面點云數據進行排序，提取葉片型面關鍵特征參數。

在提取葉片前后緣參數時，傳統的方法是按前后緣是圓弧的形式進行提取，但是近年來采用非圓弧形如橢圓形前緣可以明顯改善葉片的氣動性能[5]。

本文基于最小二乘法和隨機原理，提出了一種用橢圓一般方程擬合葉片前后緣的抗噪聲的高精度算法。最后，應用MATLAB軟件開發了用于葉片型面特征參數提取的軟件模塊。

1 凸包算法的改進

凸包問題是計算幾何的基本問題之一，在計算機圖形學、圖像處理、模式識別以及CAD、CAM等領域，都有著廣泛的應用[6]。

構建凸包通常要解決兩個問題：其一是要從大量的離散點中判斷出哪些是符合要求的凸包頂點；其二是要解決這些點的拓撲連接關系。

文獻[2]求取凸包最耗時的地方，是對點集快速排序與生成凸包。如果能夠減少參與排序和生成凸包的點集數量n，就可大大提高算法的效率，本文基于該思想，提出用矩形區域腐蝕法確定凸包邊界區域，以刪除矩形區域內的冗余點。

1.1 針對文獻[2]算法做出的改進

本文根據凸多邊形的邊同側特性[7]，針對文獻[2]的算法做出如下改進：

（1）用矩形區域腐蝕法確定凸包邊界區域，刪除矩形區域內大量的冗余點，形成初始的平面點集；

（2）找出初始平面點集中的極值點構成初始凸包；

（3）將初始平面點集中的點進行掃描得到按x、y坐標值從小到大排序的橫向排序表H（x值相同，按y值從小到大排序）和縱向排序表V（y值相同，按x值從小到大排序）；

（4）在每次調入H、V點集中第一點和最后一點進行判別后，要重新生成H、V點集，并判斷這兩個點集是否為空，若全為空時，結束判斷，生成凸包點集。

矩形區域腐蝕法原理如下：

先找出點集中的4個邊角點，即

A點（x+y值最小的點），

B點（x–y值最大的點），

C點（x+y值最大的點），

D點（y–x值最大的點）

然后用

y ＝ max（Ay，By），

y ＝ min（Cy，Dy），

x ＝ max（Ax，Dx），

x ＝ min（Bx，Cx）

畫直線圍成矩形區域，最后用腐蝕法將點集中落于這個矩形區域的點腐蝕掉，剩下凸包邊界區域點，即初始平面點集，如圖1所示。

圖1 矩形區域腐蝕法示意圖

1.2 初始凸包生成方法

（1）找出點集中的極值點。

A1為x坐標等于x min的點中y坐標最小的點；

A2為x坐標等于x min的點中y坐標最大的點；

A3為y坐標等于y max的點中x坐標最小的點；

A4為y坐標等于y max的點中x坐標最大的點；

A5為x坐標等于x max的點中y坐標最大的點；

A6為x坐標等于x max的點中y坐標最小的點；

A7為y坐標等于y min的點中x坐標最大的點；

A8為y坐標等于y min的點中x坐標最小的點。

（2）將這些極值點按順時針方向的順序，寫入數組A中，即

A＝[A1；A2；A3；A4；A5；A6；A7；A8]

(A為8×2的數組)。

（3）刪除A中相同的點并將A中的第一點復制增加為A的最后一點，使其首尾相連，形成初始凸包如圖2所示。

圖2 初始凸包生成示意圖

1.3 定義點與有向直線的關系判別函數

定義點與有向直線的關系判別函數S如下：設

A(xa、ya)、B(xb、yb)和 C(xc、yc)

為XY平面內任意不同的3點，記

L（AB）為A指向B的有向直線，用三點面積S判斷點C與L(AB)的側向關系，其判別函數為

S為正時，表示A、B、C是逆時針；

S為負時，表示A、B、C是順時針。

（1）如果 S＞0，點 C 在 L（AB）的左側；

（2）如果S＝0，點C在L(AB)上；

（3）如果S＜0，點C在L(AB)的右側。

如果當前判斷點在當前凸包外(存在一條邊使得該點在該邊的左邊)，則將該點插入凸包表中該邊的兩端點之間[2]。這種說法欠妥，如果這個點在當前凸包相鄰兩邊交點的外頂點區域，增加到其中任一邊兩端點之間后，會使凸包發生畸變，改變了這個多邊形的凸包屬性。

改進如下：

設凸包兩相鄰邊為AB、BD，當前判斷點C，

若 S（A，B，C）＞0

且S（B，D，C）≤0，C點增加到邊AB 兩端點之間；

若 S（A，B，C）＞0

且 S（B，D，C）＞0，去掉 B 點，C 點增加到邊 AD兩端點之間。

1.4 凸包算法流程及其效率分析

凸包算法流程如圖3。

圖3 凸包算法流程圖

下面進行該算法效率的分析：

該算法主要分4步：

第一步是用矩形區域腐蝕法求原始平面點集（設點數為n）得初始平面點集（設點數為n1），n1大約是n的1/10，其時間復雜度為O（n）；

第二步是提取初始平面點集中的極值點生成初始凸包，其時間復雜度為O（n1）；

第三步是對去除初始凸包點后的初始平面點集（設點數為n2）進行雙向排序生成H、V集，若采用時間復雜度逼近于O（n）的排序方法時[8]，其時間排序的時間 O（n2）；

第四步根據判別函數S，在初始平面點集中提取凸包點生成凸包，其時間復雜度為O（n2）。

綜合起來，該算法的時間復雜度為O（n），這與目前平面點集凸包問題的漸進最好算法的時間復雜度為O（n log h）相比，下降了一個數量級，具有很好的時間效率[9]。

在2G內存2.5 GHz主頻的電腦上，利用MATLAB軟件編程對上述算法進行實驗分析。

在用矩形區域腐蝕法刪除矩形區域內部點效率分析的實驗中，點集的坐標由隨機函數（rand）產生，在實驗中分別設置點集數量為103、104、105、106、1.5×106和1.5×107，實驗結果如表1所示。實驗表明，大約有90%以上的冗余點集，可以通過本算法刪除。

表1 用矩形區域腐蝕法刪除矩形區域內部點效率分析

在凸包算法運行結果比較實驗中，點集的坐標由隨機函數（rand）產生，文獻[2]和本文的算法運行結果如表2，可以看出本文算法所費時間比文獻[2]少一個數量級，原因在于本文在進行初始平面點集排序前運行矩形區域腐蝕法去除了矩形區域內的大量冗余點，使其點數銳減為原來點數的3‰左右，從而大大提高了程序的運行效率。

表2 凸包算法運行結果

2 葉片型截面凸包點云數據的提取

為了體現出葉片型面特征，不能隨意用一個平面去截取葉片型面，必須使截平面垂直于葉片葉身的軸線，即截面的法矢與軸線方向一致，以確保準確提取出葉片型面的特征參數。

在接觸式測量中，可以精確調整坐標軸，使其Z軸與葉片葉身軸線一致（如圖4），這樣在測量的葉片點云數據中鎖定Z軸，根據具體的精度要求，以一定的間隔△Z移動垂直于Z軸的截面去截葉片型面，即可得相應的葉片型截面點云數據（如圖5)。

圖4 CMM測量葉片

圖5 葉片截面示意圖

在非接觸式測量中，由于測量儀坐標軸方向隨機，葉片葉身軸線與坐標軸不一致，我們可以通過對測量的葉片點云數據進行主成份分析找到葉片葉身的軸線方向[10]，并將葉片點云數據通過坐標轉換，使其新坐標系的Z新與葉片葉身軸向一致，然后根據精度要求以一定的間隔△Z新移動垂直于Z新軸的截面去截葉片型面得到相應的葉片型截面點云數據。下面介紹一下基于主成份分析法提取葉片葉身軸向的方法。

主成份分析法是一種降維的統計方法（Principal Components Analysis，PCA），也是一種數學變換的方法，其借助于一個正交變換，將其分量相關的原隨機向量，轉化成分量不相關的新隨機向量，這在代數上，表現為將原隨機向量的協方差矩陣變換成對角形陣；在幾何上，表現為將原坐標系變換成新的正交坐標系。

經過特征矢量矩陣的轉置矩陣作為變換矩陣變換后，新的坐標軸是沿數據分布的主要分量方向即特征矢量的方向。對于二維空間，如果數據分布是一橢圓，那么特征矢量正好在橢圓的兩個軸上，U軸方向是第一主成份，V軸是第二主成份方向，如圖6所示。

圖6 橢圓主成份分析示意圖

這樣，將原相關變量集在特征向量空間的投影，即得不相關變量集，其特征向量矩陣為投影變換矩陣，也為XOY坐標系到UO1V坐標系的變換矩陣。

利用主成份分析法提取葉片葉身軸向的一般步驟為：

（1）根據葉片點云數據

計算其數學期望

協方差矩陣

（2）對協方差矩陣ΣX用特征值分解。即

ΣX× V=V×D，

D為ΣX的特征值對角陣；

V為ΣX的特征向量矩陣；

從而計算變換矩陣A=VT和主成份矢量方向。

根據葉片葉身長度與葉片型截面的長寬粗略比較：

若葉身長度是最長的，則選第一主成份（特征值最大）矢量方向為葉身軸向；

若為第二長，則選第二主成份（特征值第二大）矢量方向為葉身軸向；

若為第三長，則選第三主成份（特征值最小）矢量方向為葉身軸向。

（3）用A×X即得葉片點云數據在新坐標系下的點云數據X新。

在MATLAB中，改進的凸包算法，能高效處理測得的葉片點云數據：獲取按逆時針排序的葉片型截面凸包點云數據；將余下的凹點按x坐標值的升序排列求得葉片凹弧點云數據。這樣就將葉片型截面點云數據就按逆時針有序排列了（見圖7），為后面的葉片型截面數據處理作好了準備。

圖7 葉片型面凸包點云數據

3 葉片型面特征參數的提取

葉片的型面檢測通過型面特征參數來表達,這些參數主要包括弦長及弦傾角、前后緣點、中心、半徑及厚度、型面最大厚度、中弧線等，具體定義如下（圖8所示）。

圖8 葉片型面特征參數示意圖

弦線：葉片型截面上與葉盆的前后部型線同時相切的切線。

弦長：整個葉片截面輪廓在弦線上投影的長度。

軸弦：整個葉片截面輪廓在發動機軸心線上投影的長度。

弦傾角：弦線與發動機軸心線的夾角即為弦傾角。

前后緣點：分別是中弧線延拓后與前緣曲線和后緣曲線的交點。

前后緣中心：前后緣擬合后圖形的幾何中心，若前后緣是圓形則為圓心。

前后緣半徑：對于圓弧形緣頭，即為緣頭所在圓的半徑；其他類型的緣頭不存在此參數，如橢圓就是長短軸。

前后緣厚度：葉片型面中弧線上距離前后緣點1 mm處點為圓心的葉型內切圓直徑。

中弧線：由一系列葉身截面內切圓圓心定義的曲線。

最大厚度：一系列葉身截面內切圓中的最大內切圓的直徑。

限于篇幅，本文僅對葉片型面弦線及前后緣關鍵特征參數提取算法進行詳細介紹，對其他葉片型面關鍵特征參數，根據其定義在MATLAB中也可方便求解。

3.1 弦線及其參數提取

弦線為葉片型截面上與葉盆的前后部型線同時相切的切線。弦線提取技術的關鍵，是獲取弦線與葉盆前、后緣相切的切點。根據前面提取的葉片型截面凸包數據特點，可以認為凹線的兩端點即葉盆與葉背的交點為待求的兩切點。弦線提取方法如下：

（1）采用凸包算法獲取輪廓數據點的凸包；

（2）逐個計算凸包中相鄰兩點的距離，獲取距離最大的相鄰兩點；

（3）過該兩點作直線，即為弦線，該弦線斜率的反正切值，就是弦線與X軸（設X軸方向為發動機軸心線）的夾角，即弦傾角。

弦長為整個葉片截面輪廓在弦線上投影的長度。對于前后緣不是圓形，是其他圖形如橢圓，采用俞學蘭法會過于復雜，因此研究一種通用的計算弦線參數的方法具有重要意義。

根據弦長的定義，弦長的計算可以將測得的點云數據往弦線方向投影，求其在新坐標系下的坐標，然后求在新坐標系X坐標最大的點與X坐標最小的點的距離即為弦長（簡稱投影變換法）。投影變換法的原理如下：

假設弦線L通過平移與X軸交于原點O，夾角為弦切角α。通過繞Z軸逆時針旋轉α，使X軸與弦線L重合如圖9所示。

圖9 坐標轉換示意圖

同一點P在原坐標系OXY下的坐標為P0（x，y），在新坐標系O1X1Y1的坐標為P1（x1，y1），則P0和P1具有如下變換關系：

弦長 L＝max（P1x）– min（P1x）

(即葉片點云數據轉化到新坐標系后的x坐標最大值與x坐標最小值之差)。

幾種方法計算弦線參數結果如表3所示，由表可知，本文提出的投影變換法計算弦線參數可靠，且能通用于任意形狀的葉片截面。

表3 幾種方法弦線參數計算結果 單位：mm

3.2 前后緣參數提取及葉片型面精確分割

由于圓是橢圓的特殊情況，所以可以直接設葉片前后緣類型為橢圓弧，根據橢圓的一般方程建立模型

（如果A＝B則為圓；如果A≠B則為橢圓）。

基于代數距離最小二乘法的常規橢圓擬合法將所有樣本點都當作準確值來擬合[11~13]，但當樣本點中出現雜質點時，擬合出的橢圓誤差較大，不能滿足精確度要求較高的系統需求，因此在擬合前后緣參數前，要精確提取出前后緣上的點。

（1）粗略確定前后緣的范圍。在前面提取凸包后，將弦線兩端點沿前后緣方向上取相對精確的6個點，代入方程（1）得到擬合橢圓的相關參數，計算葉片型截面點云數據到這個橢圓方程的歐式距離d，若d≤△d（設定的閥值），則認為該點在橢圓上，從而得到葉片前后緣的粗略范圍（記為樣本點云）。

（2）精確確定前后緣上的點云數據。在第一步求取的前后緣點中，有可能含有噪聲，如何除去前后緣上的噪聲，即不在橢圓上的點，本文基于最小二乘法和隨機原理，提出了一種能有效去除前后緣噪聲確定前后緣精確數據點的方法。

其算法原理如下：

其一，在初步確定的前后緣范圍內的樣本點云（已編號）隨機取6個點。

其二，將這6個點代入

Ax2+By2+Cxy+Dx+Ey+F=0中，

求得 A、B、C、D、E、F。

其三，計算樣本點云中所有點到前面6個點擬合成的橢圓的代數距離，即

D＝｜Ax2+By2+Cxy+Dx+Ey+F｜。

若D≤△D（設定的閥值），將這個點的編號記入數組Dianhao中，然后將樣本點云中滿足要求的點的總個數記入Dianshu中。

其四，比較匹配點總個數（Dianshu）與匹配點數最大值（Maxdianshu），當前者大于后者，將橢圓參數(A、B、C、D、E、F)和記錄匹配點編號的數組（Dianhao）保存下來,分別拷貝至數組Canshu和數組Maxdianhao，最后將Dianshu賦值給Maxdianshu。

其五，循環執行步驟其一到其四一定次數（根據運行時間、需要結果的準確度以及樣本點總個數適當定義），最后在Canshu保存了最優橢圓參數，在數組Maxdianhao保存了在所有樣本點中匹配點的編號，也就可以相應地得到不匹配點的編號。

算法流程圖如圖10所示。根據Maxdianhao中保存所有樣本點匹配的編號，求得前后緣精確的點。

圖10 前后緣參數算法流程圖

算法實驗：將方程為(x2)/9+(y2)/25＝1的橢圓離散成63個點，設定循環次數為N次，D的閥值為△D，在不同的高斯噪聲（0，σ2）下，本算法提取到符合要求的點數為n1，橢圓的5個參數及運行時間。

在不同的噪聲（即改變σ2）下，擬合結果如表4、表5所示，從表可以看出在N＝200，△D＝0.001情況下，除噪聲效果最好。

經過實踐證明，在求取最優橢圓的過程中，剔除了絕大多數的雜質點，相比傳統方法，準確度得到一定提高，系統魯棒性得到增強，且所需運行時間較短。

表4 本文方法對含有噪聲數據的擬合結果(N=100,△D=0.01)

表5 本文方法對含有噪聲數據的擬合結果（N＝200，△D＝0.001）

（3）精確確定前后緣的參數。在確定前后緣數據點后運用橢圓擬合的程序確定前后緣的精確參數，若橢圓的長短軸相差的絕對值不大于error(設定的閥值)，再運用圓擬合程序進行擬合得到圓弧形前（后）緣的精確參數（提取前后緣結果如圖11）。

圖11 葉片截面前后緣提取結果

（4）葉片型截面點云的精確分割。葉片型截面點云數據可分為如下4部分：葉盆、葉背、前緣、后緣。在精確確定前后緣點后，凸包上的點集減去在其上的前后緣點即得葉背部分點云，凹弧部分的點集減去在其上的前后緣點即得葉盆部分點云。

4 葉片型面特征參數提取的GUI設計

在研究葉片型面關鍵特征參數提取算法的基礎上，利用MATLAB軟件的GUI模塊，編寫了一個提取葉片型面特征參數的用戶界面，便于實現人機交互。

該界面主要包括葉片型面測量點云的輸入（*.dat或*.txt格式）、圖形顯示設置（如線型、線寬、網格、坐標軸等）、處理后結果的顯示（如文本輸出、各特征參數在圖形窗口中的顯示等）。

將某型號航空發動機葉片型面的激光掃描點云數據，導入到葉片型面特征參數提取的軟件模塊，提取結果如圖12所示。

圖12 葉片截面特征參數提取界面

本文方法提取的結果與Geomagic Qualify軟件提取結果（如圖13）的比較見表6。本文采用了時間效率高的凸包算法，并改進了葉片型面特征參數的提取算法，從表6給出的數據可以看出，本方法在精度上比這款軟件高。

另外，本方法在提取葉片型面特征參數方面比Qualify軟件有更好的操作性。

圖13 Qualify軟件提取葉片型面參數結果

表6 葉片截面特征參數提取結果比較

5 結束語

本文提出了基于矩形區域腐蝕法的改進凸包算法，用于提取葉片型面特征參數，如弦線、前后緣參數等，實驗結果表明，改進后的算法具有更高的計算效率和計算精度，通用性好，且具有一定的抗噪能力。主成份分析法用于分析和提取葉片截平面法矢，為準確提取葉片型面特征參數提供方向矢量。在Matlab軟件平臺下開發了面向葉片測量點云數據葉片型面特征參數提取的人機交互界面，并通過實驗數據驗證了該軟件模塊的可行性與有效性，為航空發動機葉片在機非接觸式檢測應用提供了基礎。

[1]陸佳艷，熊昌友，何小妹，等.航空發動機葉片型面測量方法評述[J].計測技術，2009，29(3)：1-3.

[2]劉人午，楊德宏，李 燕，等.一種改進的最小凸包生成算法[J].大地測量與地球動力學，2011，31(3)：130-133.

[3]王文軍.基于Surfacer葉片型面檢測模塊的開發[D].無錫:江南大學,2008.

[4]俞學蘭，葉佩青.基于MATLAB的葉片參數辨識[J].航空維修與工程，2009，(4)：56-58.

[5]張力寧，張定華.葉片前緣高精度重建方法研究[J].航空動力學報，2006，21(4)：722-726.

[6]李光軍，鄭軍紅，張光忠.一種凸包的改進算法設計與實現[J].現代計算機(專業版)，2010，(6)：92-94.

[7]劉光惠，陳傳波.求解簡單多邊形和平面點集凸包的新算法[J].計算機科學，2007，34(12)：222-226.

[8]Pigeon S.Quicksort and radix sorts on lists[J].Dr.Dobb's Journal，2002，27(5)：89-89.

[9]Kirkpatrick D G,Seidel R.Ultimate planar convex hull algorithm[J].SIAM Journal on Computing，1986，15(1)：287-299.

[10]白困利.基于主成份分析法的絕緣子外絕緣特性的綜合評判[J].華中電力，2010，23(2)：9-12.

[11]Fitzgibbon A，Pilu M，Fisher R B.Direct least square fitting of ellipses[J].IEEE Transactions on Pattern Analysis and Machine Intelligence，1999，21(5)：476-480.

[12]雷志術，張雁波.橢圓定形曲線擬合問題若干新型算法[J].上海交通大學學報，2002，36(8)：1210-1213.

[13]Halí R，Flusser J.Numerically stable direct least squares fitting of ellipses[C].The sixth International Conference in Central Europe on Computer Graphics and Visualization[A].1998，1(1)：125-132.