基于塔克分解的濃度重建方法

王才謙 劉石

1.華北電力大學(xué)能源動(dòng)力與機(jī)械工程學(xué)院；2.華北電力大學(xué)控制與計(jì)算機(jī)工程學(xué)院

針對(duì)流場(chǎng)中存在復(fù)雜機(jī)械結(jié)構(gòu)而導(dǎo)致得濃度難以測(cè)量的問題，提出了基于塔克分解的濃度重建算法。采用塔克分解對(duì)樣本數(shù)據(jù)庫進(jìn)行特征提取，將其運(yùn)用于濃度分布重建。為了驗(yàn)證算法的可行性，進(jìn)行了泄漏仿真實(shí)驗(yàn)。仿真實(shí)驗(yàn)結(jié)果表明濃度重建算法可以快速重建濃度分布。最后將本方法與基于主成分分析的重建方法進(jìn)行比較，實(shí)驗(yàn)數(shù)據(jù)表明本方法在三維濃度重建上具有更好的重建精度。

流場(chǎng)中的濃度信息可以通過實(shí)驗(yàn)測(cè)量和數(shù)值計(jì)算獲得。其中實(shí)驗(yàn)測(cè)量主要可分為接觸式與非接觸式兩大類。接觸式測(cè)量通過將傳感器放置在流場(chǎng)中進(jìn)行采樣測(cè)量，只能獲取流場(chǎng)少數(shù)位置的濃度值，無法獲取全場(chǎng)的瞬時(shí)濃度。若要進(jìn)行多點(diǎn)同時(shí)采樣，需要大量的采樣儀器與分析設(shè)備，測(cè)量成本高。非接觸式測(cè)量利用流體的光學(xué)[1]或者聲學(xué)[2]特性進(jìn)行濃度測(cè)量，但是這些設(shè)備也會(huì)有結(jié)構(gòu)復(fù)雜、布置困難、易受環(huán)境影響等問題，一定程度上制約了非接觸時(shí)測(cè)量的應(yīng)用。近年來，隨著計(jì)算流體力學(xué)的發(fā)展與計(jì)算機(jī)水平的進(jìn)步，數(shù)值計(jì)算廣泛地運(yùn)用于流場(chǎng)分析中。但是要想求解流場(chǎng)的濃度值信息，需要獲取流場(chǎng)的邊界條件，因此在邊界條件未知時(shí)，就難以使用數(shù)值計(jì)算方法。而且使用數(shù)值計(jì)算方法需要大量的計(jì)算資源和時(shí)間，很難實(shí)時(shí)獲取流場(chǎng)的濃度值信息。

本文的主要工作是將數(shù)值計(jì)算方法和測(cè)量方法結(jié)合起來，提出一種基于塔克分解的濃度重建算法。該算法將CFD 數(shù)值模擬得到的流場(chǎng)濃度值作為先驗(yàn)信息，然后將先驗(yàn)信息降維提取主要特征后，結(jié)合測(cè)量的少量濃度值，重建出全場(chǎng)的濃度分布。

1 算法原理

基于特征提取的數(shù)據(jù)處理方法最早運(yùn)用于屋蓋風(fēng)壓場(chǎng)的重建[3]，后來又廣泛運(yùn)用于其他流場(chǎng)參數(shù)的重建過程中：QIN 等人[4]將主成分分析、CFD 計(jì)算數(shù)據(jù)結(jié)合起來，實(shí)現(xiàn)了對(duì)風(fēng)場(chǎng)分布的快速預(yù)測(cè)；Sun 等人[5]提出一種基于本征正交分解的風(fēng)速重建方法，并且對(duì)測(cè)量點(diǎn)位置進(jìn)行了優(yōu)化，減少了重建誤差；陳敏鑫等人[6]將主成分分析引入到溫度分布重建工作中，利用經(jīng)驗(yàn)數(shù)據(jù)和少量測(cè)量數(shù)據(jù)快速準(zhǔn)確地完成了溫度分布重建。但是以上方法研究的均是二維平面上流場(chǎng)參數(shù)的重建，所使用的均是基于向量或矩陣的數(shù)據(jù)降維方法。如果將這些方法運(yùn)用于三維流場(chǎng)的重建，會(huì)在一定程度上造成三維空間特征的損失，為了避免這一情況出現(xiàn)，本文引入了塔克分解作為降維算法對(duì)先驗(yàn)信息進(jìn)行特征提取。

1.1 塔克分解

隨著對(duì)數(shù)據(jù)的深入研究，矩陣已經(jīng)很難完整表述出數(shù)據(jù)的時(shí)空關(guān)系。而張量作為可以表征高階數(shù)據(jù)的多維數(shù)組，可以完整地保留各維度的信息。張量的構(gòu)成如圖1所示。

圖1 張量示意圖Fig.1 Tensor diagram

1963 年Tucker 首次提出了塔克分解，將一個(gè)n 階高維的張量，得到n 個(gè)分解因子矩陣和一個(gè)n 階低維的核心張量。n 階高維張量的分解過程的公式如式（1）所示：

圖2 三階張量塔克分解示意圖Fig.2 Schematic diagram of the third-order tensor Tucker decomposition

1.2 濃度分布重建算法

利用經(jīng)驗(yàn)數(shù)據(jù)或者測(cè)量數(shù)據(jù)獲得n 個(gè)工況下的三維濃度分布，組成一份先驗(yàn)數(shù)據(jù)集XT={X1，X2，…，。對(duì)先驗(yàn)數(shù)據(jù)集XT 進(jìn)行塔克分解將其分解成核心張量與分解因子矩陣模態(tài)積的形式，如式（2）所示。

將XT 與F(4)內(nèi)的元素展開，公式（2）可轉(zhuǎn)化為如式（3）所示形式：

式中，bi為4 模態(tài)方向上的分解因子矩陣中的第i 個(gè)因子向量，

由公式（3）可得：

那么對(duì)于待重建的三維濃度分布XR,如式（5）所示：

式中，b 為4 模態(tài)方向上待重建的三維濃度分布XR所對(duì)應(yīng)的分解因子向量。

根據(jù)模態(tài)積的定義將式（5）從張量形式轉(zhuǎn)化為矩陣與向量相乘的形式，如式（6）所示：

式中，xR為將三維濃度分布XR中的元素轉(zhuǎn)換為向量形式，維度為(n1×n2×n3)×1；Aunfold為將重建系數(shù)張量沿著模態(tài)4 方向展開的矩陣形式，維度為(n1×n2×n3)×m4。

少量的測(cè)量數(shù)據(jù)與整個(gè)流場(chǎng)的濃度分布的位置關(guān)系可由式（7）表示。

式中，y 為實(shí)際測(cè)量得到的數(shù)據(jù)，維數(shù)為r×1，其中r為測(cè)點(diǎn)數(shù)目；M 為濃度測(cè)量矩陣，維數(shù)為r×(n1×n2×n3)。濃度測(cè)量矩陣M 為0-1 矩陣，每一行只有一個(gè)元素的值為1，其余元素的值為0。在第n 行（n 取1 到r），值為1 的元素所在的位置，代表了第n 個(gè)測(cè)量點(diǎn)與全部濃度點(diǎn)的相對(duì)位置。

為了使重建的三維濃度分布與實(shí)際的三維濃度分布相等，即x=xR，把式（6）代入到式（7）中，可以得到如式（8）所示：

在式（8）中，y 通過測(cè)量獲得，M 為測(cè)量矩陣，在決定好測(cè)點(diǎn)位置后為已知量，通過對(duì)樣本數(shù)據(jù)張量XT 進(jìn)行塔克分解后獲得，因此可以先用最小二乘法求解出式（8）中的b，再代入式（6）即可求解出待重建的三維濃度分布。

最后為了檢驗(yàn)重建效果，定義相對(duì)重建誤差：

基于塔克分解的濃度重建算法如圖3 所示。從圖3中可以看出，完成重建算法的第一步是利用CFD 進(jìn)行數(shù)值模擬，得到邊界條件下濃度場(chǎng)分布數(shù)據(jù)集；第二步是對(duì)樣本張量進(jìn)行塔克分解，求解其核心張量與各個(gè)模態(tài)方向上的分解因子矩陣；第三步是求解重建系數(shù)張量，并將其展開為重建系數(shù)矩陣；第四步是得到各測(cè)點(diǎn)的位置與測(cè)量的濃度值；最后計(jì)算待重建的濃度分布向量所對(duì)應(yīng)的分解因子向量b，進(jìn)而求出重建結(jié)果。

圖3 基于塔克分解的濃度重建算法流程圖Fig.3 Flow chart of concentration reconstruction algorithm based on tak decomposition

2 案例驗(yàn)證

為了驗(yàn)證算法的可靠性，建立泄漏模型。使用前處理軟件ICEM 建立幾何結(jié)構(gòu)并劃分網(wǎng)格，使用Fluent17.0進(jìn)行仿真求解。

2.1 模型幾何信息

泄漏模型如圖4所示，整個(gè)容器的計(jì)算尺寸為200mm×200mm×875mm，水流入口取水在容器中經(jīng)過穩(wěn)流裝置后的截面，金屬棒下端固定在距水流入口50mm 高的位置上，每個(gè)金屬棒長565mm，為了簡化計(jì)算，固定裝置忽略不計(jì)，水流出口的截面尺寸為60mm×40mm。冷卻水從下方的水流入口進(jìn)入，經(jīng)過金屬棒然后從水流出口流出。中間的金屬棒在高度為420mm 的位置上存在一個(gè)直徑為3mm 的圓形破口，方向朝著x 軸正方向。泄漏物會(huì)由于濃度差以及流體攜帶等機(jī)制向水中擴(kuò)散和傳輸，最終隨水流從出口排出。

圖4 流動(dòng)模型結(jié)構(gòu)圖Fig.4 Flow model structure diagram

從圖4 中可以看出，9 根金屬棒遮擋了整個(gè)容器下半部分的濃度分布，很難直接測(cè)量這些區(qū)域的濃度值，而金屬棒上方的濃度分布是可以利用前文中提到的測(cè)量方式測(cè)量的。因此泄漏模型驗(yàn)證的目的是利用金屬棒上方區(qū)域可測(cè)量的濃度分布來重建出金屬棒遮擋區(qū)域不可測(cè)量的濃度分布。

2.2 網(wǎng)格劃分

結(jié)構(gòu)化網(wǎng)格具有生成速度快、網(wǎng)格質(zhì)量高、數(shù)據(jù)結(jié)構(gòu)簡單等優(yōu)點(diǎn)，但是對(duì)于不規(guī)則區(qū)域的適應(yīng)性較差。由于前文中對(duì)本模型進(jìn)行了簡化，幾何結(jié)構(gòu)并不復(fù)雜，因此使用結(jié)構(gòu)化網(wǎng)格進(jìn)行網(wǎng)格劃分。

一般來說網(wǎng)格數(shù)量越多，計(jì)算結(jié)果就能更加精確，但是計(jì)算耗時(shí)也會(huì)變長。為了在從稀疏到密集的5 組網(wǎng)格中選取一個(gè)合適的網(wǎng)格節(jié)點(diǎn)數(shù)，在保證精度的前提下減少計(jì)算耗時(shí)，對(duì)這5 組網(wǎng)格進(jìn)行網(wǎng)格無關(guān)性檢驗(yàn)。5組網(wǎng)格的網(wǎng)格數(shù)分別為1131852、1729692、2658972、4526216、6142286。選取中間金屬棒頂端中心到水流出口中心所連成的線段上100 個(gè)點(diǎn)的流速作為檢驗(yàn)計(jì)算結(jié)果的判斷標(biāo)準(zhǔn)，不同網(wǎng)格數(shù)目下樣本點(diǎn)流速如圖5 所示。

圖5 不同網(wǎng)格數(shù)下各樣本點(diǎn)流速Fig.5 Flow velocity of each sample point under different grid numbers

觀察樣本點(diǎn)折線可得，使用以上網(wǎng)格所計(jì)算出的結(jié)果差別不大，考慮到在網(wǎng)格數(shù)為1131852 時(shí)計(jì)算速度最快，因此在本文中仿真使用網(wǎng)格數(shù)為1131852 的網(wǎng)格。

2.3 計(jì)算模型

泄漏實(shí)驗(yàn)的泄漏物為羅丹明B 水溶液，由于溶液的濃度不高，且不會(huì)與水發(fā)生化學(xué)反應(yīng)，其物理性質(zhì)仍按照水的性質(zhì)設(shè)置。泄漏物與水都視為不可壓縮流體，對(duì)于水流的湍流流動(dòng)，采用k-ε 湍流模型方程并結(jié)合連續(xù)性方程和動(dòng)量守恒方程進(jìn)行三維仿真。選擇速度、壓力的耦合求解算法，動(dòng)量和湍動(dòng)能采用二階迎風(fēng)格式求解。水流入口和泄漏口設(shè)定為速度入口，水流出口設(shè)定為壓力出口，出口壓力設(shè)為0pa。殘差的收斂條件設(shè)為10e-4，根據(jù)殘差和出口、入口質(zhì)量守恒來判斷計(jì)算是否收斂。

3 重建結(jié)果

本文將重建區(qū)域的三維空間用幾個(gè)堆疊的平面來表示。由于計(jì)算結(jié)果中越靠近z=0m 處的濃度值較高，越靠近z=±0.1m 處的濃度值較低，而濃度值較低的點(diǎn)對(duì)濃度重建的意義不大，因此本文中選取z=-0.04m，z=-0.02m，z=0m，z=0.02m，z=0.04m 五個(gè)平面來代表整個(gè)三維流體域，如圖6 陰影部分所示。

圖6 堆疊平面示意圖Fig.6 Stacked floor plan

每個(gè)平面中濃度點(diǎn)的排列方式為：在x 方向上從x=-0.1m 開始到x=0.1m 結(jié)束，每隔0.002m 設(shè)置一個(gè)濃度點(diǎn)，共設(shè)置101 個(gè)；在y 方向上從y=0.42m 開始到y(tǒng)=0.82m 結(jié)束，每隔0.002m 設(shè)置一個(gè)濃度點(diǎn)，共設(shè)置201個(gè)，每一個(gè)平面的濃度點(diǎn)數(shù)目為20301 個(gè)。本文根據(jù)不同的水流速與泄漏物流速設(shè)置36 個(gè)邊界條件作為樣本工況，水流速分別為0.01、0.02、0.03、0.04、0.05 和0.06m/s，泄漏物流速分別為0.01、0.02、0.03、0.04、0.05 和0.06m/s。將樣本工況計(jì)算完成后，建立一個(gè)101×201×5×36 的四維張量樣本數(shù)據(jù)庫。

另外，在實(shí)際測(cè)量的過程中，只知道測(cè)量點(diǎn)的濃度分布，不知道水流速度與泄漏物泄漏速度這些邊界條件，測(cè)量的工況基本上都是在樣本工況之外。因此，想要驗(yàn)證濃度重建算法的實(shí)用性，還要驗(yàn)證待重建的濃度分布工況位于樣本工況外的重建結(jié)果，因此除了36 個(gè)樣本工況外還設(shè)置了5 個(gè)測(cè)試工況，構(gòu)成了一個(gè)維數(shù)為101×201×5×5 的四維張量測(cè)試數(shù)據(jù)集。

從高度在0.62m 以上的濃度點(diǎn)（即沒有被金屬棒遮擋的區(qū)域的濃度點(diǎn)）中隨機(jī)選取400 個(gè)點(diǎn)作為測(cè)量點(diǎn)數(shù)據(jù)帶入到重建算法中，三維濃度重建所使用的測(cè)試工況如表1 所示。

表1 測(cè)試工況參數(shù)Tab.1 Test condition parameters

經(jīng)過濃度重建后，5 個(gè)測(cè)試工況的重建誤差如表2所示。

表2 測(cè)試工況重建誤差Tab.2 Reconstruction error of test conditions

從表2 可得，各測(cè)試工況的重建誤差均在10%以下，重建時(shí)間在1s內(nèi)，因此只要能夠及時(shí)獲取測(cè)量點(diǎn)的濃度值，就可以獲得全場(chǎng)的瞬時(shí)濃度。重建所得測(cè)試工況5的濃度分布如圖7 所示，由于各平面的濃度值的數(shù)量級(jí)差距較大，故將各平面的濃度分布分開展示。

由圖7 可以看出，在測(cè)試工況不在樣本工況范圍內(nèi)的情況下，使用基于塔克分解的濃度重建分布算法在z=±0.02m 和z=0m 這三個(gè)平面有較好的重建結(jié)果。對(duì)于z=±0.04m 兩個(gè)平面來說，雖然在濃度分布圖上的數(shù)值計(jì)算和重建的結(jié)果差距較大，但是這兩個(gè)平面的濃度值很低，為10-9的數(shù)量級(jí)，雖然相對(duì)重建誤差較大，但是絕對(duì)誤差仍然在一個(gè)較小的水平。因此可以認(rèn)為本算法可以實(shí)現(xiàn)對(duì)模型濃度分布的有效重建。

圖7 各平面數(shù)值計(jì)算與重建濃度分布對(duì)比Fig.7 Comparison of numerical calculation and reconstructed concentration distribution in each plane

此外，在三維濃度分布重建算法的研究中，可以將一個(gè)工況下的三維濃度點(diǎn)進(jìn)行編號(hào)，將其儲(chǔ)存在一個(gè)向量中，計(jì)算多個(gè)工況的濃度值就形成一個(gè)樣本矩陣，再采用基于主成分分析的方法來提取特征向量[7]。為了比較該方法與基于塔克分解方法的優(yōu)劣，另外采用了基于主成分分析的濃度重建算法與基于塔克分解的重建算法進(jìn)行對(duì)比。

如圖8 所示，相較于基于塔克分解的重建算法，使用基于主成分分析的重建算法的重建誤差在各測(cè)試工況下都要更高。而且使用主成分分析的過程中對(duì)先驗(yàn)信息的處理會(huì)涉及到協(xié)方差矩陣的計(jì)算，在先驗(yàn)信息維數(shù)較高時(shí)，計(jì)算機(jī)甚至?xí)霈F(xiàn)內(nèi)存不足無法處理的情況。因此使用塔克分解的重建算法相較于使用主成分分析的重建算法在三維濃度重建方面具有更好的效果。

圖8 各工況下基于主成分分析與塔克分解的重建算法的重建誤差對(duì)比Fig.8 Comparison of reconstruction errors of reconstruction algorithms based on principal component analysis and Tucker decomposition under various working conditions

4 結(jié)語

本文針對(duì)由于技術(shù)條件和環(huán)境影響難以獲取流體域內(nèi)部濃度信息的問題，將特征提取引入濃度分布重建的研究中，結(jié)合了實(shí)驗(yàn)測(cè)量與數(shù)值計(jì)算方法，提出了基于塔克分解的三維濃度重建算法，并從原理上論證了算法的可行性，并總結(jié)了進(jìn)行濃度分布重建的計(jì)算流程。然后對(duì)于泄漏過程進(jìn)行數(shù)值模擬，將數(shù)值模擬的結(jié)果與重建算法得到的結(jié)果進(jìn)行對(duì)比分析，重建結(jié)果表明本文的算法能夠較準(zhǔn)確地重建數(shù)值模擬的數(shù)據(jù)，初步證明了本算法的可行性。最后將本算法與基于主成分分析的濃度重建算法的重建誤差進(jìn)行對(duì)比，結(jié)果表明了本算法相對(duì)于基于主成分分析的濃度重建算法在三維濃度場(chǎng)重建上具有一定的優(yōu)越性。

本文中的重建誤差較大，所使用的測(cè)點(diǎn)數(shù)目也較多，后期研究將對(duì)提升重建精度和降低測(cè)量點(diǎn)數(shù)目進(jìn)行進(jìn)一步探索。

引用

[1] 王雪梅,劉石.基于TDLAS技術(shù)的CO_2濃度測(cè)量[J].化工自動(dòng)化及儀表,2016,43(11):1158-1161.

[2] 丁喜波,陳晨,張任,等.基于超聲波相位差的氣體濃度測(cè)量方法[J].高技術(shù)通訊,2014,24(2):189-192.

[3] 李元齊,沈祖炎.本征正交分解法在曲面模型風(fēng)場(chǎng)重構(gòu)中的應(yīng)用[J].同濟(jì)大學(xué)學(xué)報(bào)(自然科學(xué)版),2006(1):22-26.

[4] QIN L,LIU S,LONG T,et al.Wind Field Reconstruction Using Dimension-reduction of CFD Data with Experimental Validation[J].Energy,2018,151(May 15):272-288.

[5] SUN S X,LIU S,LIU J,et al.Wind Field Reconstruction Using Inverse Process with Optimal Sensor Placement[J].IEEE Transactions on Sustainable Energy,2018,10(3):1290-1299.

[6] 陳敏鑫,劉石,孫單勛,等.機(jī)器學(xué)習(xí)算法在溫度分布重建中的應(yīng)用[J].自動(dòng)化儀表,2020,41(1):95-100+105.

[7] 陳敏鑫.基于降維和深度學(xué)習(xí)方法的溫度分布重建[D].北京:華北電力大學(xué)(北京),2021.