基于蒙特卡洛法的葉片型面參數測量不確定度分析

張學儀，何小妹，劉峻峰，王卓然，王一璋

（航空工業北京長城計量測試技術研究所，北京 100095）

葉片是航空發動機的關鍵零件，具有種類多、數量大、曲面形狀多樣、加工工藝復雜、加工難度大等特點，其質量直接影響到發動機的氣動性能與壽命[1–2]。隨著新型號發動機葉片型面在彎、扭和掠的增強以及前后緣曲率半徑的小型化，對葉片的幾何尺寸和形狀測量精度要求也越來越嚴格。由于葉片的尺寸、形狀、物理和力學等因素決定了葉片的整體性能[3–5]，葉片生產的不同階段均需要對葉片的幾何尺寸和形狀進行檢測[6]。

葉片型面參數描述了葉身幾何尺寸和各部分相互位置關系，參數種類主要包括葉型形狀、位置、扭轉、前后緣、長度等。總體來說葉片型面參數可以分為4類：尺寸類、位置類、角度類和輪廓度類。

航空發動機葉片型面的測量精度要求較高[5]，在對坐標測量機的葉片測量結果的準確性上有需求，這就需要對檢測結果進行不確定度評定[7–10]。現階段對于航空發動機葉片葉型測量不確定度評定尚未有統一的評定標準。一方面，由于坐標測量機在測量過程中影響測量結果不確定度的誤差來源復雜多樣，導致分析困難，各誤差源對測量結果的影響難以量化[11]，誤差模型建立難度大[12]。另一方面，對葉片這類自由曲面的復雜測量任務，葉片點云數據處理方法、參數評定算法的不同，都會對測量不確定度評定產生影響[13–15]。

而在現階段的研究中各國學者對于坐標測量機對特定測量任務的不確定度方面的研究，選擇采用了蒙特卡洛法，如Miura等[16]用蒙特卡羅模擬法估算基于坐標測量機的階梯規測量不確定度；宋愛國等[17]提出基于一種種群優化的微分進化算法評定了圓度測量誤差，并采用蒙特卡洛法評定出圓度測量的測量不確定度；意大利Ruffa等[18]研究了 CMM 進行圓度測量時，采樣點數和測量不確定度之間的關系；埃及Ali[19]討論了掃描式 CMM 在不同掃描速度和擬合算法的條件下，對直徑和圓度測量的影響；美國 Ramaswami等[20]建立了圓柱度測量時測量精度和測量策略之間的雙向估計模型。

可以看出，國內外對坐標測量機一些常規形位參數的不確定度測量研究成果較多且研究更有深度，在測量點數、測量策略等對形位測量不確定度的影響上做了大量的研究，但對于葉片這種自由曲面尚未有太多的不確定度評定方法。

因此針對三坐標測量機的葉片測量過程，本文提出一種基于蒙特卡洛法評定葉片型面參數測量結果不確定度的方法，并將本文的研究成果與基于GUM （Guide to the expression of uncertainty in measurement）法的部分參數評定結果進行了比較，證明了蒙特卡洛法葉片型面參數測量不確定評定方法的準確性和有效性。

1 基于蒙特卡洛法的葉片型面參數測量不確定度評定理論

應用蒙特卡洛法的測量不確定度評定方法是基于概率理論。利用蒙特卡洛法評定葉片型面參數的不確定度計算過程如圖1所示，x1，x2，…，xn為影響葉片參數誤差的主要因素，蒙特卡洛法首先要分析這些誤差源的測量估計值和測量不確定度，產生符合各參數分布狀態的隨機數；然后代入各參數誤差模型，構成f的概率分布；最后根據f的概率分布，求出其期望值和標準差，即為所求的葉片參數值和測量標準不確定度。

圖1 蒙特卡洛法評定參數不確定度Fig.1 Monte Carlo method for parameter uncertainty evaluation

針對坐標測量機測量葉片的過程，應用模型模擬測量任務并估計坐標測量機測量葉片的不確定度。模型的建立需要考慮測量點的分布、測頭的選擇和配置、測量系統所在環境等在模型中的影響表現。首先需要實際測量得到的點云數據，考慮在仿真過程中各誤差源不確定度和采用的其分布參數[12]，為所有測量點指定了1組仿真三坐標測量機的誤差特性[16]，通過選擇各誤差源的分布對實際測量數據的每個點進行擾動，生成多組測量值，并輸入到如尺寸參數（長度、厚度、角度等）、輪廓參數等實際參數數學模型計算表示，得到多組模擬測量參數值，并且由于選擇算法不同，通過對特定算法進行仿真可以得到同一參數不同算法產生的不確定度。

坐標測量機空間測量數學模型為：

根據誤差傳播定律，空間長度d的測量誤差可表達為：

其中，測量時，坐標測量系統返回的數據格式為（x，y，z，α，β，I，J，K），其中（x，y，z）為被測點的坐標值；α、β分別為測量時A軸與B軸的角度；（I，J，K）為測量時的探測矢量。坐標測量系統返回的是合成后的數據，并不能從控制系統中獲得各主軸的光柵值、測桿變形量以及測頭A、B軸運動產生的位移量。

通過對坐標測量機測量葉片測量精度的評估與坐標測量機精度的評估問題的分析，假設坐標測量系統模型的實質是分量誤差相互重疊，因此可以通過分析提取出測量過程中每個點的相互獨立的各分量誤差矢量。以試驗采用的測量范圍為2000mm×1200mm×1000mm的坐標測量機對某發動機葉片的測量過程為例，通過分析測量設備的不確定度、探觸形狀誤差引入的測量不確定度、重復性、溫度影響的特性指標，反映了葉片測量過程中不確定性的特征，反映各種誤差對于測量結果不確定度的影響，進行以下分析。

1.1 測量設備的不確定度

通過分析坐標測量機單軸測量不確定度、坐標垂直度誤差、坐標直線度誤差和旋轉誤差，最終合成得到測量設備的不確定度，最大示值誤差UMPE=（1.2+l/350）μm。假設在這個范圍內測量設備不確定度任何結果都是可能的，那么測量設備的不確定度服從矩形分布。

1.2 溫度變化對測量結果的影響

環境條件影響的估算在模擬測量模型中很重要。對于幾何量，機械零件的熱梯度會顯著影響后續的測量，影響測量尺寸的其中一個因素即為溫度變化。溫度的變化對測量結果的影響主要考慮三坐標測量機與工件之間的膨脹系數的差異及三坐標測量機標尺膨脹系數不確定度，因此有必要在考慮熱膨脹系數的情況下提供熱補償。根據ISO標準，幾何測量的參考溫度為20℃。與該溫度的任何偏差均應根據以下公式進行校正：

其中，lc為校正后長度；l為測得的長度；t為溫度；α為對于給定的校準元件材料采取的熱膨脹系數。

被測件和設備的膨脹系數差在半寬為2×10–6℃–1的區間內溫度Δt=0.3℃，服從矩形分布，

1.3 算法誤差

葉片的評價結果是由坐標測量機測得的點云數據通過葉片參數評價算法給出。就葉片來說，其最終測量結果的誤差源勢必包括了測量設備硬件系統本身誤差和測量軟件處理算法帶入的誤差因素，因此這部分主要研究對軟件評價算法的誤差估計。在算法誤差中，不考慮算法數學模型建立不同的影響，主要來源有方法誤差（截斷誤差），即算法中包含的計算公式（如泰勒公式等）本身是一種求解的近似（連續的離散化處理，無窮的有限話處理）；舍入誤差，即計算機中的數（機器數）是具有有限精度的實數的有限子集，由于計算時的四舍五入，或者因計算機的字長有限而使原始數據只能用有限位數表示，由此產生的誤差。通過合成，算法誤差可以記為0.05μm，服從矩形分布。

1.4 探測誤差

探測誤差是使用坐標測量機測量標準球半徑的示值變化范圍而確定的誤差，主要反映了測頭的各向異性、瞄準誤差和作用直徑的影響，屬于坐標測量機的方向特性參數。探測誤差是影響測量不確定度的重要因素，對于不同的測頭，探測誤差也不同。對葉片測量過程中采用的測頭，通過坐標測量機對標準球進行測量，得到探測誤差為2μm，服從均勻分布。

通過以上分析，論文中對MCM模型是基于以下參數的試驗計算和估計的：三坐標測量機運動誤差、系統的和隨機的三坐標測量機探頭誤差、算法的誤差分析以及環境條件的變化。對于葉片各參數相應的數學模型，為每個輸入量分配適當的概率密度函數。由于蒙特卡羅方法（MCM）是基于對與參考對象的校準值有關的測量誤差的統計評估，因此試驗次數需要足夠多才能正確確定輸出量的分布，因此對于概率p=0.95，確定了M=106的試驗次數。最終得到M個模擬參數值，構成了輸出量的分布，可以確定期望值和標準偏差。

應用蒙特卡洛法對坐標測量機測量葉片型面參數不確定度評定的具體步驟如下：

（1）分析葉片型面參數測量不確定度的來源，對坐標測量機葉片測量過程中的不確定因素進行定量評估，判斷其分布類型及分布區間。通常不確定度來源主要有坐標點測量不確定度、示值誤差引起的不確定度、環境因素等引起的不確定度[15]；

（2）針對葉片型面不同類型的參數，模擬坐標測量機對點云數據進行處理運算，確定誤差模型式中輸入變量的范圍和分布；

（3）考慮到點云數據中每個測量點相互獨立，對點云數據中的每個點，都隨機加入按各個輸入量的分布生成的誤差，來模擬葉片測量中得到的點云數據。通過采用大樣本來進行葉片型面參數測量不確定度評定，樣本容量越大，越能真實地模擬出實際測量值；

（4）根據以上得到的模擬測量數據，代入葉片各參數數學模型中，求得試驗結果值。根據該組值，構造一個概率分布，判斷其分布類型，計算其期望值、標準差和95%置信區間，期望值即為葉片型面參數值，標準差即為所要求的各葉片型面參數測量的標準不確定度，置信區間上下限差值表示葉片型面參數的變化范圍。

2 基于蒙特卡洛法的葉片型面參數測量不確定度評定試驗

為了更全面地表示葉片型面測量過程的不確定度，綜合評價葉片型面的測量結果，選取了葉片型面4類參數中的典型參數，分別為葉型最大厚度、弦長、位置度、扭轉角和輪廓度。

2.1 最大厚度

最大厚度是指葉片截面的最大厚度值[16]，是對葉型輪廓度的補充，為控制葉片的最大厚度分布，通常是對封閉葉型進行內切圓計算，得到內切圓的最大半徑，將最大內切圓的直徑作為葉型的最大厚度值。

葉片測量數據通常是等高截面上離散的點云數據，首先將葉片點云數據粗分割成前緣、后緣、葉盆和葉背4部分，然后對粗分割獲取得到的葉片前后緣數據進一步精確提取。

（1）計算葉片點云數據中距離最遠的兩點，分別為葉片截面的前緣極值點和后緣極值點。

（2）從前后緣極值點分別向葉盆、葉背方向各延伸數點。

（3）對選取的點基于最小二乘法擬合圓，該圓的半徑作為初始前（后）緣半徑。

（4）計算半徑誤差，設待擬合的數據點為（xi，yi），擬合的圓心坐標為（xc，yc），擬合圓半徑為R。

（5）基于初始擬合點往葉盆葉背方向逐步增加數據點，當半徑誤差超過閾值時，到達前后緣與葉盆葉背連接處。

（6）獲取葉盆葉背與前后緣的連接點，實現葉盆葉背和前后緣的分割，并可以擬合出前后緣圓弧對應圓，計算出前后緣圓弧半徑和圓心。

（7）對精確分割得到的葉盆、葉背、前緣和后緣點云數據重新排序、加密來重構葉片截面型線。由于葉盆葉背處點云較稀疏，而且分布疏密不均，利用B樣條曲線插值對葉盆葉背上點進行加密。

（8）求取加密后的葉盆點Ai到葉背Bi之間距離最小值的合集，對葉盆上點Ai，

（9）在最小值合集中找到最大的值，即為求得的最大厚度距離。

對葉盆上所有點A，求得的li對應的合集L，最大厚度h=max{li}。

選取基于三坐標測量機對葉片截面實測的點云數據，對測得的1836個點，每個點隨機加入服從均勻分布空間坐標測量不確定度（1.2+l/350）μm；服從等概率分布的環境因素對測量結果的影響（1.15×10–3lΔt）μm，其中假設Δt=0.3℃；服從均勻分布的校準裝置探觸形狀誤差引入的不確定度2μm。輸入最大厚度的數學模型中模擬106次，得到對應截面最大厚度的期望值，標準差及置信區間為95%的區間范圍，最大厚度直方圖如圖2所示。

由圖2可以看出，數據在[30.358mm，30.372mm]區間內變化，且在中間值區域集中，計算出的數組的期望值作為測量值為30.365mm，標準差作為標準不確定度值為1.76μm，95%置信區間為[30.362mm，30.369mm]，上下限差值計算得到3.43μm。

圖2 蒙特卡洛法最大厚度直方圖Fig.2 Monte Carlo method maximum thickness histogram

2.2 弦長

葉型的弦線是將葉型投影在葉盆的公切線上，葉型弦長即為投影點最小值與最大值之間的距離。弦線定義如圖3所示。

圖3 弦線定義Fig.3 Definition of chord

葉型弦線的計算方法包括前后緣點法、公切線法和最大值法等。試驗中選擇前后緣點法計算弦長，選取基于三坐標測量機對葉片5個不同截面高度的實測的點云數據，對每個截面測得的點云數據，每個點隨機加入服從均勻分布最大示值誤差（1.2+l/350） μm；服從等概率分布的環境因素對測量結果的影響（1.15×10–3lΔt）μm，其中假設Δt=0.3℃；服從均勻分布的測量設備探觸形狀誤差引入的不確定度2μm；服從矩形分布的算法誤差0.05μm。由2.1節中最大厚度數學模型中前后緣點的求取，通過計算前、后緣兩點連線的長度，模擬106次，對得到的多組數據應用蒙特卡洛法評定測量不確定度，結果如表1所示。

通過表1中數據可以看到，對不同高度的不確定度評定，不同高度的測量不確定度差值在1μm以內。

表1 不同高度的弦長不確定度及測量值Table 1 Uncertainty and measurement value of chord length at different heights

最大厚度和弦長參數的評定證明了葉片尺寸類參數的測量不確定度是基本一致的。

2.3 位置度

葉型位置度是葉片實測數據相對于理論葉型坐標點進行最佳擬合時的平移量。葉型在軸向與周向上的平移量分別稱為軸向位置度和周向位置度。

在采用最小二乘法進行葉型擬合時，除了最小化理論數據與檢測數據的方差，也最小化兩者的均方差。在這個過程中，葉型在軸向與周向上的平移量分別稱為軸向位置度（x_offset）和周向位置度（y_offset），葉型在平面內繞積迭點旋轉的角度稱為扭轉角，通常情況下，以檢測截面平面內以逆時針方向為負。

有些情況下，也采用T_offset對位置度進行控制。T_offset的計算如下：

將實測的數據值與理論數據ICP配準后，獲得最優旋轉矩陣R和平移向量T。其中最優旋轉矩陣可以計算得到θ；通過平移向量可以求得位置度參數。葉片位置度和扭轉度評價算法流程圖如圖4所示。

圖4 葉片位置度和扭轉度評價算法流程圖Fig.4 Flow chart of blade position and torsion evaluation algorithm

選取基于三坐標測量機對葉片中間高度的截面實測點云數據，對這1832個點，每個點隨機加入服從均勻最大示值誤差（1.2+l/350） μm；服從等概率分布的環境因素對測量結果的影響（1.15×10–3lΔt）μm，其中假設Δt=0.3℃；服從均勻分布的校準裝置探觸形狀誤差引入的不確定度2μm；服從矩形分布的算法誤差0.05μm。將模擬受到擾動生成的106組數據輸入到位置度數學模型中，通過試驗結果計算這組值的平均值、標準差及95%置信區間，可以得到軸向位置度測量值為– 0.1013mm，標準不確定度為0.089μm；周向位置度測量值為0.0538mm，標準不確定度為0.045μm。結合圖5的直方圖可以看出，不確定度分布上和葉片型面尺寸類參數類似，但由于評定過程中還需要與理論數據進行對齊，因此必須單獨分析。

圖5 軸向和周向位置度直方圖Fig.5 Histograms of axial and circumferential position

2.4 扭轉角

葉型扭轉角是葉片實測數據相對于理論葉型坐標點進行最佳擬合時，葉型在平面內繞積疊點旋轉的角度。

通過對2.3節中計算位置度及扭轉角的數學模型和試驗，對計算得到的模擬參數值進行分析，可得扭轉角測量值為0.0512°，標準不確定度值為0.00031°，其中95%的置信區間對應的不確定度為0.00021°，由圖6扭轉角的直方圖可以看出，葉片型面角度類參數受誤差變化影響較大，導致數據集中偏左。

圖6 扭轉角直方圖Fig.6 Torsion angle histogram

2.5 輪廓度

輪廓度一般指葉片型面上不同部位輪廓上最大值與最小值的差，包括前緣輪廓度、后緣輪廓度、葉盆輪廓度及葉背輪廓度。輪廓度定義如圖7所示。

圖7 輪廓度定義Fig.7 Definition of profile

葉型參數評價最基本的形狀誤差是輪廓度，而輪廓度的評價又需要在曲線曲面測量數據與理論模型的最佳擬合這一前提下實現，因此葉型輪廓度的計算涉及實測數據與理論數據的對齊，有以下計算流程。

（1）對葉型截面的測量數據和截面對應的理論數據進行ICP算法配準。

（2）按照2.1節中最大厚度數學模型中葉盆、葉背、前后緣分割，對葉型測量點集進行分割，得到4個區域的點集。

（3） 對理論葉型型值點進行NURBS曲線插值擬合，形成理論型線。

（4）基于圓裁剪和牛頓迭代算法實現點到NURBS曲線距離計算。

（5）最后根據向量叉積法解決輪廓偏差方向的判定問題。

選取基于三坐標測量機對葉片中間高度的截面實測點云數據，對這1832個點，每個點隨機加入服從均勻分布最大示值誤差（1.2+l/350） μm；服從等概率分布的環境因素對測量結果的影響（1.15×10–3lΔt）μm，其中假設Δt=0.3℃；服從均勻分布的校準裝置探觸形狀誤差引入的不確定度2μm；服從矩形分布的算法誤差0.05μm。基于蒙特卡洛法對測量得到的點云數據進行模擬106次，輸入到輪廓度測量的數學模型中，得到多組輪廓度模擬測量值，通過模擬測量值計算各位置輪廓度的平均值、標準差及95%置信區間，得到測量值和標準不確定度。輪廓度的評定涉及葉片截面區域分割的問題，因此各個位置的輪廓度在不確定度評估中也得到了不同的值，如表2所示。葉背的最小輪廓度分布如圖8所示。

結合表2和圖8可以分析得到前、后緣部分的輪廓度比葉盆、葉背輪廓度不確定度評定結果大，這符合前后緣部分曲率變化大的實際測量情況。葉背輪廓度概率分布與弦長、最大厚度這類僅需要實測數據計算的尺寸類參數的相比，中間區域出現明顯下降和上升，這是由于輪廓度的計算涉及到與理論數據的對齊算法的影響，也從不同角度說明了葉片型面參數的測量不確定度各有不同。

表2 輪廓度不確定度Table 2 Uncertainty of profile

圖8 葉背輪廓度直方圖Fig.8 Leaf profile histogram

3 基于GUM法的葉片弦線參數測量不確定度評定

通過對某一葉片同一高度的截面在相同條件下重復測量多次，應用GUM法對葉片弦長參數的不確定度評定。首先定量分析影響不確定度的各分量。

（1）坐標測量機尺寸最大允許示值誤差引入的測量不確定度u1。

坐標測量機尺寸測量的最大允許示值誤差為（1.2+l/350） μm，按均勻分布，u1=（1.2+l/350）/（μm）。

（2）被測件和設備的熱膨脹系數引入的不確定度分量u2。

被測件和設備的膨脹系數差在半寬為2×10–6℃–1的區間內以等概率分布，Δt=0.3℃ 時，

（3）葉片測量重復性引入的測量不確定度u3。

對葉片在重復條件下連續測量多次，得到的標準偏差：

u3= 2.087（μm）

最后通過分量合成的方法得到葉片型面弦長的測量不確定度為：

按包含概率ρ=95%，則k=1.96，擴展不確定度為U95=1.96×uc=4.686（μm）

4 兩種方法的比較和分析

基于GUM法計算得到的擴展不確定度和基于蒙特卡洛法的葉片弦長測量不確定度評定結果的范圍量級基本一致，同時，基于蒙特卡洛法葉片弦長的分布結果符合概率統計分布規律，充分說明基于蒙特卡洛法的葉片型面參數測量不確定度評定模型的正確性。

基于蒙特卡洛法的葉片弦長測量不確定度評定結果與GUM法計算得到的測量合成標準不確定度相比，蒙特卡洛法可以給出測量不確定度的包含概率，可以更好地表示測量結果的特征。而計算擴展不確定度時，GUM法通常假設模型輸出量屬于正態分布，采用估算的方式，在實際模型不一定是正態分布的情況下，具有一定的局限性；基于蒙特卡洛法評定的不確定度，通過大量模擬測量過程，結果包括了測量值的概率分布圖形，沒有對輸出量分布類型假設所造成的誤差影響，更加的合理、準確。

5 結論

本文提出了一種基于蒙特卡洛法評定坐標測量機測量葉片型面參數不確定度的計算方法。該方法充分利用葉片實測點云數據和理論模型，基于葉片型面參數評定算法及蒙特卡洛模擬試驗，實現對葉片型面參數不確定度的準確評估。基于該方法進行的葉片型面參數測量不確定度評定結果與傳統GUM法評定處于一個量級范圍；利用該方法初步掌握了葉片尺寸類、位置類、角度類和輪廓類典型型面參數的測量不確定度評定規律。同傳統GUM評定方法相比，基于蒙特卡洛模擬試驗統計計算不確定度的方法在理論上具有優勢性、可靠性和在線性，值得在葉片型面參數測量及評價領域進行推廣應用。