張 巍,李雨成,張 歡,李俊橋,張 靜,李博倫
(太原理工大學 安全與應急管理工程學院,山西 太原 030032)
監測數據[1]是礦井通風系統智能化實現的基石,結構清晰、純度較高的通風基礎數據是風網實時解算、系統優化調節等技術的關鍵。礦井通風智能化系統建設中信息技術的快速發展,使得井下監測監控系統產生數據的規模急劇增長,數據種類不斷豐富。井下傳感器傳回的數據為結構化數據,但礦井通風網絡是一個動態平衡系統,受采掘、運輸、人員活動、地質條件等因素變化的影響。在生產實際中,風速、風質、風壓等類型傳感器監測到的數據污染程度較強[2],這就需要對數據進行挖掘和清洗。數據處理的合理與否,既關系到井下作業環境的優劣,又對提前預防災害的發生起到關鍵作用[3]。
國內外在監測數據處理上取得了一定的研究成果,但尚有很大提高空間。王其軍等[4]通過對多組傳感器數據進行組合,進而提出了評估全部傳感器數據的方法;付華等[5]通過使用神經網絡模型,建立了1種預測模型,給出了基于序列恢復信號理論,處理監測數據異常的方法;馬明煥等[6]通過滑動數據窗模型、加權等價類變換算法等數據挖掘技術構建了一種隱患預警方法;付華等[7]通過集合神經網絡與主元分析方法,正確定位并分離出失效傳感器;王軍號等[8]將傳感器檢測到的異常數據分為偏置型、沖擊型、漂移型和周期型等4種;趙金憲等[9]基于小波分析理論,通過拆解各階段的能量譜,識別監測數據的異常;黃序楨[10]采用均值濾波的方法對傳感器數據進行處理,再使用最小二乘法對數據進行修正處理。能量矢量特征提取主要用于解決大數據條件下監測數據需要分類壓縮等問題[11]。
井下風速傳感器監測到的風速數據具有復雜性、關聯性和層次性。首先,煤礦井下風流在短期內沿著一定的數值上下浮動,數據處理過程既要考慮因運輸、風流短路等原因產生的噪聲,又要對正常波動數值進行標記。其次,多個傳感器回傳的數據具有一定的關聯屬性,降噪算法既要保留數據的原有特征,又要保持2(多)個傳感器數據間的邏輯性。再次,如何開發具有層次降噪能力的通風數據降噪算法也很重要。一些傳感器所處的位置較為核心,其降噪算法以準確性為主,而一些傳感器所處的位置并不核心,其降噪算法以迭代速度等其他特性為主。
本文將模糊C均值聚類算法、魯棒局部加權回歸法和滑動平均法應用到煤礦井下風速傳感器監測數據的處理中。依據算法原理,編寫計算程序,研究不同原理降噪方法對于監測數據處理的特點,分析不同降噪方法對同一對象數據處理異常的原因,得到每一種算法在風速傳感器數據處理上的適用條件和適用場景。研究結果可為礦井通風的異常診斷、災變識別等研究提供合理的基礎數據參數。
模糊C均值(Fuzzy C-means,FCM)算法用于數據聚類分析,其原理是基于特定的目標函數,將風速傳感器Ti周期內監測到的數據集合Xi劃分為c個類,則每個樣本xj屬于某一類i的隸屬度為uij[12-13]。FCM聚類降噪算法的目標函數及約束條件分別為式(1)和式(2):
(1)
(2)
式中:J為目標函數;ci為第i類樣本數據中心;m為隸屬度因子,表示樣本的輕緩程度,一般取2;xj-ci為樣本xj到中心點ci的歐式距離。
目標函數J越小越好,結合約束條件,首先采用Lagrange乘數法建立式(1)~(2)的Lagrange函數,之后對函數中uij,ci,λj等變量依次求偏導數,并使偏導數為0,最終得到變量uij和ci的迭代公式,如式(3)~(4)所示:
(3)
(4)
計算開始時,在開區間(0,1)內隨機生成一uij值,通過uij值計算出ci值,ci值進一步計算uij值,此過程反復迭代,直到目標函數J小于預設精度ε,迭代過程停止,得到最終結果。
1.2.1 算法原理
魯棒局部加權回歸(Robust locally weighted regression,簡稱Rloess)是一種用于局部回歸分析的非參數方法,算法直接從數據特征出發,在回歸擬合之前不指定各變量之間所滿足的函數關系,因此,Rloess做局部降噪處理時具有更明顯的適用性和靈活性。
Rloess降噪算法的原理是把樣本劃分成一個個小區間,對區間中的樣本進行加權多項式擬合,在擬合過程中加入魯棒性的過程,利用絕對中位差MAD賦予數據魯棒權重[14-15],從而剔除離群值,不斷重復這個過程得到回歸與魯棒雙重平滑的曲線,最后再把這些回歸曲線的中心連在一起合成完整的回歸曲線。
1.2.2 計算流程
步驟1:計算區間中每個數據點的回歸權重,權重由式(5)給出:
(5)
式中:x為需要平滑的值;xi為x兩側的第i個值;d(x)為區間長度的2范數。
擬合鄰近點的誤差對擬合效果影響較大,而擬合點較遠處的數值對結果影響最小。此權重函數可以根據實際情況不同而調整,但其應具有2個特征:需要平滑的數據點權重最大,并且對擬合影響最大;區間外的數據點權重為零,且對擬合沒有影響。
步驟2:加權最小二乘回歸,求得x的平滑值。
步驟3:計算上述平滑過程中殘差,在范圍內計算每個數據點的魯棒權重。權重由bi-square函數給出,如式(6)所示:
(6)
式中:ri為通過平滑過程生成的第i個數據點的殘差;MAD=median(|r|),是數據點與樣本中位數偏差的絕對值的中位數。如果ri<6MAD,則魯棒權重接近1;如果ri>6MAD,則魯棒權重為0,并且相關數據點從平滑計算中排除。
步驟4:使用魯棒權重使得數據再次平滑,利用局部回歸權重和魯棒權重兩者來計算最終的平滑值。
步驟5:重復步驟3、步驟4,共迭代5次。
1.3.1 算法原理
Savitzky-Golay(以下簡稱“S-G”)平滑去噪是時域內低通濾波預處理算法,在滑動平均法的基礎上改進而成,其基本思想是利用多項式卷積濾波系數,如以P點為中心,取鄰近的n個點做多項式擬合,利用擬合得到的多項式求點P的平滑值P1,之后將濾波窗口移動1個樣本單位,如此,將所有數據依次遍歷[16]。
1.3.2 計算流程
Savitzky-Golay算法濾波效果與選取的窗口寬度有關,算法的關鍵在于矩陣算子的求解。設x(n)中的1組數據為x(i)(濾波窗口),i=(-m,-m+1,…-1,0,1,…,m-1,m),濾波窗口寬度為n=2m+1,如將窗口內的數據進行k-1次多項式擬合,即可得到n個k元線性方程組y(i),如式(7)所示:
(7)
式中:y-m,y-m-1,…,ym代表擬合多項式的結果;1,-m,…,(-m)k-1中m代表擬合多項式中的未知數;k-1代表擬合多項式的最高次數;a0,a1,…,ak-1代表擬合多項式的系數;e-m,e-m-1,…,em代表擬合值與實際值之間的誤差。
將式(7)用矩陣形式可表示為式(8):
Y(2m+1)×1=X(2m+1)×k·AK×1+E(2m+1)×1
(8)
(9)
(10)
式中:B=X·(XT·X)-1·XT
為驗證本文提出的3種算法在處理風速傳感器監測數據的適用特征,實驗對某一連續時段內300個風速組成的數據集D(見表1)進行處理。
表1 煤礦井下某一時段風速傳感器監測數據
使用FCM算法時,預先設定聚類中心數c=3,分別表示風速正常波動、風速過低和風速過高3種情況。將風速傳感器監測到的數據進行100次迭代,所得結果如圖1所示。
圖1 FCM算法聚類圖
由圖1和圖2可知,FCM算法將風速傳感器監測到的數據按照預設種類依次分為了3類,樣本中心分別為:2.423 618,0.573 616,2.644 809。巷道風速具有一定的容差性,一段時間內,風速處于動態穩定狀態。雖然第1類和第3類聚類中心不同,但根據現場經驗,這兩類數據性質類似。編號16,26,47,72,107,134,274,284和296共9個樣本數據有微弱的突起,因此,并不能將第1類集合中的數據全部判定為風速正常集合,也不能將第3類集合中的樣本全部判定為風速過大。第2類隸屬度圖中的數據為風速異常數據集,算法對樣本180~208聚集性噪聲去除效果較為明顯。
圖2 FCM算法隸屬度
Rloess算法降噪處理前,需要優選窗口寬度,窗口寬度的確定與數據均方誤差有關。因為樣本數據中含有離群值,利用Rloess算法計算的預測值與樣本數據求均方誤差來優選參數是不科學的。因此,使用MAD法去除離群值后的均方誤差來優選Rloess的窗寬。
圖3為窗框寬度3~14條件下,風速監測樣本均方誤差。當窗口寬度為7時,誤差最小,因此,窗口長度設為7。使用二階多項式回歸,采用三角函數作為范圍內數據點的權值函數,采用6倍的中值絕對偏差MAD數進行魯棒權重分配,得到的平滑結果如圖4所示。
圖3 窗寬與均方誤差關系
圖4 Rloess回歸降噪結果
Rloess算法整體上對樣本進行了平滑降噪處理。與FCM聚類降噪算法不同,編號為16,26,47,72,107,134,274,284和296共9個樣本數據變化較為突出,Rloess算法將其直接剔除掉,噪聲去除后,并沒有對鄰近點數據的平滑產生影響。141號樣本原屬于常規噪聲,但算法將其識別為離群值,在計算中被剔除。Rloess算法對聚集性噪聲樣本180~208沒能作出識別,而是在內部進行了噪聲優化處理。樣本數據209~250小幅震蕩波動,Rloess算法對數據進行了較強的平滑處理,將其擬合成了1條光滑的曲線,并基本保持了原數據的變化趨勢。
采用S-G濾波器對監測數據進行降噪處理。由于風速監測樣本量較大,為盡可能減少數據失真,滑動窗口寬度不宜過大。理論上,階數的取值范圍是從(0,n-1),如圖5所示,窗寬從5逐次遞增到25時,均方誤差整體呈增大趨勢,窗寬為7的4階擬合比窗寬為5的2階擬合的計算結果僅僅降低了0.5%。需要降噪的數據樣本風速浮動不大,考慮到計算速度,最終采用窗寬為5的二次多項式對數據進行擬合,最終結果如圖6所示。
圖5 1~4階相關窗寬與均方誤差計算序列
圖6 S-G算法去噪結果
S-G算法對風速樣本數據整體上進行了平滑降噪處理。編號為16,26,47,72,107,134,274,284和296等9個離群樣本數據參與了平滑,并對鄰近點數據的平滑效果產生了一定影響。與Rloess降噪算法相同,對集群噪聲樣本180~208沒能做出識別,而是在內部進行了噪聲優化處理。對于反復波動樣本數據209~250,S-G算法進行了一定的平滑處理,最大程度地保持了原數據的特性。
短時間內,影響井下風流穩定性的變量可分為過程變量和狀態變量2類。過程變量改變一般由臨近風路發生改變、瓦斯突出、礦井突水、巷道不可逆變形或損壞等因素引起,風速數據體現為聚集性、呈現周期性變化規律,一般持續一定的時間。狀態變量改變一般由人車行駛、采煤機切割、罐籠提升等因素引起,對風速樣本影響相對短時,3種算法對于井下過程變量和狀態變量變化引起的數據噪聲處理效果不盡相同,具體見表2。
表2 3種方法優缺點對比
2.5.1 FCM聚類算法效果分析
FCM對過程變量變化引起的聚集性噪聲處理較為優越。在對風速樣本處理前根據巷道斷面風速分布特點及現場通風管理經驗,指定數據分類數。根據不同的聚類中心,通過隸屬度以及各狀態間銜接和離散的情況來判斷和識別風速異常。
FCM聚類降噪算法對人員流動、車輛駛出等引起的偶發性噪聲識別能力較差,僅能從隸屬度角度分析出此時風速值出現了異常波動,但波動趨勢卻無法合理解釋。此外當監測數據量較大、風速正常波動范圍較大時,預先給定的數據分類數較難確定,此種情況下將很難達到滿意的降噪結果。
2.5.2 Rloess與S-G算法效果對比
Rloess與S-G算法對狀態變量變化引起的數據噪聲處理較為優越,但二者側重點不同。Rloess算法側重于對于偶發性數據噪聲的處理,從第16,26,47,72,107,134,274,284和296樣本處可以看出,由于Rloess算法引入了魯棒性的過程,Rloess算法將這9個離群值剔除處理,且會對局部范圍內參與擬合的樣本利用權值函數進行權值分配,很大程度上降低了某些非必要風速狀態變量對客觀分析造成的影響。S-G算法側重于對趨勢性數據噪聲的處理,在第16,26,47,72,107,134,274,284和296樣本處,S-G法基于時間域上的多項式擬合,對上述9個樣本數據進行了平滑處理。平滑處理過程中需將所有數據遍歷,所有數據會百分之百參與運算,并對最后的降噪結果產生明顯影響,因此S-G算法在去除噪聲的同時可以很好地保持原樣本的形狀。
Rloess與S-G算法均從局部回歸分析出發,平滑過程依賴于周邊數據,往往忽略了全局意義,因此,處理過程變量引起的聚集性噪聲較為乏力,直觀表現為在異常過程變量中去除異常狀態變量,并不能很好地處理異常過程變量。二者對于噪聲的處理較大地依賴于參數的選取,Rloess在剔除離群值的同時,會錯誤的剔除常規值;而S-G方法易受離群值的影響,使得平滑后的數據失真。
2.5.3 3種方法適用場景
針對3種降噪方法的特性,給出它們的適用場景。當礦井發生風流短路、礦井涌水、巷道變形等時,通風狀態會發生改變并維持一段時間,傳感器的風速數據表現出明顯的分類特征,利用FCM降噪算法可以將此類聚集性噪聲識別出來。聚集性噪聲往往不是隨機出現,其背后存在著一定的必然性,還需要現場技術人員根據所收集樣本的聚類結果有針對性地排查通風設施、現場風質或傳感器設備運行狀態。
由于生產需求,巷道中不可避免地存在車輛行駛和人員流動,此時傳感器監測到的通風狀態會發生暫時變化,風速數據會產生較大地波動即產生離群值,如果降噪時離群值參與運算,其平滑結果將會失真,并對通風狀態的識別和實時解算造成影響。此時需要利用Rloess降噪方法對離群值剔除處理后再進行平滑降噪。
當采煤機在工作面割煤、罐籠提升時,通風狀態發生了趨勢性的改變,平滑后的數據要盡量保持原有的數據特征,便于判斷井下工作狀態、排查設備隱患,此時使用S-G算法可以較好地處理此類數據的噪聲。
1)處理由過程變量變化引起的風速異常數據時選用FCM較為優越。使用該方法時,要在分析現場風速波動范圍及引起風速異常原因基礎上,給定合理的聚類中心數目。
2)處理由狀態變量變化引起的風速異常數據時,選用Rloess算法或S-G算法較為合理。S-G算法側重于保持通風數據的特性但容易受異常風速的影響,Rloess側重于去除異常風速,但有時會將正常數據錯誤識別為異常數據。經過測試Rloess算法取窗寬為7時誤差最小,而S-G算法取窗寬為5擬合階數為2階時誤差最小。
3)由過程變量和狀態變量同時引起的風速異常數據處理時,可結合使用FCM-Rloess或FCM-SG算法。首先剔除掉過程變量噪聲,之后再根據數據精度、特點及平衡效果的要求,選用Rloess或S-G算法進行局部優化。