李文軍,孫宏健,鄭永軍

(中國計量大學(xué)計量測試工程學(xué)院,杭州 310018)

熱電偶被廣泛應(yīng)用于工業(yè)領(lǐng)域。當(dāng)熱電偶用于測量動態(tài)溫度時,需要對熱電偶動態(tài)響應(yīng)特性進(jìn)行評價[1-2]。熱電偶動態(tài)響應(yīng)特性評價方法通常遵循溫度傳感器動態(tài)特性評價方法,簡化為熱電偶時間常數(shù)的測試[3]。

熱電偶時間常數(shù)的測試方法有多種。根據(jù)溫度激勵產(chǎn)生方式的不同,常用方法有熱風(fēng)洞法、電加熱法、激波管法、投入法和激光加熱法[4]。在這些測試方法中,都存在一個比較困難的環(huán)節(jié),即如何確定溫度激勵起始時間點(diǎn)。對于投入法,通過激光對射開關(guān)等輔助測量手段,測量投入介質(zhì)時的時間點(diǎn)并作為起始時間[5]。對于激光加熱法,是用響應(yīng)更快的紅外探測設(shè)備輔助測量確定時間點(diǎn)[6]。但是這兩種輔助測量手段本身又引入了誤差,當(dāng)評價小時間常數(shù)快響應(yīng)熱電偶時,這種誤差影響變大[7-8]。

因此需要尋找其他的熱電偶動態(tài)特性評價方法,比如采用神經(jīng)網(wǎng)絡(luò)技術(shù)對校準(zhǔn)數(shù)據(jù)進(jìn)行建模[9]。類似地,如果熱電偶測溫視為一個動態(tài)過程,則可以采用辨識和回歸技術(shù)建模以考察熱電偶的動態(tài)響應(yīng)。建模可以從理論建模入手,也可以從實(shí)驗(yàn)建模入手,也可以從兩者的結(jié)合入手。本文從實(shí)驗(yàn)建模的角度,建立了熱電偶動態(tài)響應(yīng)實(shí)驗(yàn)系統(tǒng),采集了熱電偶動態(tài)響應(yīng)數(shù)據(jù)。把熱電偶響應(yīng)數(shù)據(jù)視為時間序列,采用貝葉斯回歸方法估計參數(shù)分布?？紤]到實(shí)驗(yàn)數(shù)據(jù)樣本容量的不足,采用了馬爾科夫鏈蒙特卡羅法MCMC(Markov Chain Monte Carlo)抽樣。

1 熱電偶動態(tài)響應(yīng)及其實(shí)驗(yàn)評價

1.1 熱電偶動態(tài)響應(yīng)的影響因素

熱電偶動態(tài)響應(yīng)受到多種因素影響。從實(shí)驗(yàn)現(xiàn)象看,熱電偶動態(tài)特性受熱流密度影響較大[10-11]。另外,熱電偶在使用中本身也存在溫度漂移,以裸露式K型熱電偶為例,在室溫到600 ℃區(qū)間內(nèi),塞貝克系數(shù)會產(chǎn)生變化,變化引起的溫度誤差在0.5%[12]。以鎧裝式K型熱電偶為例,在室溫到600 ℃區(qū)間內(nèi),溫度的漂移在2.5 ℃左右[13]。

從物理過程分析,當(dāng)熱電偶測量流體溫度時,影響熱電偶與流體之間非穩(wěn)態(tài)換熱的因素有:對流以及流動狀態(tài)、流體的物理性質(zhì)以及換熱面形狀、幾何尺寸和相對位置[14]。如果流體為氣體,則還存在熱電偶與氣體之間的輻射換熱,以及熱電偶的熱端指向冷端的熱傳導(dǎo)[15],如圖1所示。

圖1 熱電偶測氣流溫度時的換熱模型

當(dāng)熱電偶測量固體表面溫度時,影響熱電偶與固體表面之間非穩(wěn)態(tài)換熱的因素有:接觸壓力、接觸面間隙介質(zhì)、接觸面表面粗糙度和接觸面形變[16]。

因此,無論是測量流體溫度或者測量固體溫度,都存在多種因素影響熱電偶與被測物體之間的非穩(wěn)態(tài)換熱過程,進(jìn)而影響到熱電偶的動態(tài)響應(yīng)。從物理過程建模時,選取和舍棄影響因素以及建立條件假設(shè)都比較困難。于是考慮從測量數(shù)據(jù)入手建模,為此建立了熱電偶測量氣體溫度時的實(shí)驗(yàn)評價系統(tǒng)。

圖2 實(shí)驗(yàn)評價系統(tǒng)組成示意圖

1.2 熱電偶動態(tài)響應(yīng)的實(shí)驗(yàn)評價系統(tǒng)

熱電偶測量氣體溫度的實(shí)驗(yàn)評價系統(tǒng)如圖2所示,主要包括計算機(jī)、數(shù)據(jù)采集卡、可編程控制器、往復(fù)運(yùn)動驅(qū)動機(jī)構(gòu)、風(fēng)道和工業(yè)調(diào)溫風(fēng)機(jī)。實(shí)驗(yàn)中,利用工業(yè)調(diào)溫風(fēng)機(jī)產(chǎn)生熱氣流,利用可編程控制器控制往復(fù)運(yùn)動驅(qū)動機(jī)構(gòu),驅(qū)動熱電偶進(jìn)出風(fēng)道,以獲得溫度激勵。熱電偶的輸出信號由數(shù)據(jù)采集卡返回給計算機(jī)。

1.3 鎳鉻鎳硅熱電偶的響應(yīng)數(shù)據(jù)集

以一種K型鎳鉻鎳硅熱電偶為對象,利用實(shí)驗(yàn)評價系統(tǒng)取得了幾種溫度激勵下的響應(yīng)數(shù)據(jù),表1給出了一組實(shí)驗(yàn)評價時的設(shè)置參數(shù)。

表1 實(shí)驗(yàn)評價的參數(shù)設(shè)置

按照上述設(shè)置中激勵溫度進(jìn)行實(shí)驗(yàn),獲得了熱電偶響應(yīng)數(shù)據(jù)集。圖3給出了鎳鉻鎳硅熱電偶在3種溫度激勵下的響應(yīng)。

圖3 鎳鉻鎳硅熱電偶溫度響應(yīng)散點(diǎn)圖

由于通過實(shí)驗(yàn)獲得各種溫度階躍下的測試數(shù)據(jù)并不容易,因此評價熱電偶響應(yīng)能力要在測試數(shù)據(jù)不充分的條件下進(jìn)行。當(dāng)測試數(shù)據(jù)數(shù)量充分時,一般可采用矩法、極大似然法等經(jīng)典推斷方法,當(dāng)測試數(shù)據(jù)不充分時,可采用小樣本的推斷方法[17]。把圖3中熱電偶溫度響應(yīng)視為時間序列,可采用貝葉斯推斷方法,分析和評價熱電偶響應(yīng)能力[18]。

2 熱電偶動態(tài)響應(yīng)模型的貝葉斯推斷

2.1 貝葉斯推斷

熱電偶的動態(tài)響應(yīng)可以視為一個動態(tài)系統(tǒng)。對于動態(tài)系統(tǒng),可以基于一組測量值或觀察值,用推斷方法來對參數(shù)做估計。推斷的具體描述是,把實(shí)驗(yàn)中獲得的熱電偶動態(tài)響應(yīng)數(shù)據(jù)T視為觀測樣本,對熱電偶動態(tài)響應(yīng)的模型參數(shù)θ進(jìn)行估計。經(jīng)典參數(shù)估計方法如矩法估計和最大似然估計,把參數(shù)θ視為參數(shù)空間Θ的一個未知但是存在的固定值。而從隨機(jī)過程視角看,參數(shù)θ是參數(shù)空間Θ中取值的隨機(jī)變量,可以采用貝葉斯推斷方法來估計參數(shù)θ的分布。

貝葉斯推斷的核心工具是貝葉斯準(zhǔn)則,按照貝葉斯準(zhǔn)則,參數(shù)和溫度響應(yīng)數(shù)據(jù)服從條件概率:

p(θ,T)=p(T|θ)p(θ)

(1)

p(θ,T)=p(θ|T)p(T)

(2)

按照貝葉斯準(zhǔn)則有:

(3)

(4)

式(1)～(4)中:θ為模型參數(shù);T為熱電偶響應(yīng)數(shù)據(jù);p(θ,T)為模型參數(shù)和數(shù)據(jù)的聯(lián)合概率密度函數(shù);p(T|θ)為似然函數(shù),即在給定一個模型時描述這個數(shù)據(jù)的似然;p(θ)為模型參數(shù)的先驗(yàn)分布,可以視為看到數(shù)據(jù)之前對模型參數(shù)分布的一種描述;p(θ|T)為后驗(yàn)分布;p(T)為邊界似然,是對所有似然概率的積分。

根據(jù)貝葉斯準(zhǔn)則,式(3)和式(4)給出了用熱電偶動態(tài)響應(yīng)模型的先驗(yàn)知識對模型進(jìn)行推斷的一種框架。按照這種框架,不直接估計參數(shù)值,而是考慮參數(shù)服從一定概率分布。貝葉斯推斷熱電偶響應(yīng)模型參數(shù)的過程可以概括如下:

首先確定模型參數(shù)的先驗(yàn)分布p(θ);其次獲取一系列觀測T;再次求出參數(shù)θ的后驗(yàn)分布p(θ|T),求出θ的期望值作為其最終值。推斷過程中,還需要定義參數(shù)的一個方差來評估推斷的準(zhǔn)確程度。

2.2 熱電偶動態(tài)響應(yīng)的回歸模型

對于圖3所示的熱電偶動態(tài)響應(yīng)數(shù)據(jù)集,為了計算方便采用過余溫度進(jìn)行轉(zhuǎn)換,令:

(5)

圖4 鎳鉻鎳硅熱電偶無量綱溫度響應(yīng)散點(diǎn)圖

圖4給出了通過轉(zhuǎn)換后鎳鉻鎳硅熱電偶無量綱溫度響應(yīng)的散點(diǎn)圖,以時間序列的形式列出了3種激勵溫度下的響應(yīng)。

圖4中的無量綱溫度時間序列包含了溫度激勵未產(chǎn)生時即熱電偶處于環(huán)境溫度時的狀態(tài)值。當(dāng)評價的數(shù)據(jù)集包含這部分狀態(tài)值時,數(shù)據(jù)處理中可以避免確定溫度激勵產(chǎn)生時間點(diǎn)。根據(jù)圖4中散點(diǎn)的分布形狀,采用如下形式的函數(shù)做回歸:

(6)

式中:t為時間,單位為s;α,β為回歸系數(shù)。

式(6)是一個非線性回歸模型,可以先做線性化變化,令:

(7)

t*=t

(8)

則式(6)轉(zhuǎn)化為:

T*=α+βt*

(9)

式(9)是一個線性回歸函數(shù),其中α為截距,β為回歸系數(shù)。為簡便起見,下文中省略*號,仍用T表示經(jīng)過式(5)無量綱化和經(jīng)過式(6)線性化后的溫度,于是回歸函數(shù)表示為:

T=α+βt

(10)

再設(shè)噪聲方差為σ2,并假設(shè)溫度為獨(dú)立高斯分布[19],即:

Ti|θ～N(μi(θ),σ2)

(11)

式中:μi表示高斯分布均值。

把μi表示為時間預(yù)報值ti和參數(shù)α和β的函數(shù),有:

(12)

假設(shè)α和β的先驗(yàn)分布都是高斯分布,即:

(13)

(14)

為方便,記噪聲方差σ2的對數(shù)為待估計參數(shù),并假設(shè)其為高斯分布,即:

logσ2～N(κ,ω2)

(15)

于是要推斷的參數(shù)可以表示為:

θ=θ(α,β,logσ2)

(16)

式(10)確定了回歸函數(shù)形式和待估計參數(shù),式(13)～式(15)確定了待估計參數(shù)的先驗(yàn)。

由于后驗(yàn)分布通常不可解,所以一般采用各種近似貝葉斯推理方法,比如變分方法和MCMC方法。變分方法的思想是把一個待解問題,引入一個變量轉(zhuǎn)變成一個優(yōu)化問題。變分方法在優(yōu)化問題中引入了約束,讓問題簡化以達(dá)到快速求解,但是如果引入的約束比較嚴(yán)格,近似的程度就會變差。而MCMC方法的思想是基于隨機(jī)模擬,即構(gòu)造一個隨機(jī)過程來逐漸逼近所求的分布,通過隨機(jī)過程不斷采樣達(dá)到刻畫目標(biāo)分布的目的。本文選擇了MCMC方法來估計后驗(yàn)分布。

2.3 雜交蒙特卡羅抽樣方法

在貝葉斯估計中,可以采用MCMC方法來估計后驗(yàn)分布。MCMC方法的思想是通過生成隨機(jī)樣本并用樣本來估計后驗(yàn)分布。MCMC生成樣本的有多種,其中雜交蒙特卡羅HMC(Hamiltonian Monte Carlo)方法是一種快速抽樣方法。雜交蒙特卡羅方法的基本思想是在MCMC抽樣中利用分子動力學(xué)模擬MD(Molecular Dynamics)方法來產(chǎn)生馬爾科夫鏈移動[20]。用這個方法所產(chǎn)生的建議分布,往往與目標(biāo)分布的動力學(xué)機(jī)理更接近[21]。

從方法上看,分子動力學(xué)模擬是對漢密爾頓方程做積分求解來估計復(fù)雜物理系統(tǒng)的典型特征,這種典型特征通常是組態(tài)函數(shù)(一般是時間的函數(shù))關(guān)于時間的平均。蒙特卡羅模擬則是從由系統(tǒng)勢能和溫度所決定的玻爾茲曼分布中抽取典型組態(tài)的樣本。分子動力學(xué)模擬和蒙特卡羅模擬之間的聯(lián)系在于,物理系統(tǒng)中時間的平均會收斂到狀態(tài)的平均。因此,分子動力學(xué)模擬產(chǎn)生的結(jié)果往往和蒙特卡洛模擬產(chǎn)生的估計通常會有較好的吻合。

把雜交蒙特卡羅抽樣方法應(yīng)用于熱電偶響應(yīng)分布的抽樣,其具體過程描述如下。

假設(shè)從分布π(T)∝exp[-U(T)]中抽取樣本,首先引入和T共軛的虛擬動量變量p,其次引入一個引導(dǎo)漢密爾頓函數(shù)H′(T,p):

H′(T,p)=U′(T)+k(p)

(17)

它用來產(chǎn)生一個建議狀態(tài)。

再定義接收漢密爾頓函數(shù)H(T,p):

H(T,p)=U(T)+k(p)

(18)

它用來決定是接受還是拒絕。

在定義了引導(dǎo)漢密爾頓函數(shù)和接收漢密爾頓函數(shù)后,雜交蒙特卡羅抽樣通過迭代過程實(shí)現(xiàn),迭代過程用算法1表示。

算法1:

假設(shè)t時刻在組態(tài)空間處于T即T(t)=T,然后在t+1時刻執(zhí)行以下步驟:

第1步 從高斯分布中產(chǎn)生一個新動量p,即p～?(p)∝exp{-k(p)};

第2步 從狀態(tài)(T,p)出發(fā),運(yùn)行Leap-frog算法L步,在相空間中得到一個新狀態(tài)(T′,p′);

第3步 以概率min{1,exp{-H(T′,p′)+H(T,p)}}接收狀態(tài)(T′,p′),以剩余概率令T(t+1)=T。用上述過程可獲得樣本。具體抽樣過程中,還需要確定抽樣的起始點(diǎn)和步長。

3 算例

3.1 雜交蒙特卡羅法抽樣的起始點(diǎn)

對于圖4中的無量綱溫度時間序列,選擇式(10)為回歸公式,選擇參數(shù)為式(13)、式(14)和式(15)形式的先驗(yàn)分布,并預(yù)置參數(shù)的均值和標(biāo)準(zhǔn)差,預(yù)置值如表2所示。

表2 參數(shù)先驗(yàn)分布均值和標(biāo)準(zhǔn)差的預(yù)置值

為計算方便,對式(3)兩邊取自然對數(shù),有:

logp(θ|T)=const.+logp(T|θ)+logp(θ)

(19)

再引入輔助函數(shù)g(θ):

g(θ)=logp(T|θ)+logp(θ)

(20)

對于所建立的線性回歸方程(10),利用上述輔助函數(shù)g(θ)計算似然函數(shù)與先驗(yàn)乘積,并求乘積的對數(shù)。

在進(jìn)行雜交蒙特卡羅抽樣之前,需要確定抽樣起始點(diǎn)?？紤]到抽樣的效率,通常使用最大后驗(yàn)估計值MAP(Maximum a Posteriori Estimation)為起始點(diǎn)。為此,還需要先計算θmap,即:

(21)

計算θmap結(jié)合了式(20)和式(21),即通過計算g(θ)的梯度得到,同步得到其所對應(yīng)的參數(shù)估計值。表3給出了算例中3個參數(shù)的最大后驗(yàn)估計值所對應(yīng)的參數(shù)值。

表3 參數(shù)最大后驗(yàn)估計值對應(yīng)的參數(shù)值

在確定了起始點(diǎn)后,抽樣過程中還需要調(diào)整步長。調(diào)整步長通過一個自適應(yīng)抽樣過程完成。

3.2 雜交蒙特卡羅法抽樣的步長調(diào)整

在熱電偶響應(yīng)的參數(shù)后驗(yàn)分布抽樣過程中,為了能提高抽樣效率,通過自適應(yīng)方法調(diào)整雜交蒙特卡羅抽樣的步長以及相應(yīng)的虛擬動量向量p。通常把算法1中第3步的目標(biāo)接受率初始值設(shè)為0.65,并把調(diào)整步長的起始點(diǎn)也定在最大后驗(yàn)估計值這點(diǎn)。當(dāng)?shù)^程收斂時可以得到目標(biāo)步長。圖5給出了算例中通過一個自適應(yīng)過程確定步長的收斂過程。其中接受率終值為0.64。

圖5 自適應(yīng)抽樣中步長的收斂過程

3.3 熱電偶動態(tài)響應(yīng)模型的回歸結(jié)果

在確定了選擇起始點(diǎn)和步長的方法后,使用獨(dú)立的馬爾科夫鏈從后驗(yàn)分布中抽取樣本。抽取樣本時選擇互異的馬爾科夫鏈初始點(diǎn),且初始點(diǎn)隨機(jī)分布在最大后驗(yàn)估計值θmap周圍。利用老化(burn-in)過程,舍棄開始階段的老化樣本,獲得所使用樣本。在得到熱電偶響應(yīng)模型的參數(shù)后驗(yàn)分布后,再利用標(biāo)準(zhǔn)MCMC診斷給出評價[22]。表4給出了算例中熱電偶響應(yīng)參數(shù)回歸結(jié)果和結(jié)果檢驗(yàn)數(shù)據(jù)。

表4 回歸結(jié)果和結(jié)果檢驗(yàn)

在表4中,Mean為參數(shù)的后驗(yàn)均值;MCSE為Monte Carlo standard error,即蒙特卡羅標(biāo)準(zhǔn)誤,是后驗(yàn)均值的標(biāo)準(zhǔn)差;SD為Standard error,即后驗(yàn)標(biāo)準(zhǔn)差;Q5為邊際后驗(yàn)分布的5分位;Q95為邊際后驗(yàn)分布的95分位;RHat為Gelman-Rubin 收斂診斷值,用來刻畫平行的馬爾科夫鏈的鏈內(nèi)與鏈間方差的差異,當(dāng)收斂診斷值小于等于1.1則說明抽樣算法收斂[23]。在表5中看到,算例中的RHat值都小于1.1,可認(rèn)為收斂到了期望的分布。

圖6 參數(shù)的軌跡圖

3.4 抽樣樣本的分布

在回歸收斂的情形下,考察算例中的抽樣樣本是否構(gòu)成一個目標(biāo)分布上的合理集合。圖6分別給出了算例中3個參數(shù)的第1個馬爾科夫鏈軌跡,其中已經(jīng)舍棄了老化樣本。從圖6可以看出,樣本與樣本之間沒有長期相關(guān)性,樣本混合比較合理。

圖7給出了算例中3個參數(shù)的馬爾科夫鏈混合后的散點(diǎn)圖和直方圖。其中直方圖給出了參數(shù)的邊際后驗(yàn)分布。

圖7 參數(shù)的邊際后驗(yàn)分布

如果采用參數(shù)的后驗(yàn)均值來表示,算例中的回歸公式可以表示如下:

(22)

式(22)是在氣流環(huán)境下,環(huán)境溫度為29 ℃,激勵溫度區(qū)間為95 ℃～200 ℃時,一種K型鎳鉻鎳硅熱電偶無量綱溫度動態(tài)響應(yīng)的回歸方程。

4 結(jié)論

針對熱電偶動態(tài)響應(yīng)能力的評價,建立了用貝葉斯推斷方法計算參數(shù)分布的一種框架。對一種鎳鉻鎳硅熱電偶,用貝葉斯推斷方法計算了其動態(tài)響應(yīng)模型的參數(shù)分布。上述工作為熱電偶動態(tài)響應(yīng)能力的評價提供了一種測試方法。同時,利用上述框架,可以根據(jù)熱電偶先驗(yàn)信息和所建回歸模型進(jìn)行預(yù)測,然后利用測量值和模型的后驗(yàn)信息更新預(yù)測值,再把后驗(yàn)信息置為新的先驗(yàn)信息再次進(jìn)行預(yù)測。因此,這種推斷過程還可應(yīng)用于熱電偶測量動態(tài)溫度的遞推估計中。

[1] 吳朋,林濤. 基于QGA_SVM的鎧裝熱電偶傳感器辨識建模研究[J]. 儀器儀表學(xué)報,2014,35(2):343-349.

[2] 吳德會,趙偉,黃松嶺,等. 傳感器動態(tài)建模FLANN方法改進(jìn)研究[J]. 儀器儀表學(xué)報,2009,30(2):362-367.

[3] Doghmane M Y,Lanzetta F,Gavignet E. Dynamic Characterization of a Transient Surface Temperature Sensor[J]. Procedia Engineering,2015,120:1245-1248.

[4] 馬勤弟,雷敏. 薄膜熱電偶的動態(tài)校準(zhǔn)及辨識建模[J]. 儀器儀表學(xué)報,1999,20(3):300-302.

[5] 趙學(xué)敏,王文廉,李巖峰,等. 火焰溫度場測試中的傳感器動態(tài)響應(yīng)研究[J]. 傳感技術(shù)學(xué)報,2016,29(3):368-372.

[6] 郝曉劍,李科杰,劉健,等. 基于CO2激光器的溫度傳感器可溯源動態(tài)校準(zhǔn)[J]. 兵工學(xué)報,2009,30(2):156-159.

[7] Diniz A C G C A,Almeida M E K D,Vianna J N S,et al. Methodology for Estimating Measurement Uncertainty in the Dynamic Calibration of Industrial Temperature Sensors[J]. Journal of the Brazilian Society of Mechanical Sciences and Engineering,2017,39:1053-1060.

[8] 趙學(xué)敏,王文廉,李巖峰,等. 火焰溫度場測試中的傳感器動態(tài)特性補(bǔ)償[J]. 傳感技術(shù)學(xué)報,2017,30(5):735-741.

[9] Danisman K,Dalkiran I,Celebi F V. Design of A High Precision Temperature Measurement System Based on Artificial Neural Network for Different Thermocouple Types[J]. Measurement,2006,39:695-700.

[10] 成紅剛,陳雄,鞠玉濤,等. 基于熱電偶動態(tài)特性的溫度預(yù)估方法實(shí)驗(yàn)研究[J]. 固體火箭技術(shù),2012,35(6):842-846.

[11] 曾磊,桂業(yè)偉,賀立新,等. 鍍層式同軸熱電偶數(shù)據(jù)處理方法研究[J]. 工程熱物理學(xué)報,2009,30(4):661-664.

[12] Webster E S. Drift in Type K Bare-Wire Thermocouples from Different Manufacturers[J]. International Journal of Thermophysics,2017,38(5):1-14.

[13] Webster E S. Thermal Preconditioning of MIMS Type K Thermocouples to Reduce Drift[J]. International Journal of Thermophysics,2017,38(1):1-14.

[14] Jaremkiewicz M. Accurate Measurement of Unsteady State Fluid Temperature[J]. Heat Mass Transfer,2017,53:887-897.

[15] 楊永軍,蔡靜. 特殊條件下的溫度測量[M]. 北京:中國計量出版社,2009:130-132.

[16] Verma N N,Mazumder S. Extraction of Thermal Contact Conductance of Metal-Metal Contacts from Scale-Resolved Direct Numerical Simulation[J]. International Journal of Heat and Mass Transfer,2016,94:164-173.

[17] 姚繼濤,王旭東. 小樣本條件下可變作用代表值的貝葉斯推斷方法[J]. 西南交通大學(xué)學(xué)報,2014,49(6):995-1001.

[18] Kantz H,Schreiber T. Nonlinear Time Series Analysis[M]. UK:Cambridge University Press,2004:22-145.

[19] Chakraborty S,Das P K. Application of Bayesian Inference Technique for the Reconstruction of An Isothermal Hot Spot inside A Circular Disc from Peripheral Temperature Measurement—A Critical Assess-ment[J]. International Journal of Heat and Mass Transfer,2015,88:456-469.

[20] 劉軍. 科學(xué)計算中的蒙特卡羅策略[M]. 北京:高等教育出版社,2009:140-141.

[21] 康崇祿. 蒙特卡羅方法理論和應(yīng)用[M]. 北京:科學(xué)出版社,2015:115-116.

[22] 劉金山. 基于MCMC算法的貝葉斯統(tǒng)計方法[M]. 北京:科學(xué)出版社,2016:49-53.

[23] Brooks S P,Gelman A. General Methods for Monitoring Convergence of Iterative Simulations[J]. Journal of Computational and Graphical Statistics,1998,7(4):434-455.