基于角聯(lián)子網(wǎng)的風量反演風阻病態(tài)改良算法

李雨成，李俊橋，鄧存寶，劉蓉蒸

(1.太原理工大學 安全與應急管理工程學院，山西 太原 030024; 2.遼寧工程技術大學 安全科學與工程學院，遼寧 阜新 123000)

摩擦風阻是通風系統(tǒng)重要參數(shù)，是網(wǎng)絡解算的基礎，由于井下大部分故障均可歸結(jié)于風阻改變，其變化也能夠反映某些巷道是否發(fā)生故障[1-2]。目前，風阻主要通過人工測定的方式獲取，測定過程不僅費時費力，也存在諸多弊端。一方面，限于測試條件，某些巷道風阻測不準，甚至無法測定[3-4]。例如:井下巷道風流存在湍流脈動現(xiàn)象[5]，讀數(shù)時不可避免地產(chǎn)生觀測誤差;測定人員技術不熟練，也會影響測試準確度;又如，回風井為立井或預留回風小井的礦井，由于風硐無檢查門或安全問題人員無法測定。另一方面，井下采掘是一個動態(tài)變化的過程，風門開啟、礦車運行、提升設備升降、巷道冒落等都會引起風阻的變化，難以做到風阻的實時更新。相比于人工測試，使用監(jiān)測監(jiān)控數(shù)據(jù)(風量、風壓)對風阻進行反演，不僅能夠?qū)⒓夹g人員從繁重的阻力測定任務中解放出來，還能實現(xiàn)風阻的實時更新，更能提高測試的安全性與準確度，因此是礦井風阻參數(shù)智能采集技術的發(fā)展方向[6-7]。

近年來，許多學者對測風求阻算法進行了研究。1991年，劉澤功[8]首次提出利用通風系統(tǒng)調(diào)風和阻力測定計算通風風阻，并給出了測算步驟與計算實例;2004年，周麗紅等[9]提出使用風量數(shù)據(jù)能夠計算出等同于獨立回路個數(shù)的風阻，降低了風阻參數(shù)獲取的難度;2012年，司俊鴻等[10]提出使用Tikhonov正則化法對測風求阻病態(tài)模型進行修正，可以減輕模型的病態(tài)程度;2015年，李雨成等[11]提出已知兩組風量及部分節(jié)點壓能數(shù)據(jù)，能夠反演未知巷道風阻。至此，前人對測風求阻算法的研究已經(jīng)十分成熟，但應用于大規(guī)模通風網(wǎng)絡時，病態(tài)問題仍然沒有得到徹底的解決，不能滿足精度要求，導致無法實際應用。基于此，筆者引入文獻[11]中的風量反演風阻算法，針對該算法的病態(tài)性質(zhì)，提出將大規(guī)模通風網(wǎng)絡劃分為若干角聯(lián)子網(wǎng)分別求解風阻，解決了算法的病態(tài)問題。同時給出了風壓傳感器優(yōu)化布置算法，結(jié)合井下監(jiān)測監(jiān)控系統(tǒng)，能夠?qū)崟r反演風阻，不僅使礦井風阻參數(shù)智能采集成為了可能，更是未來實現(xiàn)風阻異常分支智能診斷的基礎理論。

1 風量反演風阻算法原理

已知巷道風阻求解風量具有唯一解，而已知巷道風量求解風阻則屬于多解問題，為了使風量反演風阻具有唯一解，可以限定部分節(jié)點壓能。基于節(jié)點壓能的風量反演風阻算法是已知兩組巷道風量及部分節(jié)點壓能，通過構(gòu)建平衡方程組求解未知節(jié)點壓能，根據(jù)通風阻力定律就能實現(xiàn)反演風阻。

對于通風系統(tǒng)任一巷道，能夠得到平衡方程

hi-hj+Gij=rijQijQij

(1)

式中，hi(hj)為巷道節(jié)點靜壓能，Pa;Gij為巷道節(jié)點間位壓差，Pa;rij為巷道風阻，N·s2/m8;Qij為巷道風量，m3/s。

通過調(diào)節(jié)某些構(gòu)筑物的方式，獲取兩組風量。對于風阻不發(fā)生改變的巷道，將風量調(diào)節(jié)前后的平衡方程作比，消去不變量rij，得到風量反演風阻計算式(2)，解出未知節(jié)點壓能，就能根據(jù)通風阻力定律實現(xiàn)風阻反演。

(2)

理論上，對于m條分支，n個節(jié)點的礦井通風系統(tǒng)，設已知壓能節(jié)點個數(shù)為y，風阻發(fā)生改變巷道(一般為獲取兩組風量時調(diào)節(jié)的構(gòu)筑物巷道)個數(shù)為z。為使方程定解，方程個數(shù)應不少于未知數(shù)個數(shù)，即

(3)

已知壓能節(jié)點個數(shù)y滿足式(3)時，風量反演風阻算法能夠得到唯一解。但在算法的實際應用中，筆者發(fā)現(xiàn)該算法具有明顯的病態(tài)特征，即使已知節(jié)點壓能個數(shù)滿足條件，也得不到正確的結(jié)果，計算風阻與實際風阻常常相去甚遠。

2 算法病態(tài)性質(zhì)分析

2.1 算法病態(tài)判別方法

對于線性方程組Ax=b，若矩陣A或向量b有小的擾動，引起解向量x很大的誤差，即解向量x對A與b高度敏感，稱方程組Ax=b為病態(tài)方程組[12]。算法實現(xiàn)過程中，數(shù)據(jù)擾動是不可避免的，因此病態(tài)問題的數(shù)值解是極其不穩(wěn)定的。病態(tài)的界限比較模糊，條件數(shù)是衡量方程組病態(tài)程度的指標，條件數(shù)越大，病態(tài)越嚴重，解的誤差也就越大。為了判斷風量反演風阻算法的誤差是否能夠接受，需要設置條件數(shù)閾值，小于該閾值時，認為計算結(jié)果符合精度。精度要求高時，閾值需要設置的小一些，精度要求不高則可以設置的大一些。本文設置條件數(shù)的數(shù)量級小于103時，計算結(jié)果是非病態(tài)的。條件數(shù)計算方法見式(4)。

Cond(A)=‖A‖·‖A-1‖

(4)

式中，Cond(A)為矩陣A的條件數(shù)，無量綱;‖A‖(‖A-1‖)為矩陣A(A-1)的范數(shù)，無量綱。

2.2 算法病態(tài)原因分析

想要改良風量反演風阻算法病態(tài)問題，需要了解算法病態(tài)的原因。筆者經(jīng)過大量的嘗試與分析，總結(jié)出兩點主要原因:

(1)系數(shù)矩陣稀疏性對算法病態(tài)的影響。

對于大規(guī)模礦井通風系統(tǒng)，分支通常多達數(shù)百條。在構(gòu)建反演風阻方程時，任意風阻不發(fā)生改變的、包含未知壓能節(jié)點的分支都能作為方程組的一個等式;未知節(jié)點壓能構(gòu)成了方程組的未知數(shù)，而同一條分支只能關聯(lián)兩個節(jié)點，即最多只能有4個未知節(jié)點壓能。這說明無論系統(tǒng)規(guī)模多大，系數(shù)矩陣每一行至多只能有4個系數(shù)是非零的。這導致隨著系統(tǒng)規(guī)模的增大，系數(shù)矩陣的稀疏性也將隨之增大，而高維矩陣中元素的不確定性往往會導致矩陣的病態(tài)。可以從降低系統(tǒng)規(guī)模的角度，降低矩陣的稀疏性，進而改善算法的病態(tài)程度。

(2)通風系統(tǒng)拓撲結(jié)構(gòu)對算法病態(tài)的影響。

為了不影響井下實際生產(chǎn)，調(diào)節(jié)構(gòu)筑物的位置和數(shù)量會受到限制，導致獲取兩組風量的手段有局限性。對于不存在調(diào)節(jié)構(gòu)筑物的并聯(lián)巷道，其風量是大致成相同比例變化的，其風比平方系數(shù)Kij非常接近。這會使方程組的行向量過于線性相關，系數(shù)矩陣接近奇異矩陣，最終導致測風求阻算法的病態(tài)。由于系統(tǒng)整體拓撲結(jié)構(gòu)已經(jīng)固定，為了方便求解風阻而調(diào)整系統(tǒng)拓撲結(jié)構(gòu)顯然是不現(xiàn)實的。

3 風量反演風阻病態(tài)改良及風壓傳感器優(yōu)化布置算法

目前，求解病態(tài)方程的手段有很多種，如正則化解法、奇異值分解法[13-14]等。但這些方法都是通過改善條件數(shù)的方法降低算法的病態(tài)程度，不能從根本上解決算法的病態(tài)問題，實際效果并不理想。應用于本文的測風求阻算法時，仍然難以滿足求解精度要求。因此，筆者考慮從算法模型的建立上降低算法的病態(tài)程度。

3.1 基于角聯(lián)子網(wǎng)的風量反演風阻病態(tài)改良算法

角聯(lián)結(jié)構(gòu)普遍存在于通風系統(tǒng)中，具有風流不穩(wěn)定的特點[15-17]。圖1為簡單角聯(lián)結(jié)構(gòu)，記作D={e1(v1,v2),e2(v1,v3),e3(v2,v3),e4(v2,v4),e5(v3,v4)}。

圖1 簡單角聯(lián)結(jié)構(gòu)Fig.1 A simple diagonal structure

將大型礦井通風系統(tǒng)劃分為若干小規(guī)模的、算法模型良態(tài)的角聯(lián)子網(wǎng)，進而對各角聯(lián)子網(wǎng)分別求解風阻。這樣不僅降低了系統(tǒng)規(guī)模，系數(shù)矩陣的稀疏性大大降低，而且角聯(lián)子網(wǎng)構(gòu)建的方程組行向量不具有相關性，因此是解決算法病態(tài)問題的有效手段。算法流程如下。

步驟1:計算通風網(wǎng)絡中所有的角聯(lián)結(jié)構(gòu)，優(yōu)先選擇分支個數(shù)少的角聯(lián)子網(wǎng)作為待求子網(wǎng)，選擇的角聯(lián)子網(wǎng)D1,D2,…,Dn應盡可能覆蓋通風網(wǎng)絡除進、回風井外的所有分支(進、回風井由于其拓撲性質(zhì)，不會包含在角聯(lián)結(jié)構(gòu)中)，即

(5)

式中，EDi為角聯(lián)結(jié)構(gòu)Di的分支集合;∪為計算各集合的并集;E為通風網(wǎng)絡分支集合;Ef為進風井分支集合;Et為回風井分支集合。

步驟2:賦予各節(jié)點權重，根據(jù)節(jié)點權重計算風壓傳感器最優(yōu)布置方案，在3.2節(jié)中將詳細介紹傳感器優(yōu)化布置算法。布置風壓傳感器的節(jié)點作為已知壓能節(jié)點參與后續(xù)風量反演風阻計算。

步驟3:使用風量反演風阻算法計算各角聯(lián)子網(wǎng)分支風阻。對于各子網(wǎng)，當參與計算的分支數(shù)等于未知壓能節(jié)點數(shù)的2倍時，具有唯一解;分支數(shù)大于未知壓能節(jié)點數(shù)的2倍時，具有最小二乘解;分支數(shù)小于未知壓能節(jié)點數(shù)的2倍時，不具有唯一解。

步驟4:檢查、篩選計算結(jié)果。① 舍去條件數(shù)≥103的角聯(lián)子網(wǎng)的病態(tài)計算結(jié)果。② 檢查方程系數(shù)矩陣任一列的非零元個數(shù)。若存在某一列非零元個數(shù)<2，一般存在于包含風阻改變分支的角聯(lián)結(jié)構(gòu)，由于該分支不參與運算，導致該列中非零元對應的節(jié)點壓能是無解的，即使條件數(shù)符合要求，計算結(jié)果也是不可靠的。③ 一條分支存在于多個角聯(lián)子網(wǎng)時，選擇條件數(shù)較小的計算結(jié)果。

理論上，只要風量值和壓能值是準確的，計算風阻可以達到很高的精度。但實際上傳感器觀測數(shù)據(jù)可能不準確，導致誤差增大。由于有些分支存在于多個角聯(lián)子網(wǎng)中，可以利用該特性進行檢查，若某一傳感器故障，會出現(xiàn)各角聯(lián)子網(wǎng)計算結(jié)果不一致的現(xiàn)象。因此，該算法具有一定的自我糾錯的能力。另外，為了盡可能減小誤差，對傳感器數(shù)據(jù)進行預處理也是必要的，如節(jié)點風量平衡修正、數(shù)據(jù)概率平均修正等[5]。

3.2 基于貪心策略的風壓傳感器優(yōu)化布置算法

不合理的風壓傳感器布置方案既不經(jīng)濟、也無法滿足3.1節(jié)中風量反演風阻病態(tài)改良算法的要求。因此，風壓傳感器的優(yōu)化布置是值得研究的。為避免風量反演風阻算法病態(tài)，從控制系統(tǒng)規(guī)模的角度考慮，需要控制方程的系數(shù)矩陣的列數(shù)不超過4，即每一個角聯(lián)子網(wǎng)中的未覆蓋傳感器節(jié)點不能超過兩個。為了使風壓傳感器的布置數(shù)量最少，而且滿足各角聯(lián)子網(wǎng)的計算要求，筆者給出了基于貪心策略的風壓傳感器優(yōu)化布置算法。首先設置節(jié)點權重為節(jié)點在各角聯(lián)子網(wǎng)中出現(xiàn)的次數(shù)，節(jié)點權重越大，說明有更多的子網(wǎng)包含這個節(jié)點。權重可以根據(jù)礦井節(jié)點實際情況進行調(diào)整，能夠使風壓傳感器布置更加合理。算法思路為:逐個選擇節(jié)點，每一步選擇節(jié)點時，優(yōu)先選擇權重最大的節(jié)點;權重相同時，優(yōu)先選擇能夠使布點均勻分布于各角聯(lián)子網(wǎng)的節(jié)點。圖2為基于貪心策略的風壓傳感器優(yōu)化布置算法程序框圖。

圖2 傳感器最優(yōu)布置算法Fig.2 Algorithm of sensor optimized layout

步驟1:初始化角聯(lián)子網(wǎng)節(jié)點集合VD1,VD2,…,VDn;角聯(lián)子網(wǎng)中未覆蓋傳感器節(jié)點最多個數(shù)N=2;待選節(jié)點集合V，備選節(jié)點集合V′，布點集合V″，節(jié)點權重Weight。

步驟2:在V中尋找最大權重節(jié)點集合va。

步驟3:判斷va是否為空集。若不是空集，轉(zhuǎn)步驟4;若是空集，轉(zhuǎn)步驟5。

步驟4:計算va中使varVDj-V″-va[t]最小的節(jié)點，記作va[t]，將其加入布點集合V′，加入該節(jié)點后，可能存在布點個數(shù)滿足要求的子網(wǎng)，將這些子網(wǎng)中未覆蓋傳感器的節(jié)點加入到備選節(jié)點集合V′。其中，var為計算方差，varVDj-V″-va[t] 為衡量布點是否分布均勻的指標，方差越小，說明布點越均勻。在節(jié)點權重相同的節(jié)點中優(yōu)先選擇使布點均勻分布于角聯(lián)子網(wǎng)的節(jié)點。

步驟5:判斷是否存在傳感器個數(shù)未達到要求的子網(wǎng)，其節(jié)點集合記作VDk。若存在，轉(zhuǎn)步驟6;若不存在，轉(zhuǎn)步驟7。

步驟6:在V′與VDk的公共節(jié)點中尋找最大權重節(jié)點集合，若計算結(jié)果包括多個節(jié)點，優(yōu)選使varVDj-V″-va[t] 最小的節(jié)點，將其加入布點集合V″，并從備選節(jié)點V′集合中刪除。

步驟7:返回V″，程序結(jié)束。

4 實例分析

如圖3所示的通風系統(tǒng)包括23條分支、16個節(jié)點、6個構(gòu)筑物、1個進風井、1個回風井。通風方法為抽出式，通風機特性曲線方程為H=-0.052 3Q2+8.560 8Q+3 061.6。為便于計算，假設各節(jié)點標高相同，即不考慮位壓能。調(diào)節(jié)分支e3，e19上的構(gòu)筑物獲得兩組風量，分支e3風阻由30.058 640變?yōu)?.058 640，分支e19風阻由18.028 300變?yōu)?.028 300。分支實際風阻及調(diào)節(jié)構(gòu)筑物前后兩組風量見表1，節(jié)點壓能見表2。限于條件，兩組風量及壓能均是根據(jù)實際風阻進行網(wǎng)絡解算得到，雖然不是傳感器數(shù)據(jù)，但同樣能夠?qū)λ惴ㄟM行驗證。

應用Matlab軟件編寫了計算程序calculResis.m，將兩組風量及部分節(jié)點壓能作為已知條件，分別使用整體直接計算風阻和劃分子網(wǎng)計算風阻兩種算法求解。通過對比計算風阻與實際風阻，驗證直接求解風阻的病態(tài)特征及劃分角聯(lián)子網(wǎng)求解風阻的可行性。

4.1 整體直接求解

根據(jù)式(3)，計算得到滿足求解條件的已知壓能節(jié)點個數(shù)為

圖3 通風網(wǎng)絡Fig.3 Ventilation network graph

表1 分支實際風阻及風量Table 1 Real resistance and air quantity of branches

表2 節(jié)點壓能Table 2 Node pressure

y≥0.5(2n-m+z)=0.5×(2×16-23+2)=5.5

4.2 劃分角聯(lián)子網(wǎng)求解

(1)角聯(lián)子網(wǎng)劃分。

計算通風網(wǎng)絡中所有的角聯(lián)結(jié)構(gòu)，從中選擇了8個分支數(shù)最少的角聯(lián)結(jié)構(gòu)作為待求子網(wǎng)，其中子網(wǎng)分支數(shù)最多不超過9條，覆蓋了通風網(wǎng)絡除進、回風井外的所有分支。選擇的各角聯(lián)子網(wǎng)如下:

D1={e15(v9,v11),e16(v9,v12),e18(v11,v13),

e19(v12,v13),e20(v12,v14),e21(v13,v15),e22(v14,v15)}

D2={e5(v3,v4),e6(v3,v8),e7(v4,v5),e9(v5,v8),

e10(v5,v9),e13(v8,v10),e15(v9,v11),e17(v10,v11)}

D3={e5(v3,v4),e6(v3,v8),e7(v4,v5),e9(v5,v8),

e10(v5,v9),e14(v8,v14),e16(v9,v12),e20(v12,v14)}

D4={e2(v2,v3),e3(v2,v6),e5(v3,v4),e6(v3,v8),

e7(v4,v5),e8(v4,v6),e9(v5,v8)}

D5={e2(v2,v3),e4(v2,v7),e5(v3,v4),e6(v3,v8),

e7(v4,v5),e8(v4,v6),e9(v5,v8),e11(v6,v7)}

D6={e9(v5,v8),e10(v5,v9),e13(v8,v10),e15(v9,v11),

e16(v9,v12),e17(v10,v11),e18(v11,v13),e19(v12,v13)}

D7={e9(v5,v8),e10(v5,v9),e14(v8,v14),e16(v9,v12),

e19(v12,v13),e20(v12,v14),e21(v13,v15),e22(v14,v15)}

D8={e3(v2,v6),e4(v4,v7),e7(v4,v5),e8(v6,v4),

e10(v5,v9),e11(v7,v6),e12(v7,v10),e15(v9,v11),

e17(v10,v11)}

(2)風壓傳感器優(yōu)化布置方案。

計算各子網(wǎng)節(jié)點集合:

VD1={v9,v11,v12,v13,v14,v15}

VD2={v3,v4,v5,v8,v9,v10,v11}

VD3={v3,v4,v5,v8,v9,v12,v14}

VD4={v2,v3,v4,v5,v6,v8}

VD5={v2,v3,v4,v5,v6,v7,v8}

VD6={v5,v8,v9,v10,v11,v12,v13}

VD7={v5,v8,v9,v12,v13,v14,v15}

VD8={v2,v4,v5,v6,v7,v9,v10,v11}

不考慮節(jié)點在通風系統(tǒng)中的實際情況，僅根據(jù)節(jié)點在角聯(lián)子網(wǎng)中出現(xiàn)的次數(shù)，設定節(jié)點權重。v2～v15節(jié)點權重分別為w2=3，w3=4，w4=5，w5=7，w6=3，w7=2，w8=6，w9=6，w10=3，w11=4，w12=4，w13=3，w14=3，w15=2。使用基于貪心策略的風壓傳感器布置算法計算傳感器布置方案，選擇節(jié)點v2，v4，v5，v7，v8，v9，v11，v12，v15為布置風壓傳感器節(jié)點，計算過程見表3，序號0代表初始化。傳感器布置方案完成后，傳感器節(jié)點的壓能作為已知條件參與風阻計算。

表3 傳感器布置方案計算過程Table 3 Computing process of sensor layout

(3)反演風阻。

對各角聯(lián)子網(wǎng)進行風阻求解。計算完成后，檢查計算結(jié)果的可靠性。D4，D7條件數(shù)太大，D6雖然條件數(shù)較小，但方程系數(shù)矩陣結(jié)構(gòu)不合理，因此舍去這3組計算結(jié)果。最終整理計算結(jié)果見表4，可以直觀地看到計算風阻與實際風阻非常接近，誤差控制在10-4以內(nèi)，可以忽略不計。計算結(jié)果不包括進、回風井分支，若要實現(xiàn)進、回風井分支反演，可以在進、回風井口節(jié)點v1，v16增加風壓傳感器，就能實現(xiàn)所有分支風阻的反演。

表4 反演風阻結(jié)果Table 4 Calculation of resistanceN·s2/m8

5 結(jié) 論

(1)系數(shù)矩陣的稀疏性及系統(tǒng)拓撲結(jié)構(gòu)是導致風量反演風阻算法病態(tài)的主要原因。實例中，整體直接求解風阻的8 008種傳感器布置方案中，僅有34.8%的計算結(jié)果非病態(tài)。而對于劃分子網(wǎng)求解，8個子網(wǎng)中僅有2個是病態(tài)的。這也說明算法病態(tài)原因分析是合理的。

(2)針對基于節(jié)點壓能的風量反演風阻算法病態(tài)特征，提出了基于角聯(lián)子網(wǎng)的風量反演風阻病態(tài)改良算法，并給出了計算實例。與整體直接求解風阻相比，雖然增加了風壓傳感器布置數(shù)目，但卻有效解決了病態(tài)問題，誤差控制在10-4以內(nèi)，準確地實現(xiàn)了風阻反演。

(3)提出了基于貪心策略的風壓傳感器優(yōu)化布置算法，能夠制定滿足風量反演風阻計算要求的傳感器最優(yōu)布置方案，為風量反演風阻病態(tài)改良算法與監(jiān)測監(jiān)控系統(tǒng)的結(jié)合提供了經(jīng)濟、有效的方案。