遼寧工業(yè)大學(xué)機(jī)械工程與自動(dòng)化學(xué)院，遼寧錦州 121001

一、引言

目前在儲(chǔ)糧測(cè)溫、爐膛測(cè)溫通常采用聲學(xué)法測(cè)溫方法，但是由于在大部分溫度場(chǎng)中布置傳感器存在邊界位置以及數(shù)目有限的問(wèn)題，導(dǎo)致的弊端是測(cè)溫點(diǎn)密度低，邊界溫度值無(wú)法得出。如果熱點(diǎn)不是恰好出現(xiàn)在內(nèi)部一定范圍中，霉變、蟲(chóng)害通常就會(huì)蔓延到較大區(qū)域才能被發(fā)現(xiàn)，而控制處理措施不及時(shí)將導(dǎo)致儲(chǔ)糧損失增大。

王明吉等[1]以最小二乘方法構(gòu)建了三維溫度場(chǎng)聲學(xué)測(cè)量重建算法，對(duì)球?qū)ΨQ(chēng)型模型溫度場(chǎng)進(jìn)行了仿真重建；白燕[2]基于最小二乘法與傅里葉正則化方法進(jìn)行了重建仿真，研究結(jié)果表明，當(dāng)傳感器數(shù)量增多時(shí)，其重建精度較高；王交峰等人[3]用最小二乘方法對(duì)航空發(fā)動(dòng)機(jī)燃燒室環(huán)形出口溫度場(chǎng)進(jìn)行了重建。這些研究所提算法都要求被測(cè)區(qū)域劃分的單元數(shù)小于系統(tǒng)所獲得的投影數(shù)據(jù)數(shù)，因此原始重建出的溫度點(diǎn)非常少，所帶來(lái)的信息缺失無(wú)法通過(guò)插值運(yùn)算彌補(bǔ)。

顏華等[4]在測(cè)量聲波飛行時(shí)間時(shí)，將采樣信號(hào)的互相關(guān)函數(shù)在峰值附近做三次樣條插值，該方法對(duì)提高聲學(xué)法溫度場(chǎng)檢測(cè)精度具有實(shí)際意義。后來(lái)，顏華又提出了一種徑向基函數(shù)逼近和正則化的溫度場(chǎng)重建算法[5,6]，結(jié)果表明此算法熱點(diǎn)定位精度高，重建算法誤差較小；

鄭永駿[7]運(yùn)用對(duì)偶Kriging模型插值分析氣象資料，擬合4類(lèi)半變異函數(shù)模型對(duì)降水分析精度較高；曾懷恩等[8]提出了Kriging方法運(yùn)用在空間數(shù)據(jù)插值中，通過(guò)實(shí)例與反距離加權(quán)法相比較，證實(shí)了Kriging插值的優(yōu)越性；王婧雅等[9]運(yùn)用改進(jìn)的克里金插值法對(duì)大尺度農(nóng)田土壤墑情分布進(jìn)行了插值，結(jié)果表明，模擬值與實(shí)測(cè)值具有較好的吻合度，兩者的相對(duì)誤差和相關(guān)系數(shù)均接近最優(yōu)值，驗(yàn)證了修正克里金插值法具有可行性和可靠性。綜上，Kriging模型雖應(yīng)用廣泛，但還沒(méi)有應(yīng)用到溫度場(chǎng)重建的內(nèi)外推中。

本文基于現(xiàn)有重建算法弊端及不足，提出了一種基于指數(shù)奇異值分解（(Singular Value Decomposition，SVD）的三維溫度場(chǎng)重建算法，并將Kriging模型引入到溫度場(chǎng)重建的內(nèi)外推中，以正方體溫度場(chǎng)空間為例，采用典型單峰、雙峰、四峰模型溫度場(chǎng)進(jìn)行了重建及插值研究。結(jié)果表明，指數(shù)SVD算法和Kriging模型結(jié)合的重建溫度場(chǎng)達(dá)到較高的精度指標(biāo)。

二、溫度場(chǎng)重建

聲學(xué)法測(cè)溫的主要過(guò)程，首先是在被測(cè)區(qū)域周?chē)M可能均勻的設(shè)置聲波收發(fā)器，任意一個(gè)聲波收發(fā)器發(fā)射信號(hào)，所有聲波收發(fā)器接收信號(hào)，這樣就形成穿過(guò)該被測(cè)區(qū)域的多條均勻分布的聲波飛行路徑，測(cè)量這些穿過(guò)被測(cè)區(qū)域的聲波飛行時(shí)間；然后再利用合適的溫度場(chǎng)重建算法和聲速與溫度的關(guān)系，重建出被測(cè)區(qū)域的溫度場(chǎng)分布[10]。

氣體介質(zhì)中的聲速、氣體介質(zhì)的絕對(duì)溫度和氣體介質(zhì)的聲音常數(shù)之間的關(guān)系為[1,6]：

其中，c—體介質(zhì)中的聲速，單位：m/s；

T—?dú)怏w介質(zhì)的絕對(duì)溫度，單位：K；B—聲音常數(shù)，由氣體組成成分決定的。

在氣體組成成分一定的條件下，B是一個(gè)常量，如被測(cè)氣體為煙道混合氣體時(shí)，B取19.08，被測(cè)氣體為空氣時(shí)，B取20.05[11]。

三、指數(shù)SVD和Kriging模型算法

1、指數(shù)SVD算法

假設(shè)被測(cè)空間溫度場(chǎng)為T(mén)(x,y,z)，聲速的倒數(shù)為f(x,y,z)，由公式（1）可得：

任一聲波路徑pk上，聲波的飛行時(shí)間gk可以表示為：

其中，k—穿過(guò)三維溫度場(chǎng)的有效聲波飛行時(shí)間總數(shù)。

將被測(cè)空間劃分成M個(gè)立方體網(wǎng)格，設(shè)第m個(gè)網(wǎng)格中心點(diǎn)坐標(biāo)表示為(xm,ym,zm)，m∈[1,M]。

將f(x,y,z)離散為M個(gè)基函數(shù)的線性組合：

其中，εm—待定系數(shù)；

φm(x,y,z)—徑向基函數(shù)[12]，有：

其中，a—徑向基函數(shù)的形狀參數(shù)，被測(cè)空間的大小和聲波收發(fā)器的位置都影響形狀參數(shù)的選取。

為了求解式公式（4）中的待定系數(shù)εm，合并式（3）、（4）和（5）得到：

定義：A=(akm)k=1,…,K;m=a,…M；

則公式（6）可寫(xiě)成：

對(duì)重建矩陣A作SVD分解：

其中，σ1,σ2,…,σγ—重建矩陣A的γ個(gè)非零奇異值，σ1≥σ2≥ … ≥σγ≥ 0；

γ—重建矩陣A的秩；

U—列向量，正交矩陣AAT的特征向量；V—列向量，正交矩陣ATA的特征向量。

由此，可以推出A的偽逆為：

通過(guò)奇異值分解，式（7）可以寫(xiě)成

為了增加系統(tǒng)的抗噪聲能力，用一個(gè)濾波函數(shù)對(duì)ε進(jìn)行濾波：

其中，g、a—濾波參數(shù)。

由此，通過(guò)合并（3）式（11）式，可以求的被測(cè)區(qū)域的聲波傳播速度，進(jìn)而代入式（2）就可以求出M個(gè)空間網(wǎng)格的中心點(diǎn)溫度值。

2、Kriging模型插值算法

Kriging模型[12]是一種基于統(tǒng)計(jì)理論的插值技術(shù)。通常Kriging模型變量x=[x1,...,xw]與真實(shí)響應(yīng)y間的關(guān)系可表示為：

式中，f(x) —回歸函數(shù)（一般采用多項(xiàng)式形式）；

λ—回歸系數(shù)；

μ(x) —均值為0、方差為σ2的隨機(jī)函數(shù)，μ(x)的協(xié)方差矩陣為：

式中，ns—采樣點(diǎn)數(shù)；

R—沿對(duì)角線對(duì)稱(chēng)的相關(guān)矩陣；

R(x(i),x(j)) —采樣點(diǎn)x(i)與x(j)的相關(guān)函數(shù)，相關(guān)函數(shù)常用平穩(wěn)高斯函數(shù)來(lái)表述：

其中，θ=(θ1,θ2,…,θw)T—相關(guān)函數(shù)參數(shù)；

w—變量的維數(shù)。

位置x處的響應(yīng)值y(x)的預(yù)測(cè)估計(jì)值為：

其中，λ的估計(jì)值；

y的長(zhǎng)度為ns，包含樣本數(shù)據(jù)的響應(yīng)值；

rT(x)的長(zhǎng)度為ns，是位置x和樣本數(shù)據(jù)(x(1),)間的相關(guān)向量：

設(shè)λ的最優(yōu)估計(jì)值為λ*，全局模型的方差估計(jì)值由λ*和y給出：

相關(guān)函數(shù)參數(shù)θ由極大似然估計(jì)給出，即在θd＞0時(shí)使下式最大：

四、實(shí)驗(yàn)結(jié)果

本文利用聲學(xué)法進(jìn)行溫度場(chǎng)的重建，將重建溫度點(diǎn)進(jìn)行內(nèi)外推的插值，使得整個(gè)被測(cè)區(qū)域的三維溫度場(chǎng)得到更加細(xì)化和邊緣部分無(wú)遺漏的重建推算。

被測(cè)的三維空間為10m × 10m × 10m的正方體，如圖1所示，在每個(gè)頂點(diǎn)及每條棱邊三等分的間隔位置都放置一個(gè)聲波發(fā)生器和一個(gè)聲波接收器，這樣在被測(cè)正方體的周?chē)卜胖?2組聲波收發(fā)器。飛行有效路徑按照剔除棱邊的選取方式，可以得到424條有效聲波路徑（如圖2所示），計(jì)算有效聲波路徑上的聲波飛行時(shí)間。將被測(cè)空間劃分為1000個(gè)均勻的網(wǎng)格，計(jì)算重建矩陣。利用指數(shù)SVD法，重建出1000個(gè)網(wǎng)格的中心點(diǎn)溫度值。

式中，n—被測(cè)區(qū)域所劃分的網(wǎng)格（像素）的總數(shù)；

T(j)和—模型溫度場(chǎng)和重建溫度場(chǎng)第j個(gè)網(wǎng)格中心點(diǎn)的溫度；

Tave和—模型溫度場(chǎng)和重建溫度場(chǎng)的平均溫度。

表1為4種典型溫度場(chǎng)重建后的最大相對(duì)誤差、平均相對(duì)誤差和均方根誤差。雙峰模型的重建溫度場(chǎng)與模型溫度場(chǎng)的三維展示圖如圖3所示。

根據(jù)表1與圖3可知，通過(guò)指數(shù)SVD算法對(duì)4種溫度場(chǎng)重建出的1000個(gè)網(wǎng)格的中心點(diǎn)溫度值與模型中1000個(gè)相應(yīng)位置真實(shí)值的三種誤差結(jié)果比較小。但是被測(cè)溫度場(chǎng)10m × 10m × 10m的范圍內(nèi)，所有邊緣部分的剖分網(wǎng)格外部一半沒(méi)有重建出來(lái)，內(nèi)部的溫度值也不夠細(xì)化分布。

本文運(yùn)用Kriging模型進(jìn)行擬合插值，對(duì)以上指數(shù)SVD算法重建后待解決的填充與細(xì)致化問(wèn)題進(jìn)行研究。對(duì)重建出來(lái)的范圍內(nèi)的溫度值進(jìn)行內(nèi)插與外推，剖分網(wǎng)格由原來(lái)的10 × 10 × 10個(gè)（即1000個(gè)）變?yōu)?1 × 41 × 41個(gè)（即68921個(gè)）。

表2為4種典型溫度場(chǎng)Kriging模型插值后的最大相對(duì)誤差、平均相對(duì)誤差和均方根誤差。圖4所示為經(jīng)過(guò)Kriging模型插值后的雙峰模型的重建溫度場(chǎng)與模型溫度場(chǎng)的三維展示圖。

由表2和圖4可知，指數(shù)SVD算法和Kriging擬合插值后對(duì)4種溫度場(chǎng)重建出的41 × 41 × 41個(gè)（即68921個(gè)）網(wǎng)格的中心點(diǎn)溫度值與模型中41 × 41 ×41個(gè)（即68921個(gè)）相應(yīng)位置真實(shí)值的三種誤差結(jié)果比較小。被測(cè)溫度場(chǎng)10m × 10m × 10m的范圍內(nèi)的所有溫度值能夠較精確的重建出來(lái)，并得到了細(xì)致化展示。

五、結(jié)論

本文提出了指數(shù)SVD法與Kriging模型相結(jié)合的三維溫度場(chǎng)重建方法。采用指數(shù)SVD法對(duì)典型單峰、雙峰、四峰模型溫度場(chǎng)進(jìn)行了重建，kriging模型對(duì)重建溫度場(chǎng)進(jìn)行了內(nèi)插、外推。研究表明：采用指數(shù)SVD和Kriging法對(duì)典型單峰、雙峰、四峰模型溫度場(chǎng)重建，重建最大溫度相對(duì)誤差小于2.2%、平均溫度相對(duì)誤差小于0.2%、溫度均方根誤差小于2%，適用模型廣泛，且重建精度高。

表1 指數(shù)SVD重建誤差

表2 Kriging模型插值誤差