張 俐， 吳中林

(北京航空航天大學機械工程及自動化學院，北京 100191)

基于點到曲面投影的飛機蒙皮波紋度計算方法

張 俐， 吳中林

(北京航空航天大學機械工程及自動化學院，北京 100191)

為了精確檢測飛機蒙皮表面的波紋度，提高檢測效率，研究了一種改進的點到參數(shù)曲面投影算法，即使用橢圓環(huán)面片逼近蒙皮設計曲面的方法求解點到曲面的投影，采用該方法從曲面的測量點云中提取待測波紋度位置的測量點，并計算測量點與理論位置的偏差值。根據(jù)偏差值的大小變化趨勢提取局部極值即波峰波谷值，并計算相應的波長，結(jié)合波紋度的定義完成波紋度計算。實驗結(jié)果表明，該方法分析結(jié)果準確、效率高。

波紋度；點到曲面投影；局部極值；飛機蒙皮

飛機蒙皮表面的波紋度會對飛機氣動特性和飛行性能造成較大影響[1]。在飛機的研制、試飛取證以及量產(chǎn)過程中都需要對飛機蒙皮表面的波紋度進行檢測。因此，在飛機外形容差規(guī)范中，對波紋度有嚴格規(guī)定。

目前，由于用實物樣條方法檢驗波紋度的精度和效率低等原因，各飛機制造商基本已經(jīng)放棄該方法，但是尚未有其他規(guī)范成熟的測量工具來精確分析波紋度[1]。隨著三維測量技術的發(fā)展，為通過軟件來精確分析飛機表面波紋度提供了可能性。

1 波紋度分析方法的思路

波紋度是用來描述介于零件表面粗糙度和形狀誤差之間的物體表面輪廓不規(guī)則變化[2]。工程中常用波幅和斜率來定量表示波紋度，如圖1所示。通常將相鄰兩波峰的距離稱為波長，用L表示；波

圖1 波紋度的定義

傳統(tǒng)的波紋度檢測方法主要有3種：鋼尺法、卡板法和特征樣條法。隨著三維數(shù)字化測量技術在航空工業(yè)領域的普及，波紋度的檢測方法有了很大改進。馬明等[2]提出應用小波濾波技術分析飛機蒙皮表面的波紋度。金鼎和吳劍[1]借助Max/Exa-scan測量系統(tǒng)采集飛機表面信息，結(jié)合 CATIA軟件對三維數(shù)字采集數(shù)據(jù)進行分析，最終獲得被測區(qū)域波紋度情況。于立明和李啟明[3]利用最小二乘擬合方法獲取曲面、曲線的波紋度數(shù)據(jù)。

本文波紋度計算是在傳統(tǒng)的波紋度檢測方法基礎上，結(jié)合三維數(shù)字化測量技術，改進點到參數(shù)曲面的投影算法來提取測量點并使用新的方法計算波紋度。

1.1 基準曲線的確定

波紋度反映的是曲線波動情況，目前波紋度檢驗通常采用的是將被測外形面按橫向或縱向劃分百分線的模擬量測量方法。因此，以此為依據(jù)，在理論曲面上確定一條基準曲線，即百分線。并通過比較測量點與基準曲線的偏差波動來分析計算波紋度。

已知飛機蒙皮理論參數(shù)曲面表達式為 S(u,v)，其中 u，v為曲面的兩個參數(shù)?；鶞是€是按如下方法確定的：

(1) 若給定理論曲面縱向波紋度檢測的百分比K，即在曲面兩條橫向邊界曲線上確定兩點P0、P1，使得P0、P1點到兩條橫向邊界曲線開始位置的弧長分別占各自邊界曲線總弧長的百分比為K。如圖2所示，；若給定橫向波紋度檢測百分比，則位置點P0、P1位置也按照上述方法確定。

(2) 將連接P0、P1點的直線投影到理論曲面上，如圖3所示，投影曲線Lp即為基準曲線。

圖2 P0、P1為縱向百分比K確定的兩個邊界點

圖3 Lp為基準曲線

1.2 實現(xiàn)流程

如圖3所示，將曲面上連接點P0、P1點的直線離散化為等間距的有限個點，將這些點投影到理論曲面上求得投影點，根據(jù)投影點能夠生成基準曲線Lp。已知基準曲線上的點，結(jié)合文獻[4]中有關方法，提取測量點并計算測量點偏差值。最后根據(jù)偏差值計算波長及波峰波谷值，完成波紋度計算。整個算法的流程圖如圖4所示。

圖4 波紋度計算流程

2 測量點的提取和偏差計算

根據(jù)直線P0P1參數(shù)方程，按一定的步長取直線上的點，其中，。再用分割逼近法[5]計算測量點 P0到理論曲面的最短距離的點 Q′(u′,v′)。以 Q′(u′,v′)為初始

在求得當前投影點T0和參數(shù)(u0,v0)以后，下一個投影點計算的初始值是當前投影點，以此類推計算點iP′在曲面上的投影點Ti和參數(shù)(ui,vi)，直到i=N。

2.1 點到曲面的投影計算

計算點到參數(shù)曲面投影的方法主要有兩大類：迭代法和細分法。迭代法又分為代數(shù)求根法和幾何迭代法。代數(shù)求根法的代表算法是Mortenson[6]提取的一種數(shù)值計算法以及Piegl和Tiller[7]提出的方法。幾何迭代法有Hu和Wallner[8]提出一種二階幾何迭代法；Liu等[9]提出的圓環(huán)面片逼近法。Piegl和Tiller[10]、Ma和Hewitt[11]、Selimovic[12]都提出了細分曲面的無效曲面片排除準則。

代數(shù)求根法對初始迭代點要求高，而細分法計算復雜。在幾何迭代法中Liu等[9]的方法計算效率相對較高，但是穩(wěn)定性有待進一步提高。所以本文對Liu等[9]的方法進行改進，用橢圓環(huán)代替圓環(huán)逼近曲面，計算效率和穩(wěn)定性都有提高。

2.1.1 橢圓環(huán)面片逼近法求投影點

如果將笛卡爾坐標原點放在橢圓環(huán)面中心，則橢圓環(huán)的方程式為

其中，R0為環(huán)的大半徑；ax和az為橢圓環(huán)橫截面橢圓的長、短半徑，如圖5所示。

圖5 坐標原點在橢圓環(huán)中心的示意圖

已知參數(shù)曲面S(u,v)，給定其參數(shù)域內(nèi)的一個參數(shù)(u0,v0)，計算該參數(shù)對應的點 Q=S(u0,v0)、曲面在該點的單位法向量n、兩個主曲率k1，k2和對應的單位主方向e1，e2。

首先要對上述幾何信息進行預處理，具體如下：

(2) 如果k1＞0，將n、k1和k2都乘以–1，改變其符號。這樣的預處理保證和k1≤0。

然后構造這樣的橢圓環(huán)面S′：中心在Q+n/k1，軸向為e2。橢圓的短軸，長軸。當k2＜0時，中心圓的半徑為；當k2＞ 0時，中心圓的半徑為。橢圓環(huán)面S′和原曲面S在點Q處有相同的法向、主曲率和主方向。根據(jù)前面所做的預處理，這個橢圓環(huán)面是唯一確定的。在本文的逼近方法中，采用橢圓環(huán)的局部面片來逼近曲面。

在構建原曲面的逼近橢圓環(huán)面以后，將測試點T投影到橢圓環(huán)面片上，求得投影點T′，然后求T′在原曲面上的參數(shù)。對于上述的橢圓環(huán)面片，求點到其投影可分為以下兩步：

(1) 將測試點T投影到中心圓上，求得橫截面橢圓的位置。

(2) 將測試點T投影到橫截面橢圓上，求得投影點T′。

在上述第二步中，測試點和橫截面橢圓在同一平面上。求點到平面橢圓的投影時，建立一個關于未知數(shù)為θ的一元四次方程。由于測試點和橢圓位于一個平面上，所以方程一定有實數(shù)解，而且有兩個解，其中到測試點最近的點作為投影點T′。

分別用本文介紹的橢圓環(huán)逼近求點到曲面方法和文獻[9]中的圓環(huán)逼近法測試文獻[8]中的實例1，算法的收斂準則都采用Piegl和Tiller[7]的標準，即

本文中當 ε1≤ 10-10時，迭代終止。表1列出了當測試點為(120, 10, 100)，初始參數(shù)值為(0.9, 0.6)時，每次迭代的參數(shù)增量Δu、Δv和迭代次數(shù)。表2中列出了初始參數(shù)值為(0.4, 0.1)時的情況。

表1 測試點(120, 10, 100)初始值(0.9, 0.6)的迭代次數(shù)和參數(shù)增量表

表2 測試點(120, 10, 100)初始值(0.4, 0.1)的迭代次數(shù)和參數(shù)增量表

由表1可以看出，本文采用的橢圓環(huán)逼近法只需要迭代 4次就能達到迭代終止條件，而文獻[9]中使用的圓環(huán)面片逼近法則需要迭代5次。表2中，在初值參數(shù)值與真實參數(shù)值的差值較大情況下，本文的方法也能比文獻[9]使用的方法少迭代1次，迭代效率分別提高25%和14%。

2.1.2 基準曲線上點的求取

將通過使用分割曲面法[5]求得的點Q0(u′,v′)作為初始迭代點，在該點構造橢圓環(huán)面片，并將起始點投影到橢圓環(huán)面片上，求得點，再用Newton迭代法求點在曲面S上對應的參數(shù)(u1,v1)和坐標值T1。算法框架如下：

(1) 計算曲面在初始迭代點處的幾何信息，包括位置Q0，法向量n，兩個主曲率k1，k2和相應的主方向e1，e2。

(2) 應用這些幾何信息，在Q0處構造一個橢圓環(huán)面片ST。

(4) 用這個新的參數(shù)作為初值，重復步驟(1)～(3)，直到求得的參數(shù)滿足必要的精度。圖6是以Q0初始迭代點，點在橢圓環(huán)面上的投影點。

圖6 點在橢圓環(huán)面上的投影

2.2 測量點的提取和偏差計算

根據(jù)求得的投影點Ti以及對應的曲面參數(shù)(ui,vi)，其中 i= 0,1,… ,N ，計算曲面S在點Ti處的法向量 ni(ui,vi)。在測量點云中搜索到法向量ni距離最近的測量點，并記最短距離為di，若di小于給定的閾值則提取該測量點并計算測量點到 Ti點的距離在法向的投影值Di，把該值作為了測量點到曲面的偏差；否則排除，如圖7所示。

圖7 根據(jù)曲面上點Ti提取測量點Ci并計算偏差Di

3 波紋度計算

由于文獻[1-3]所述方法的不足，本文使用一種新的波紋度計算方法，即在完成測量點的提取以及每個測量點的偏差計算后，通過依次比較偏差值大小的方法來分析波紋度。

根據(jù)波紋度的定義，波峰、波谷為測量點到曲面距離的局部極值，即測量點偏差的局部極值；波長為相鄰波峰(或波谷)對應基準曲線的弧長。

記得在紀錄片也看到過這樣的故事，冰島國家破產(chǎn)，金融危機階段，許多人去了鄰國挪威找工作，然而過了些年，大家又回來，引起了冰島當?shù)孛襟w的好奇，接受采訪，有人說是對冰島經(jīng)濟又有了信心，克朗貶值，旅游業(yè)大力發(fā)展，也有一群人的答案整齊統(tǒng)一：為了暖氣！挪威燒煤供暖，費用驚人，難怪家家戶戶有古老的壁爐，寧可燒木柴，缺點是只有客廳熱，要是臥室睡覺，或廁所洗澡，冷得哭出聲。

在計算波紋度之前，先對測量點的偏差Di進行如下預處理：

(1) 將測量點的偏差按照之前計算順序依次進行前后比較，后者較大的值用1表示，否則以–1表示。將比較的結(jié)果存入數(shù)組COMP中。

(2) 依次統(tǒng)計數(shù)組COMP中每個1、–1的連續(xù)個數(shù)，結(jié)果保存到數(shù)組LinkNum中，并將每個連續(xù)的開始點在測量點中對應位置存在數(shù)組Index中。

波峰、波谷值提取方法流程圖如圖8所示。

已知測量點偏差值和基準曲線上的點 Ti是一一對應的關系，根據(jù)偏差值的極值就能確定極值在Ti中對應點。相鄰波峰(或波谷)對應Ti中兩點之間距離即為波長。

圖8 波峰、波谷值提取方法流程圖

4 實例驗證

根據(jù)上述方法開發(fā)波紋度檢測程序，用機翼上蒙皮進行波紋度檢測驗證。圖9給出了翼面各百分比位置的波紋度檢測實例。

理論曲面如圖9(a)所示。對原始測量點云進行去噪、簡化等預處理，并導入程序中，顯示結(jié)果如圖9(b)所示。然后對點云進行坐標變換，提取縱向62%位置的測量點，如圖9(c)所示。測量點的偏差曲線如圖9(d)所示，其中測量點橫坐標為該點在基準曲線Lp的弧長位置，縱坐標為偏差值。計算的波紋度值如圖9(e)，其中每一個矩形框代表一個波紋度值，矩形框的縱坐標為波紋度值，矩形框的長度為波長，矩形框在橫坐標中的位置為波長在基準曲線弧長中的位置，檢測波紋度的標準參考文獻[13]。

圖9 實驗過程

對比圖9(d)和(e)，圖9(e)中波紋度的分布情況符合其曲線波幅分布情況。在圖9(f)中，分別對縱向50%、75%和橫向25%、65%位置進行波紋度計算結(jié)果。通過圖9的實驗驗證，本文能快速、有效地分析曲面的波紋度，分析效率高。

5 結(jié)束語

本文在總結(jié)傳統(tǒng)波紋度檢測經(jīng)驗基礎上，結(jié)合先進的數(shù)字化檢測設備，使用投影算法和波峰、波谷提取方法分析波紋度。該方法通過實驗驗證，并實際運用于飛機蒙皮波紋度檢測。

從2.1節(jié)的實驗數(shù)據(jù)可以看出，同一測試點，在迭代效率上，本文使用橢圓環(huán)面片逼近法求點到參數(shù)曲面的投影方法，比文獻[9]圓環(huán)面片逼近法求投影的方法迭代效率分別提高25%和14%。

通過使用新的波紋度分析方法，能夠?qū)崿F(xiàn)波紋度的快速分析計算，不依賴測量儀器，也不需要人工提取波峰、波谷點。只需要給定整個曲面的測量點，就可以實現(xiàn)曲面各個位置波紋度的檢測。本文的方法比金鼎使用的方法更自動化。

需要指出的是，由于測量點的疏密直接影響波紋度計算的結(jié)果，所以為保證波紋度分析的準確性，在保證計算速度的前提下應提高測量點的密度。

[1] 金 鼎, 吳 劍. 基于三維測量技術的飛機表面波紋度數(shù)字化分析方法[J]. 航空工程進展, 2013, 4(1): 85-89.

[2] 馬 明, 曹喜鋒, 李啟明, 等. 基于小波濾波的民用飛機蒙皮表面波紋度分析方法[J]. 科學技術與工程, 2013, 13(35): 10750-10754.

[3] 于立明, 李啟明. 基于最小二乘法的飛機蒙皮波紋度處理方法研究[J]. 科技信息, 2012, (28): 424-427.

[4] 劉子路, 張 俐, 鄭國磊. 一種基于CAD的翼類曲面誤差檢測方法[J]. 圖學學報, 2012, 33(6): 136-139.

[5] 廖 平. 基于遺傳算法和分割逼近法精確計算復雜曲面輪廓度誤差[J]. 機械工程學報, 2010, 46(10): 1-7.

[6] Mortenson M. Geometric modeling [M]. New York: John Wiley & Sons, 1985: 505-525.

[7] Piegl L, Tiller W. The NURBS book [M]. 2nd ed. New York: Springer, 1997: 232-234.

[8] Hu S M, Wallner J. A second order algorithm for orthogonal projection onto curves and surfaces [J]. Computer Aided Geometric Design, 2005, 22(3): 251-260.

[9] Liu X M, Yang L, Yong J H, et al. A torus patch approximation approach for point projection on surfaces [J]. Computer Aided Geometric Design, 2009, 26(5): 593-598.

[10] Piegl L A, Tiller W. Parametrization for surface fitting in reverse engineering [J]. Computer-Aided Design, 2001, 33(8): 593-603.

[11] Ma Y L, Hewitt W T. Point inversion and projection for NURBS curve and surface: control polygon approach [J]. Computer Aided Geometric Design, 2003, 20(2): 79-99.

[12] Selimovic I. Improved algorithms for the projection of points on NURBS curves and surfaces [J]. Computer Aided Geometric Design, 2006, 23(5): 439-445.

[13] 何文治. 航空制造工程手冊:飛機裝配[M]. 北京: 航空工業(yè)出版社, 2010: 14-15.

A Calculation Method of Aircraft Surface Waviness Based on the Method of Point Projection on Curved Surface

Zhang Li, Wu Zhonglin

(School of Mechanical Engineering and Automation, Beihang University, Beijing 100191, China)

To measure the waviness of aircraft surface accurately and efficiently, an improved algorithm for the projection of points on parameter-curved surface is proposed. The point projection on curved surfaces is calculated by adopting the approach of elliptical torus patch approximation. Firstly, by using this algorithm, points are extracted from the surface’s massive measuring points and the offset value between measuring points and their projection on the theoretical positions is calculated. Secondly, the local extremum is extracted at the peak and valley of waviness according to the variation trend of the offset value. Then its corresponding wavelength is calculated. Finally, the waviness is calculated based on its definition. Experiments show that this method is highly efficient, and generates accurate results.

waviness; point projection on surfaces; local extremum; aircraft surface

TP 391

10.11996/JG.j.2095-302X.2016050626

A

2095-302X(2016)05-0626-06

2015-12-23；定稿日期：2016-03-19

基于工藝驅(qū)動的全機全三維綜合測量系統(tǒng)開發(fā)與集成技術項目(2013E11222)；飛機后段自動精準對合技術與成套設備項目(2014ZX04001-081-07)

張 俐(1961–)，女，北京人，高級工程師，學士。主要研究方向為飛機數(shù)字化裝配檢測技術。E-mail：gracejune@buaa.edu.cn