岳健國 郭瑞 傅行軍 田新啟 宗曜犇 王旭

（1.大型發(fā)電裝備安全運行與智能測控國家工程研究中心；2.東南大學能源與環(huán)境學院）

0 引言

目前，我國正朝著實現(xiàn)“雙碳”目標大步邁進，構建“以新能源為主體的新型電力系統(tǒng)”是實現(xiàn)“雙碳”目標的重要手段。新型電力系統(tǒng)中高比例可再生能源和高比例電力電子設備接入電力系統(tǒng)的特征會導致電網運行波動性升高，從而對火電機組深度調峰能力提出了更高的要求。同時，火電機組負荷的多變使得各種輔機的運行狀況發(fā)生較大變化，深入研究輔機參數間內在規(guī)律對維持其穩(wěn)定健康運行及保證發(fā)電機組運行效率有重要意義。

當前發(fā)電設備運行狀態(tài)分析方法主要分為兩類，即基于模型的方法和基于數據驅動的方法。其中，基于數據驅動的方法在信息技術和電力設備運行數據采集技術的大力推動下得到快速發(fā)展，此類方法能揭示設備運行參數間的復雜耦合關系，并能對設備運行狀態(tài)和故障進行有效監(jiān)測及預警。

目前數據分析技術主要依靠相關性概念分析設備運行數據間的關聯(lián)關系，在設備強相關關系的基礎上，若能結合參數間因果不對稱性，利用因果推斷技術研究設備運行變量間的因果關系，將能更加深入挖掘設備內在規(guī)律性，為機組輔機安全穩(wěn)定運行提供幫助，并能為現(xiàn)場運行調整提供有力指導。

因果推斷理論已在醫(yī)學、生物制藥、社會實驗、氣候變化等問題的研究上已獲得廣泛應用[1-3]。同時，已有部分將因果推斷技術應用到工業(yè)領域中去的嘗試，穆鋼[4]在信息幾何因果推理方法的基礎上提出一種逆信息熵因果推理的方法用于電力系統(tǒng)運行數據中因果關系分析，Zhang[5]將個體平均因果效應ITE用到了室內空調系統(tǒng)的故障診斷中，通過正常異常情況下數據分類對比結合因果關系判斷故障類型。Sun[6]提出了一種將因果圖和滑動窗口KPC結合的方法用于復雜工業(yè)過程的故障監(jiān)測中。

但現(xiàn)有研究中將因果發(fā)現(xiàn)模型用于實際設備中的實踐較少，其中因果推斷技術中不同因果發(fā)現(xiàn)模型所要求的不同數據生成方式和噪聲假設是其用于實際設備時進行模型選用及準確性和有效性判斷的一大障礙。

本文對因果關系發(fā)現(xiàn)及評估算法進行了改進，并將專家知識引入數據驅動的因果推斷技術中。同時，以某一次風機為研究對象驗證了文中所提方法的有效性，利用本方法構建了其合理可信的因果關系網絡、并通過計算因果效應強度反應了其狀態(tài)惡化過程和狀態(tài)突變節(jié)點。

1 基于因果函數的改進因果推斷方法

1.1 改進因果關系發(fā)現(xiàn)及評估方法

因果關系發(fā)現(xiàn)是因果推斷技術的一部分，目前常見的因果發(fā)現(xiàn)算法主要分為兩類，分別是基于約束的方法和基于因果函數模型的方法，基于約束的方法由于無法解決Markov等難題所以存在部分無法判斷的因果關系方向[7]。而基于因果函數模型的方法通常對數據的生成方式和噪聲有較強的假設性。計算機領域往往有標準集對模型因果關系識別結果有效性和準確度進行驗證，但實際生產設備中不同參數間因果關系通常只能依靠物理模型或隨機試驗進行驗證，而物理模型的復雜性和隨機試驗的高成本或不可實施性使得模型對于實際設備因果關系識別結果的驗證往往無法順利進行。

因此，為了提高從觀測數據識別因果關系的準確度，本文提出一種改進因果關系發(fā)現(xiàn)及評估方法，首先對原始數據進行分組，在分組大小保證所選模型可正常識別的情況下，分別使用所選模型對不同組的數據進行因果關系發(fā)現(xiàn)。據此，使用因果保持率CRR（Causality retention rate）描述同一因果方向在多組數據中被識別出的概率:

式中，CRRi代表第i條因果關系方向的保持率；pi代表第i條因果關系出現(xiàn)的次數；N代表組的個數。

為了校驗所選模型準確度，本文提出設定強因果關系和非因果關系。結合專家知識將已知具有明確物理意義的參數間因果方向設置為強制因果方向，已知絕對不是因果方向的路徑設置為非因果方向。據此判斷每組因果圖識別結果的置信度Si,并使用置信度修改每組數據得到的因果圖中除強因果方向和非因果方向外其他方向在總體因果圖的鄰接矩陣中的權重。因果圖置信度Si計算公式如下:

式中，Si為每組識別得到的因果圖的置信度；m為分組的總數量；hi為錯誤或未識別出的強因果方向的數量；li為錯誤識別的非因果方向的數量；M為強因果方向和非因果方向的總數量。

之后根據Si對因果關系保持率進行修改得:

除此之外，本文提出利用CRRi計算因果信息不對稱度DDAi（Degree of directional asymmetry），用以描述同一因果關系兩方向之間差異，判斷兩參數間因果方向識別結果的可靠性。

式中，CRRi+和CRRi-代表同一關系兩個方向的因果關系保持率。

基于因果函數的模型識別準確率都較為依賴不同模型所提出的數據生成方式假設和兩參數間的因果信息不對稱度。其中，數據生成方式假設可以依靠專家經驗進行辨別，而因果信息不對稱度DDAi可以通過式（4）進行計算。對CRRi和DDAi設置一定的閾值，即可最終確定基本正確的因果關系網絡。完整改進因果關系發(fā)現(xiàn)及評估方法流程圖如圖1所示。

圖1 完整改進因果關系發(fā)現(xiàn)及評估方法流程圖Fig.1 The flow diagram of improved causal relationship discovery and evaluation method

1.2 線性非高斯無環(huán)模型

因果關系發(fā)現(xiàn)作為因果推斷技術的一部分，針對基于約束的方法存在的馬爾可夫等價類問題，許多學者從因果作用機制出發(fā)，以結構因果模型為基礎結合數據分布特性提出了不同的因果函數模型。

線性非高斯無環(huán)模型LiNGAM[8]（Linear Non-Gaussian Acyclic Model）屬于基于因果函數的因果發(fā)現(xiàn)模型，此模型對數據生成方式提出了三點假設:

a）數據生成方式是線性的，變量xi是對應的原因變量xj,噪聲變量ei和常數ci的線性求和，即:

b）模型中不存在未觀測到的混雜因子；

c）噪聲變量ei獨立非高斯，即ei服從非零方差的非高斯分布。

通過數據預處理減去每個變量的均值，消去公式（5）中的常數項ci，從而得到以下矩陣方程:

式中，X為i維隨機變量；B為i×i的連接矩陣；N為i維的非高斯隨機噪聲變量。

對式（6）求解X可得:

式中，A=(I-B)-1，A矩陣可以通過適當的置換算法置換為下三角矩陣或近似為下三角矩陣。

式（7）和噪聲變量ei的獨立非高斯性構成了標準的線性非高斯無環(huán)模型。通過獨立成分分析ICA算法[9]（Independent Component Analysis）等算法可以對LiNGAM 模型進行求解,獲得包含具有所尋求的連接強度bij的矩陣B，從而即可得到具體因果方向，最后通過修剪邊枝算法[10]得到最終因果圖網絡。

1.3 因果效應強度計算

在得到各參數間的因果關系網絡后，因果推斷技術可以量化各因果方向“因”對“果”的影響程度即因果效應強度。對原因X按是否進行了干預分為x=1和x=0組，Y1(x)和Y0(x)分別代表單個個體接受與未接受x干預后的結果。雖然無法同時觀察到單個個體的Y1(x)和Y0(x)從而計算x的作用結果，但是可以識別總體的因果效應，總體的因果效應強度為進行干預與未進行干預時結果差值的期望ATE(average treatment effect)，即:

當數據維度升高時，在兩組間尋找除X外其他參數均相同的個體幾乎是不可能的，因此本文選用傾向得分[11]作為多維數據匹配依據。

傾向得分匹配的原理是，通過以所計算因果關系中的“因”為因變量，以協(xié)變量作為自變量，利用Logistic回歸建立回歸模型估計每個樣本接受實驗的可能性作為傾向性得分，然后，根據傾向得分使用近鄰匹配法在實驗組和對照組間進行匹配。由此，就可以計算具有多維可觀測特征樣本中實驗組收到參數X的干預后的因果效應ATT(x)(average treatment effect for thextreated group)，計算公式為式（9）。

最后，對式（9）所有ATT(x)進行加權平均得到ATT式（10）:

2 一次風機樣本實驗驗證

2.1 一次風機因果發(fā)現(xiàn)模型選擇

本文選取某火電廠一次風機整年運行數據為研究對象，測點包括調節(jié)擋板閥位、風機電流、軸承溫度（前、中、后）、軸承X和Y方向振動信號，根據一次風機物理性質結合專家知識，本文設置了一條強因果關系和六條非因果關系，如圖2所示。

圖2 強因果關系與非因果關系示意圖Fig.2 The diagram of strong causality and non-causality

從原始數據抽取1000 組樣本，每組樣本量10 萬個數據點，分別使用不同因果發(fā)現(xiàn)模型進行識別，由公式（1.3）計算置信度，最終計算結果如表1所示：

表1 不同模型對此一次風機因果圖識別結果置信度Tab.1 The Confidence coefficient of cause-and-effect diagram recognition result of primary air fan with different model

由表1 可見，對此一次風機，線性非高斯無環(huán)模型LiNGAM[8]的識別結果具有極高置信度，因此本文選用LiNGAM模型對一次風機進行因果關系發(fā)現(xiàn)。

2.2 因果函數模型可識別性分析

基于因果函數模型使用前需要驗證參數不對稱性和模型假設。一次風機的參數間固然具有較強的因果信息不對稱性，電流增加會導致?lián)醢彘y門開度增大，但是反過來不成立。擋板閥門開度增大會對前軸承溫度產生影響，雖然不能確定正向因果關系是否成立，但是反方向一定不成立，不對稱性顯著。

對于LiNGAM模型的三條假設:1）在一次風機穩(wěn)定運行階段，各參數間相互作用近似處于線性范圍內，運行數據基本服從線性變化。2）本文選取數據中包含一次風機主要參數測點，不存在重要的未觀測的混雜因子。3）假設風機運行過程中的噪聲項為非高斯分布。

由此可見，一次風機中參數間不對稱性清晰，參數情況滿足LiNGAM模型假設。因此，對一次風機運行變量使用LiNGAM模型進行因果關系分析是可行的。

2.3 構建一次風機參數間因果關系網絡

對此一次風機運行數據重新抽取1000 組樣本，每組樣本量10萬個數據點，經過LiNGAM模型識別，計算每兩個變量間的因果關系保持率CRR。最終結果因果方向保持率CRRi如表2所示。

表2 因果保持率CRRi 矩陣Tab.2 The causal retention rate CRRi matrix

由表2 可以看出，電流指向閥門開度的因果方向CRR為100%，對于六條非因果方向，除了后軸承溫度對擋板閥門開度出現(xiàn)兩次誤判，其他五條在所有1000 次計算中均未被識別，說明此模型成功識別出了事前設定的一組強因果關系和六組非因果關系，表明LiNGAM模型可以識別出一次風機中基于物理規(guī)律的總體框架，并能夠對其他復雜耦合方向進行識別。

因果關系網絡中某一關系具體因果方向的識別結果可靠性主要依賴于各方向間參數的因果不對稱度，在得到各方向間的因果方向保持率之后，通過計算其不對稱度來確定識別結果的可信度。不對稱度越高，因果方向識別結果可信度越高。設置CRRi的閾值為20%，對超過此閾值的因果關系方向計算因果信息不對稱度DDAi。因果信息不對稱度的計算結果如表3 所示，不對稱度DDAi的置信值選取為0.2。之后，對于因果信息不對稱度DDAi≥0.2 的邊，選取CRRi較大的方向作為因果圖中的因果方向，最終得到的因果關系網絡如圖3所示。

表3 因果關系不對稱度Tab.3 The causality asymmetry

因果關系網絡可以為更深層次地認識系統(tǒng)物理過程提供有利幫助。如圖3所示，對于溫度變化和振動的耦合問題，本系統(tǒng)中前中后軸承的溫度變化都是軸承X和Y方向振動改變的“因”。這些因果方向的CRR基本均在50%以上且后軸承溫度對軸承Y向振動的CRR達到了77.15%，說明此一次風機運行過程中溫度變化是振動情況改變的一個原因，為風機的溫度振動耦合作用的解耦指出具體的方向。

圖3 通過LiNGAM模型識別得到的一次風機因果關系網絡Fig.3 The causality network of primary air fan identified through LiNGAM model

2.4 因果效應強度計算

得到因果關系網絡后，“因”對“果”的影響程度即因果效應強度也是因果推斷的關注點之一。本文使用滑動窗口結合傾向得分匹配的方法量化不同因果關系方向之間的因果效應強度隨時間的變化情況。本文選取窗口長度為24小時，滑動距離一小時。

通過對具有因果關系兩參數的值按升高和降低以及不變進行分組處理后計算傾向得分，再根據傾向得分進行匹配。計算參數變化組即實驗組的平均干預效應即為此因果方向的因果關系強度。

為進一步驗證因果關系網絡的準確性，本文通過三組反駁實驗測試2.3 節(jié)得到的因果關系網絡的準確性，分別是隨機添加混雜因子，數據子集驗證和安慰劑干預。雖然無法通過觀測數據對因果關系網絡進行證真，但可以通過反駁測試中模型出現(xiàn)的異常行為對其證偽。若模型結構正確，則隨即添加混雜因子和數據子集測試時的因果效應強度與原始結果應近似相同，安慰劑干預測試是用獨立的隨機變量替換真正的原因變量計算因果效應，若模型正確，其結果應該接近0,三組測試結果如圖4所示。

由三組反駁測試結果圖4可以看出，隨機添加混雜因子和數據子集測試時的因果效應強度均與原始效應差別很小，兩者殘差分別在0.1和0.12以內，說明在此因果圖中隨機添加混雜因子對其無較大影響，且用隨機選擇的子集替換給定的數據集時也具有較高的準確性。安慰劑干預測試的結果表明作為原因的獨立隨機變量對結果的因果效應強度都為0。在三種測試情況下此因果關系網絡均沒有異常產生，進一步證明了此因果關系網絡擁有。

圖4 反駁測試結果Fig.4 The refuted test results

2.5 因果狀態(tài)監(jiān)測

根據2.3 節(jié)得到的因果圖和1.3 節(jié)所述的因果效應強度計算方法，計算得各因果方向間的主要因果效應變化曲線如圖5所示。

圖5 各因果方向因果效應強度隨時間窗口變化曲線Fig.5 The causal effect intensity of each causal direction varies with the time window

由圖5可以看出：（1）如圖5-1所示，預設的強因果關系的效應變化情況符合實際物理模型。此一次風機全年電流-擋板閥位開度的因果效應強度基本一直維持在0.8附近，說明電流大小的變化將在很大程度上影響擋板開度的變化，符合實際物理模型。1月100點處的突降對應為停機，但在7月上旬750點處，此值由0.8突降至0附近，觀察實際數據發(fā)現(xiàn)此時負荷增大，擋板閥門開度開至最大，電流對其作用效果達到飽和，進一步增大電流也不會引起閥門開度的變化，因此因果效應降為0。

（2）如圖5-2 所示，電流大小對軸承X向振動的因果效應變化情況大概經歷了200，750，1600三個窗口節(jié)點，其中，在窗口200 點，系統(tǒng)在長時間停機檢修后投產，此效應緩慢爬升后穩(wěn)定在0.5。在750點窗口即7月2 日附近，此值發(fā)生突變由0.5 升至0.7 高位附近保持，說明此時電流大小變化對軸承X向振動的影響增大。而在此處，軸承X向振動對軸承Y向振動的因果效應強度如圖5-4 也發(fā)生了突變，由原來在0 波動躍升到0.4附近，說明此時兩個方向的振動不再獨立，上述討論說明7 月2 日的高負荷導致系統(tǒng)振動情況發(fā)生突變性惡化。此后，1600點窗口對應11月初的停機檢修，從原始數據及圖5可以看出，此次停機檢修后系統(tǒng)振動情況發(fā)生了較大改善。翻閱檢修報告知，此風機失效軸承潤滑油脂過少，在7 月中旬后逐漸出現(xiàn)軸承振動異常，并于11月1日進行停機檢修替換失效軸承，系統(tǒng)狀態(tài)變化過程與本文結果相符。

此例證明本方法成功跟蹤了系統(tǒng)狀態(tài)的變化，并發(fā)現(xiàn)原始數據無法得到的7 月2 日高負荷運行導致的系統(tǒng)狀態(tài)突變情況。

3 結論

1）提出了一種改進因果關系發(fā)現(xiàn)和評估方法，通過引入因果關系保持率CRR、因果信息不對稱度DDA，并結合專家知識設置強因果關系和非因果關系，可以對不同的因果函數模型因果關系識別準確度給出具體的評估結果，優(yōu)化因果發(fā)現(xiàn)過程，提高因果關系網絡的可靠性。

2）以某一次風機為對象進行因果關系發(fā)現(xiàn)，證明了改進的因果發(fā)現(xiàn)評估方法能成功判斷不同模型的有效性和準確度并對因果關系發(fā)現(xiàn)過程進行修正。CRR評分表明LiNGAM模型對此一次風機強、非因果關系的識別準確率達到了99.95%。三組反駁測試結果均表明了此方法得到的因果關系網絡的正確性。此方法為對發(fā)電設備進行因果關系發(fā)現(xiàn)提供了新的模型評估方法，解決了因果模型難以選擇和準確度無法判斷的問題。

3）因果狀態(tài)監(jiān)測提供了一種發(fā)電機組輔機運行狀態(tài)分析的新角度，本文以某一次風機為例證明其不僅能夠成功跟蹤反映系統(tǒng)狀態(tài)，還能夠發(fā)現(xiàn)原始數據無法發(fā)現(xiàn)的系統(tǒng)狀態(tài)惡化突變節(jié)點從而進行精確高效的故障告警。結合因果推斷技術的數據分析方法可以成為研究系統(tǒng)更深層次本質規(guī)律和進行運行狀態(tài)分析的有效手段。