基于蒙特卡洛方法的支撐劑導(dǎo)流能力不確定度評(píng)價(jià)

陶祖文

（中國(guó)石化西南石油工程有限公司井下作業(yè)分公司，四川德陽(yáng) 618000）

基于蒙特卡洛方法的支撐劑導(dǎo)流能力不確定度評(píng)價(jià)

陶祖文

（中國(guó)石化西南石油工程有限公司井下作業(yè)分公司，四川德陽(yáng) 618000）

為解決壓裂支撐劑充填層導(dǎo)流能力測(cè)量不確定度評(píng)定時(shí)數(shù)學(xué)模型過(guò)于復(fù)雜以及部分偏導(dǎo)數(shù)難以求解的問(wèn)題，本論文采用蒙特卡洛方法對(duì)支撐劑充填層導(dǎo)流能力不確定度進(jìn)行評(píng)價(jià)，推導(dǎo)了概率密度函數(shù)對(duì)應(yīng)的隨機(jī)變量的獲取公式，研究表明試驗(yàn)次數(shù)在10萬(wàn)～100萬(wàn)次之間，采用蒙特卡洛方法計(jì)算出的估計(jì)值和標(biāo)準(zhǔn)不確定度結(jié)果穩(wěn)定，并且試驗(yàn)次數(shù)越大，精確度越高。

蒙特卡洛法；壓裂；支撐劑；導(dǎo)流能力；測(cè)量不確定度

測(cè)量不確定度是指由于測(cè)量誤差的存在，對(duì)被測(cè)量值不能肯定的程度，是表征測(cè)量質(zhì)量的重要指標(biāo)[1，2]。根據(jù)CNAS-CL07關(guān)于測(cè)量不確定度的要求[3]，在客戶有要求或不確定度影響到對(duì)結(jié)果符合性的判定時(shí)，校準(zhǔn)實(shí)驗(yàn)室必須對(duì)校準(zhǔn)數(shù)據(jù)給出不確定度及評(píng)定程序。目前，測(cè)量不確定度已在油氣計(jì)量領(lǐng)域得到廣泛的應(yīng)用[4]，但是受到各不確定度分量大小不相近、測(cè)量模型的偏導(dǎo)數(shù)難以求得、輸入量的概率密度函數(shù)不對(duì)稱等方面的限制。

但是，蒙特卡洛法為測(cè)量不確定度評(píng)定提供了一個(gè)通用的數(shù)值方法，適用于任意多個(gè)由概率密度函數(shù)表征的輸入量和單一輸出量的模型[5]。蒙特卡洛方法又稱隨機(jī)模擬方法，它以概率理論為基礎(chǔ)理論，以隨機(jī)抽樣為主要手段[6]。國(guó)外學(xué)者Willink R[7]、Esward T J等[8]研究了蒙特卡洛方法在不確定度評(píng)定中的應(yīng)用，Crowder S V等[9]、Hall B D[10]研究了二階蒙特卡洛方法在非線性方程和小樣本量測(cè)試中的不確定度評(píng)定中的應(yīng)用，F(xiàn)ernández M S等[11]研究了自適應(yīng)蒙特卡洛方法的不確定度評(píng)定方法。國(guó)內(nèi)該領(lǐng)域研究方面，祖先鋒等[12]、王珊等[13]、高玉英[14]以蒙特卡洛方法為基礎(chǔ)，在現(xiàn)代測(cè)試及儀器精度領(lǐng)域進(jìn)行了不確定度評(píng)定。本文詳細(xì)闡述了蒙特卡洛方法評(píng)定步驟及隨機(jī)數(shù)的獲取方法，并以油氣計(jì)量領(lǐng)域?qū)嵗M(jìn)行說(shuō)明。

1 蒙特卡洛方法及不確定度評(píng)定步驟

蒙特卡洛方法的基本思想是針對(duì)所要求解的問(wèn)題，對(duì)各個(gè)輸入量的概率密度函數(shù)進(jìn)行依次輸入，抽樣時(shí)調(diào)用相關(guān)函數(shù)產(chǎn)生服從相應(yīng)概率分布的隨機(jī)數(shù)，結(jié)合建立的數(shù)學(xué)模型及最大信息熵原理進(jìn)行數(shù)值仿真，從而獲得數(shù)學(xué)模型的方差、期望值、包含區(qū)間等參數(shù)。蒙特卡洛方法不確定度評(píng)定流程圖（見(jiàn)圖1）。

圖1 蒙特卡洛方法不確定度評(píng)定步驟

1.1 建立數(shù)學(xué)模型

當(dāng)輸出量 Y 是由相應(yīng)的輸入量 X=（X1，X2，…，XN）確定時(shí)，建立的數(shù)學(xué)模型即為：

式中：X=（X1，X2，…，XN）。

1.2 輸入量的確定

輸入量的確定主要是基于貝葉斯定理及最大熵原理，根據(jù)已有的實(shí)驗(yàn)數(shù)據(jù)、計(jì)量理論及經(jīng)驗(yàn)，為每個(gè)獨(dú)立的輸入量確定概率密度函數(shù)，常見(jiàn)的統(tǒng)計(jì)分布（見(jiàn)表1）。

表1 常見(jiàn)的輸入量信息及其對(duì)應(yīng)的概率密度函數(shù)

1.3 試驗(yàn)次數(shù)M的確定

蒙特卡洛方法產(chǎn)生的誤差ε[15]：

式中：ε-誤差；λa-置信度a對(duì)應(yīng)的值；σ-均方差；M-數(shù)值計(jì)算次數(shù)。

由公式（2）可知，試驗(yàn)次數(shù)M對(duì)于蒙特卡洛方法測(cè)量不確定度具有重要作用，M值越大，通過(guò)仿真模擬產(chǎn)生的誤差越小，M值越小，產(chǎn)生的誤差就會(huì)增大，輸出值就會(huì)失真。

當(dāng)測(cè)量結(jié)果的不確定度取小數(shù)點(diǎn)后兩位有效數(shù)字時(shí)：

其中，自由度的公式為：

式中：v-自由度；μ（y）-標(biāo)準(zhǔn)不確定度。

可知M≥2×104，所以M取值105～106即可滿足要求。

1.4 輸入量隨機(jī)數(shù)的獲取方法

考慮到統(tǒng)計(jì)分布函數(shù)與獲取的真實(shí)值之間的一致性，實(shí)驗(yàn)測(cè)量數(shù)據(jù)預(yù)測(cè)宜采用矩形分布、正態(tài)分布及三角分布等。這里以三角分布為例，分析如何產(chǎn)生服從三角分布的隨機(jī)數(shù)。通過(guò)計(jì)算機(jī)產(chǎn)生均勻分布的偽隨機(jī)數(shù)ur～[0，1]，三角分布的概率密度函數(shù)為：

其中：a'-三角形分布的底限，c'-眾數(shù)（最可能值），b'-上限，其分布函數(shù)為：

反函數(shù)變換為：

從而得到需要的隨機(jī)數(shù)：

其他分布形式的概率密度函數(shù)以及隨機(jī)變量產(chǎn)生方法（見(jiàn)表 2）。

表2 概率密度函數(shù)對(duì)應(yīng)的隨機(jī)變量

2 支撐劑充填層短期導(dǎo)流能力不確定度評(píng)定實(shí)例

2.1 測(cè)量流程及測(cè)量模型

圖2 測(cè)量流程圖

表3支撐劑充填層導(dǎo)流能力數(shù)據(jù)表

壓裂支撐劑充填層短期導(dǎo)流能力的測(cè)量以支撐裂縫導(dǎo)流儀為基礎(chǔ)，通過(guò)程序采集得到巖心夾持器兩端的壓差以及液體流量，即可通過(guò)公式求得支撐劑充填層導(dǎo)流能力，測(cè)量流程圖（見(jiàn)圖2）。

具體的測(cè)量模型為：

式中：KWf-支撐劑充填層的導(dǎo)流能力，μm2·cm；Q-流量，cm3/min；μ-黏度，mPa·s；ΔP-壓差，kPa；t-溫度，°C。

2.2 測(cè)量數(shù)據(jù)

根據(jù)SY/T 6302-2009《壓裂支撐劑充填層短期導(dǎo)流能力評(píng)價(jià)推薦方法》，環(huán)境溫度為24℃，在蒸餾水或去離子水溫度為24℃時(shí)，其黏度為0.911MPa·s。在閉合壓力為6.89MPa的條件下進(jìn)行8次獨(dú)立測(cè)量，觀察記錄巖心夾持器兩端的壓差及液體流量，具體實(shí)驗(yàn)數(shù)據(jù)（見(jiàn)表 3）。

2.3 采用蒙特卡洛方法進(jìn)行不確定度評(píng)定

2.3.1 建立數(shù)學(xué)模型 支撐劑充填層導(dǎo)流能力測(cè)量的數(shù)學(xué)模型為式（9）～式（11）。

2.3.2 輸入量的確定

（1）壓差的確定：壓力變送器的最大允許誤差為±0.5%，在閉合壓力6.89MPa時(shí)，壓差的平均值為0.667 kPa，下限 a為 0.664 kPa，上限 b 為 0.670 kPa。根據(jù)最大熵原理，壓差設(shè)定為矩形分布R（0.664，0.670）。

（2）流量的確定：查閱計(jì)量泵的校準(zhǔn)證書可知，流量為10 mL/min時(shí)的擴(kuò)展不確定度為1.2%，包含因子k=2，則標(biāo)準(zhǔn)不確定度為0.6%。因此，流量的最佳估計(jì)值為10 mL/min，根據(jù)最大熵原理，流量設(shè)定為正態(tài)分布 N（10，0.0062）。

（3）黏度的確定：由公式（10）可知，黏度的確定與溫度直接相關(guān)。查閱溫度傳感器的校準(zhǔn)證書可知，溫度為24℃時(shí)的擴(kuò)展不確定度為0.13℃，包含因子k=2，則標(biāo)準(zhǔn)不確定度為0.065℃。因此，溫度的最佳估計(jì)值為24℃，根據(jù)最大熵原理，溫度設(shè)定為正態(tài)分布N（24，0.0652），黏度由溫度隨機(jī)數(shù)獲得。

2.3.3 MATLAB編程及總結(jié)報(bào)告 以蒙特卡洛方法為基礎(chǔ)，利用MATLAB軟件進(jìn)行計(jì)算（代碼見(jiàn)附件1），試驗(yàn)次數(shù)分別取10萬(wàn)、50萬(wàn)、100萬(wàn)次，輸出量的估計(jì)值及標(biāo)準(zhǔn)不確定度分別按公式（12）、公式（13）求?。?/p>

式中：yr-不同輸入量對(duì)應(yīng)的支撐劑充填層導(dǎo)流能力，μm·2cm；不同支撐劑充填層導(dǎo)流能力的平均值，μm2·cm；μ（y）-支撐劑充填層導(dǎo)流能力的標(biāo)準(zhǔn)不確定度，μm2·cm。

計(jì)算結(jié)果（見(jiàn)表4、圖3～圖5）?？梢钥闯觯囼?yàn)次數(shù)在10～100萬(wàn)次時(shí)，估計(jì)值及不確定度保留小數(shù)點(diǎn)后兩位小數(shù)時(shí)，導(dǎo)流能力估計(jì)值均為75.88 μm2·cm，標(biāo)準(zhǔn)不確定度均為0.23 μm2·cm，蒙特卡洛方法計(jì)算結(jié)果穩(wěn)定。同時(shí)，試驗(yàn)次數(shù)越大，概率對(duì)稱的95%的包含區(qū)間越小，概率密度曲線越集中，精確度越高。

3 結(jié)論

（1）用蒙特卡洛方法進(jìn)行不確定度評(píng)定適用于具有單一輸出量、輸入量能夠近似為概率密度函數(shù)的任意模型，特別是測(cè)量公式過(guò)于復(fù)雜或偏導(dǎo)數(shù)難以確定的數(shù)學(xué)模型。

表4 支撐劑充填層導(dǎo)流能力數(shù)值計(jì)算結(jié)果

圖3 試驗(yàn)次數(shù)為10萬(wàn)次的導(dǎo)流能力計(jì)算結(jié)果

圖4 試驗(yàn)次數(shù)為50萬(wàn)次的導(dǎo)流能力計(jì)算結(jié)果

圖5 試驗(yàn)次數(shù)為100萬(wàn)次的導(dǎo)流能力計(jì)算結(jié)果

（2）分析了蒙特卡洛方法實(shí)現(xiàn)不確定度評(píng)定的實(shí)施步驟，并推導(dǎo)了概率密度函數(shù)對(duì)應(yīng)的隨機(jī)變量的獲取公式。

（3）支撐劑充填層短期導(dǎo)流能力測(cè)量不確定度評(píng)定結(jié)果表明，試驗(yàn)次數(shù)在10萬(wàn)～100萬(wàn)次，采用蒙特卡洛方法計(jì)算出的估計(jì)值和標(biāo)準(zhǔn)不確定度結(jié)果穩(wěn)定，試驗(yàn)次數(shù)越大，概率對(duì)稱的95%的包含區(qū)間越小，精確度越高。

［1]黃美發(fā)，景暉，匡兵，等.基于擬蒙特卡羅方法的測(cè)量不確定度評(píng)定［J］.儀器儀表學(xué)報(bào)，2009，（1）：120-125.

［2]劉建坤，朱家平，鄭榮華.測(cè)量不確定度評(píng)定研究現(xiàn)狀及進(jìn)展［J］.現(xiàn)代科學(xué)儀器，2013，10（5）：12.

［3]CNAS-CL07.測(cè)量不確定度的要求［S］.中國(guó)合格評(píng)定國(guó)家認(rèn)可委員會(huì)，2011.

［4]邱鳴霞.測(cè)量不確定度的分析及應(yīng)用［J］.黑龍江科技信息，2008，（26）：10.

［5]JJF1059.2-2012.用蒙特卡洛法評(píng)定測(cè)量不確定度［S］.國(guó)家質(zhì)量監(jiān)督檢驗(yàn)檢疫總局，2012.

［6]陳懷艷，曹蕓，韓潔，等.基于蒙特卡羅法的測(cè)量不確定度評(píng)定［J］.電子測(cè)量與儀器學(xué)報(bào)，2011，25（4）：301-308.

［7]Willink R.On using the Monte Carlo method to calculate uncertainty intervals［J］.Metrologia，2006，43（6）：L39.

［8]Esward T J，de Ginestous A，Harris P M，et al.A Monte Carlo method for uncertainty evaluation implemented on a distributed computing system［J］.Metrologia，2007，44（5）：319.

［9]Crowder S V，Moyer R D.A two-stage Monte Carlo approach to the expression of uncertainty with non-linear measurement equation and small sample size［J］.Metrologia，2005，43（1）：34.

［10]Hall B D.Monte Carlo uncertainty calculations with smallsample estimates of complex quantities［J］.Metrologia，2006，43（3）：220.

［11]Fernández M S，Calderón J A，Díez P B.Implementation in MATLAB of the adaptive Monte Carlo method for the evaluation of measurement uncertainties［J］.Accreditation and Quality Assurance，2009，14（2）：95-106.

［12]祖先鋒，毛健人，蔣志文.基于自助法仿真的不確定度抽樣反分析方法研究［J］.計(jì)算機(jī)應(yīng)用與軟件，2011，28（6）：215-217.

［13]Wang S，Chen X，Yang Q.Evaluation of measurement uncertainty based on Bayesian information fusion:Sixth International Symposium on Precision Mechanical Measurements［C］.International Society for Optics and Photonics，2013.

［14]高玉英.基于貝葉斯理論的動(dòng)態(tài)不確定度評(píng)定方法的研究［D］.合肥：合肥工業(yè)大學(xué)，2007.

［15]Li J H，Tong R S，Cao S S.Analysis of Uncertainty in Standard Metal Tank Volume Verification with Monte Carlo Method:Advanced Materials Research［C］.Trans Tech Publ，2013.

Evaluation of uncertainty in measurement based on Monte-Carlo method

TAO Zuwen
（Downhole Operation Branch of Sinopec Southwest Petroleum Engineering Company，Deyang Sichuan 618000，China）

When the measuring mathematical model and the partial derivatives are too difficult to get,Monte-Carlo method is applicable to a single output and input can be approximated to any model probability density function.To further investigate the Monte-Carlo method in the evaluation of measurement uncertainty of the application,the paper analyzes the implementation steps of uncertainty evaluation,and derive the formula of random variables corresponding to its probability density.Example show that trials between one hundred thousand to one million times,the result of using the Monte-Carlo method to calculate the estimated value and standard uncertainty are stable.When the trials number are larger,the interval containing of probability symmetrical of 95%,and the accuracy is higher.

Monte-Carlo method；hydraulic fracturing；ceramic proppant；flow conductivity；measurement uncertainty

TE312

A

1673-5285（2017）09-0067-06

10.3969/j.issn.1673-5285.2017.09.017

2017－08-04

陶祖文，男（1989-），助理工程師，畢業(yè)于西南石油大學(xué)，主要從事石油工程現(xiàn)場(chǎng)應(yīng)用及科研工作。