吳呼玲

(陜西國(guó)防工業(yè)職業(yè)技術(shù)學(xué)院 機(jī)械工程學(xué)院，西安 710300)

基于蒙特卡羅法的平面度測(cè)量不確定度評(píng)定

吳呼玲

(陜西國(guó)防工業(yè)職業(yè)技術(shù)學(xué)院 機(jī)械工程學(xué)院，西安 710300)

形位誤差的測(cè)量不確定度評(píng)定是目前測(cè)量領(lǐng)域研究的熱點(diǎn)；但由于其測(cè)量的復(fù)雜性和測(cè)量結(jié)果評(píng)定的多樣性，導(dǎo)致在實(shí)際測(cè)量結(jié)果中形位誤差測(cè)量的不確定度評(píng)定成了難題；為此，根據(jù)形狀誤差評(píng)定準(zhǔn)則，選取最小二乘法建立數(shù)學(xué)模型，確定形狀誤差數(shù)學(xué)模型中各參數(shù)值的傳遞系數(shù)和單點(diǎn)不確定度，并分析具體的測(cè)量方法和測(cè)量過(guò)程中的不確定度來(lái)源，根據(jù)傳統(tǒng)的GUM法對(duì)其進(jìn)行不確定度評(píng)定；然后采用蒙特卡羅偽隨機(jī)數(shù)的方法來(lái)模擬實(shí)際測(cè)量數(shù)據(jù)，從而得到平面度誤差的不確定度；通過(guò)設(shè)置實(shí)驗(yàn)對(duì)比，驗(yàn)證了蒙特卡羅法評(píng)定平面度不確定度的可靠性和準(zhǔn)確性；該方法不需要求出數(shù)學(xué)模型中的傳遞系數(shù)，利用MATLAB軟件很容易實(shí)現(xiàn)，為平面度誤差測(cè)量結(jié)果不確定度評(píng)定提供了更加簡(jiǎn)便的方法，值得推廣和應(yīng)用。

蒙特卡羅；平面度誤差；不確定度；最小二乘法

0 引言

平面度測(cè)量誤差的評(píng)定是測(cè)量領(lǐng)域中經(jīng)常遇到的誤差評(píng)定項(xiàng)目。因此，準(zhǔn)確、快速的評(píng)定平面度測(cè)量結(jié)果的不確定度成為一個(gè)重要的課題。隨著測(cè)量不確定度在我國(guó)受到越來(lái)越高的重視，對(duì)于從事測(cè)量的專業(yè)技術(shù)人員來(lái)講，測(cè)量不確定度成為了測(cè)量領(lǐng)域研究的熱點(diǎn)[9]。然而其復(fù)雜性和多樣性的特點(diǎn)，致使不確定度的評(píng)定難以評(píng)估。

蒙特卡羅法[1]又稱統(tǒng)計(jì)模擬法、隨機(jī)抽樣技術(shù)，是使用隨機(jī)數(shù)(或偽隨機(jī)數(shù))來(lái)解決問題的一種方法。它的基本思想是，為了所求問題的解是某個(gè)隨機(jī)變量的數(shù)學(xué)期望時(shí)，通過(guò)某種實(shí)驗(yàn)的方法，得出該隨機(jī)變量若干個(gè)具體觀察值的算術(shù)平均值，通過(guò)它得到問題的解。本文利用蒙特卡羅法進(jìn)行平面度測(cè)量結(jié)果不確定度的評(píng)定，利用MATLAB軟件產(chǎn)生一組服從測(cè)量值分布概率的隨機(jī)數(shù)組來(lái)模擬測(cè)量值，從而利用統(tǒng)計(jì)方法來(lái)進(jìn)行不確定度評(píng)定。為了驗(yàn)證此方法的準(zhǔn)確性和可靠性，根據(jù)不確定度合成公式再次進(jìn)行評(píng)定，求出傳遞系數(shù)和各參數(shù)的單點(diǎn)測(cè)量不確定度，計(jì)算出平面度誤差的測(cè)量不確定度。實(shí)驗(yàn)結(jié)果表明，蒙特卡羅法評(píng)定平面度測(cè)量不確定度準(zhǔn)確、快捷、簡(jiǎn)單、可靠，解決了形位誤差測(cè)量不確定度評(píng)定的復(fù)雜性和多樣性，為其他形位誤差簡(jiǎn)單、便捷的評(píng)定提供了評(píng)定方案。

1 平面度誤差模型

本文用最小二乘法原理建立數(shù)學(xué)模型。

(1)

(2)

(3)

(4)

2 平面度測(cè)量結(jié)果不確定度評(píng)定

2.1 蒙特卡羅法評(píng)定平面度測(cè)量結(jié)果不確定度

根據(jù)平面度最小二乘法評(píng)定的誤差數(shù)學(xué)模型，利用蒙特卡羅偽隨機(jī)數(shù)的方法產(chǎn)生服從期望為各參數(shù)的測(cè)量值，方差為各參數(shù)的單點(diǎn)標(biāo)準(zhǔn)不確定度的正太分布的隨機(jī)值序列，得出平面度誤差測(cè)量不確定度。

具體操作的方法和步驟如下[7]：

1)分析測(cè)量過(guò)程中不確定度來(lái)源。確定分布區(qū)間和分布類型，確定因子K的取值，計(jì)算出各個(gè)方面不確定度來(lái)源對(duì)應(yīng)的不確定度數(shù)值。為后續(xù)數(shù)組中的方差提供依據(jù)(方差為各參數(shù)的不確定度)。

2)確定平面度誤差模型(公式4)中的參數(shù)：xM；xL；yM；yL；zM；zL；a；b的期望和方差(期望為各參數(shù)的測(cè)量值，方差為各參數(shù)的單點(diǎn)標(biāo)準(zhǔn)不確定的度)。

3)以xM；xL；yM；yL；zM；zL；a；b參數(shù)的期望和方差生成八維隨機(jī)數(shù)來(lái)模擬平面度誤差的測(cè)量值，樣本容量為M，采用大樣本進(jìn)行平面度誤差的測(cè)量不確定度評(píng)定。生成的八維隨機(jī)數(shù)分別為：xM1，xM2，xM3，……xMM；xL1，xL2，xL3，……xLM；yM1，yM2，yM3，……yMM；yL1，yL2，yL3，……yLM；zM1，zM2，zM3，……zMM；zL1，zL2，zL3，……zLM；a1，a2，a3，……aM；b1，b2，b3，……bM；

4)根據(jù)以上隨機(jī)序列，帶入(公式4)中平面度誤差模型，求出M個(gè)平面度誤差f的值，根據(jù)這組f的值，構(gòu)造一個(gè)概率分布。判斷分布類型求出方差，即為所要求的平面度誤差測(cè)量標(biāo)準(zhǔn)不確定度。

2.2GUM法評(píng)定平面度測(cè)量結(jié)果不確定度

1)求出(4)式中平面度誤差模型各參數(shù)的傳遞系數(shù)[3]:

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

2)各參數(shù)測(cè)量不確定度分別為：

平面度測(cè)量不確定度的分析，只需考慮分析z的不確定度即可[4]。z值的測(cè)量不確定度等于單點(diǎn)測(cè)量不確定度，即u(zM)=u(zL)=u0。而a和b的不確定度比較復(fù)雜，需要進(jìn)行以下推導(dǎo)得到。

因此，

u(zM)=u(zL)=u0

(13)

(14)

(15)

將以上平面度誤差模型各參數(shù)的傳遞系數(shù)(公式9～12)和各參數(shù)單點(diǎn)測(cè)量不確定度(公式13～15)帶入不確定度合成公式16中，即可得到平面度誤差測(cè)量不確定度評(píng)定結(jié)果：

uf=

(16)

3 實(shí)驗(yàn)數(shù)據(jù)分析及應(yīng)用

使用愛德華公司的MQ686型三坐標(biāo)測(cè)量機(jī)對(duì)一塊平板進(jìn)行平面度測(cè)量和評(píng)定。將被測(cè)平板放置于三坐標(biāo)測(cè)量機(jī)的測(cè)量平臺(tái)上，并使被測(cè)工件的長(zhǎng)邊平行于X軸，然后對(duì)平面度誤差進(jìn)行測(cè)量。被測(cè)工件長(zhǎng)度80 mm，在X軸方向均勻取6個(gè)點(diǎn)，每個(gè)點(diǎn)間隔10 mm進(jìn)行坐標(biāo)測(cè)量，共測(cè)3行數(shù)據(jù)，即在被測(cè)平面上取3行6列共18個(gè)數(shù)據(jù)，布點(diǎn)情況[5]如圖1所示，測(cè)量結(jié)果如表1所示。

圖1 三坐標(biāo)測(cè)量機(jī)測(cè)量平板平面度的布點(diǎn)情況

測(cè)點(diǎn)序號(hào)X坐標(biāo)/mmY坐標(biāo)/mmZ坐標(biāo)/mm19．999710．0019-0．0079219．99979．9964-0．0066330．00299．9962-0．0082439．99989．9945-0．0064549．99809．9940-0．0060659．99879．9946-0．005779．986219．9999-0．0064819．995719．9946-0．0066930．002919．9960-0．00571039．999619．9954-0．00591149．998419．9954-0．00451259．999019．9960-0．0051139．985229．9992-0．00741419．995229．9944-0．00761530．002429．9954-0．00671639．999029．9957-0．00641749．998829．9953-0．00651859．999029．9958-0．0061

3.1 不確定度來(lái)源分析[8]

本次實(shí)驗(yàn)采用的MQ686型全封閉框架移動(dòng)橋式三坐標(biāo)測(cè)量機(jī)，其量程范圍是:X600 mm，Y800 mm，Z600 mm，該儀器處于良好的工作狀態(tài)，室內(nèi)空調(diào)已經(jīng)連續(xù)運(yùn)行一周，使被測(cè)平板和三坐標(biāo)等溫，測(cè)量人員在1小時(shí)后進(jìn)行測(cè)量工作。因此，可以忽略溫度、濕度等環(huán)境帶來(lái)的影響。其測(cè)量不確定度主要來(lái)自以下幾個(gè)方面：

1)重復(fù)性引起的不確定度分量[10]。

三坐標(biāo)測(cè)量機(jī)的重復(fù)性測(cè)量引起的不確定度，在平面度的測(cè)量不確定度的評(píng)定過(guò)程中是一個(gè)重要的因素。在平面度測(cè)量過(guò)程中，平板的長(zhǎng)邊是平行于X軸放置的，短邊是平行于Y軸的，所以影響平面度誤差的主要因素是Z軸的坐標(biāo)值。因此，在評(píng)定重復(fù)性引入的不確定度過(guò)程中，只需要考慮Z軸坐標(biāo)測(cè)量的重復(fù)性誤差。

在測(cè)量的18個(gè)點(diǎn)中，選取3個(gè)點(diǎn)(1點(diǎn)、6點(diǎn)和15點(diǎn))進(jìn)行重復(fù)性測(cè)量實(shí)驗(yàn)，利用貝塞爾公式求出各點(diǎn)的重復(fù)性誤差，取這三個(gè)點(diǎn)中重復(fù)性測(cè)量引起的不確定度的最大值，作為三坐標(biāo)測(cè)量機(jī)的重復(fù)性引起的不確定度。

選取1點(diǎn)重復(fù)測(cè)量10次的測(cè)量數(shù)據(jù)為：

-0.0082；-0.0079；-0.0081；-0.0083；-0.0081，

-0.0078；-0.0076；-0.0082；-0.0079；-0.0075

u6=σ6=0.074 μm；u15=σ15=0.067 μm；

取最大值作為重復(fù)性所引起的不確定度，因此三坐標(biāo)測(cè)量機(jī)的重復(fù)性測(cè)量引入的不確定度為u重復(fù)=0.085μm。

2)示值誤差引入的不確定度分量。

3)測(cè)量力引入的不確定度分量。

在測(cè)量中測(cè)量力調(diào)節(jié)至很小，可忽略測(cè)量力引起的零件彈性變形。因此，測(cè)量力的影響帶來(lái)的不確定度分量為u力=0。

4)三坐標(biāo)測(cè)量機(jī)分辨力引入的不確定度分量。

5)溫度引入的不確定度分量。

在測(cè)量過(guò)程中，實(shí)驗(yàn)室空調(diào)已經(jīng)連續(xù)運(yùn)行了一周，被測(cè)件、三坐標(biāo)和測(cè)量人員都在等溫后進(jìn)行測(cè)量工作。因此，測(cè)量力的影響帶來(lái)的不確定度分量為u溫度=0。

根據(jù)以上分析，直線度的單點(diǎn)測(cè)量不確定度為：

3.2 平面度誤差計(jì)算

依據(jù)表1中的測(cè)量數(shù)據(jù)，帶入公式(1～3)，借助MATLAB軟件求出a=3.551696623161069e-05、b=8.394157412702668e-07、c=-7.687578147000283e-03，可得最小二乘平面為z=(3.55e-05)x+(8.39e-07)y-(7.69e-03)。然后求出各測(cè)點(diǎn)到最小二乘平面的距離，計(jì)算結(jié)果如表2所示。

表2 各測(cè)點(diǎn)到最小二乘平面的距離

由表2可知，距離最小二乘平面最高的點(diǎn)是測(cè)點(diǎn)11(在此稱為最大距離的點(diǎn))，對(duì)應(yīng)的點(diǎn)坐標(biāo)為(xM，yM,zM)、最低的點(diǎn)是測(cè)點(diǎn)3(在此稱為最小距離的點(diǎn))，對(duì)應(yīng)的點(diǎn)坐標(biāo)為(xL，yL,zL)。將各參數(shù)帶入平面度誤差模型公式(4)中，則平面度誤差值為：

2.981427013954202e-03mm=2.98μm

3.3 蒙特卡羅法的不確定度估算

利用MATLAB軟件中的normrnd函數(shù)[6]生成服從正態(tài)分布的八維隨機(jī)數(shù),每一維數(shù)組的期望分別為各參數(shù)測(cè)量值，方差分別為各參數(shù)的不確定度，并且服從正態(tài)分布，即為：

[xM1，xM2，xM3，……xMM] ～N(49.9984mm，1.56μm)；

[xL1，xL2，xL3，……xLM] ～N(30.0029mm，1.56μm)

[yM1，yM2，yM3，……yMM] ～N(19.9954mm，1.56μm)

[yL1，yL2，yL3，……yLM] ～N(9.9962mm，1.56μm)

[zM1，zM2，zM3，……zMM] ～N(-0.0045mm，1.56μm)

[zL1，zL2，zL3，……zLM] ～N(-0.0082mm，1.56μm)

[a1，a2，a3，……aM] ～N(3.551696623161069e-05 mm，1.767929943168856e-04)

[b1，b2，b3，……bM] ～N(8.394157412702668e-07mm，4.420132407275619e-04)

樣本容量為M=10 000。由于樣本很大，這里不予列出。根據(jù)所得隨機(jī)數(shù)，求出10 000個(gè)f的值后，構(gòu)建出平面度誤差概率分布，其概率分布根據(jù)MATLAB中的直方圖(橫坐標(biāo)為平面度誤差值，縱坐標(biāo)為其分布概率)來(lái)構(gòu)建，求出這組數(shù)據(jù)的期望和方差，即為平面度誤差和平面度測(cè)量不確定度，分別求得：

平面度誤差：f=2.966479016299404e-03 mm=2.97 μm；

平面度測(cè)量不確定度：uf=6.036951149940499e-03 mm=6.04 μm；

圖2 平面度誤差的概率分布

3.4 GUM法的不確定度驗(yàn)證

q=3.599307000000000e+02；u=2.729853448666000e+04；

v= 8.397192341290000e+03；w=1.259620107763000e+04；

S=1.134251187913089e+08根據(jù)公式(16)計(jì)算出平面度測(cè)量不確定度結(jié)果：

uf=6.074399101030872e-03 mm=6.07 μm

3.5 結(jié)果比較

表3 兩種方法計(jì)算的平面度誤差和評(píng)定的測(cè)量不確定度結(jié)果

4 結(jié)語(yǔ)

通過(guò)表3數(shù)據(jù)比較可知，用蒙特卡羅法進(jìn)行平面度不確定度評(píng)定是可靠的。通過(guò)GUM法驗(yàn)證，兩者結(jié)果幾乎相同，而且蒙特卡羅法比GUM法評(píng)定結(jié)果更加準(zhǔn)確。因此，通過(guò)實(shí)驗(yàn)驗(yàn)證，蒙特卡羅法評(píng)定平面度不確定度以評(píng)定結(jié)果精度高、計(jì)算快捷，簡(jiǎn)便、不必計(jì)算偏導(dǎo)數(shù)、不受模型復(fù)雜性的影響而凸顯出它的優(yōu)勢(shì)。為以后其他形位誤差測(cè)量不確定度的評(píng)定和其他數(shù)據(jù)處理領(lǐng)域的應(yīng)用奠定了一定的基礎(chǔ)，具有非常廣泛的應(yīng)用前景和工程實(shí)用價(jià)值。

[1] 周桃庚.用蒙特卡洛法評(píng)定測(cè)量不確定度[M].北京:中國(guó)質(zhì)檢出版社 2013.

[2] 倪驍驊.形狀誤差評(píng)定和測(cè)量不確定度估計(jì)[M].北京:化學(xué)工業(yè)出版社，2008.

[3] JJF1059.1-2012《測(cè)量不確定度的評(píng)定與表示》[S].

[4] 倪驍驊，黃曉峰，葛友華 平面度誤差最小區(qū)域評(píng)定結(jié)果不確定度估計(jì)[J].鹽城工學(xué)院學(xué)報(bào)(自然科學(xué)版),2006,19(2)：5-8.

[5] 王金星，蔣向前，馬利民，等.平面度坐標(biāo)測(cè)量的不確定度計(jì)算[J].中國(guó)機(jī)械工程,2005,16(19)：1701-1703.

[6] 王 偉，宋明順、陳意華. 蒙特卡羅方法在復(fù)雜模型測(cè)量不確定度評(píng)定中的應(yīng)用[J].儀器儀表學(xué)報(bào),2008,29(7)：1446-1449.

[7] 連慧芳，陳曉懷.基于蒙特卡羅方法的圓度測(cè)量不確定度評(píng)定[J]. 工具技術(shù),2010,44(6)：82-84.

[8] 王東霞，宋愛國(guó).基于三坐標(biāo)測(cè)量機(jī)的圓度誤差不確定度評(píng)估[J].東南大學(xué)學(xué)報(bào)(自然科學(xué)版),2014，44(5)：952-956.

[9] 陳懷艷，曹蕓，韓潔.測(cè)量不確定度的發(fā)展和應(yīng)用研究[J].宇航計(jì)測(cè)技術(shù), 2014,10,34(5)：65-70.

[10] 費(fèi)業(yè)泰.誤差理論與數(shù)據(jù)處理[M].北京:機(jī)械工業(yè)出版社, 2012.

Uncertainty Evaluation of Flatness Measurement Based on Monte Carlo Method

Wu Huling

(School of Mechanical Engineering, Shaanxi Institute of Technology, Xi’an 710300, China)

Evaluation of measurement uncertainty of shape and position error is a hotspot in the field of measurement. However, due to the complexity of measurement and the diversity of measurement results, the uncertainty of the measurement of the shape and position error in the actual measurement results has become a difficult problem. Therefore, according to the evaluation criterion of shape error, establish mathematical model adopting least squares method, determining the transfer coefficient of each parameter value of the shape error mathematical model and single point uncertainty, and uncertainty analysis of measurement methods and the specific process of the source, according to the traditional GUM method for uncertainty evaluation. Then the Monte Carlo pseudo random number method is used to simulate the actual measurement data, thus the flatness error uncertainty is obtained. The reliability and accuracy of Monte Carlo method to evaluate flatness uncertainty are verified by setting experimental comparison. This method does not need to calculate the transfer coefficient of mathematical model, it is easy to realize by using Matlab software, and provides a more convenient method for the uncertainty evaluation of flatness error measurement results.

Monte Carlo method；flatness error；uncertainty；least square method

2016-12-01；

2017-01-05。

陜西國(guó)防工業(yè)職業(yè)技術(shù)學(xué)院2016年科研項(xiàng)目(Gfy16-03)。

吳呼玲(1979-)，女，陜西臨潼人，講師，碩士，主要從事機(jī)械產(chǎn)品檢驗(yàn)檢測(cè)、誤差理論與數(shù)據(jù)處理、機(jī)械設(shè)備狀態(tài)監(jiān)測(cè)等方向的研究。

1671-4598(2017)05-0262-04

10.16526/j.cnki.11-4762/tp.2017.05.072

TH124

A