金宏偉，徐云柯，邵建宇，王麗娜，苑志猛

（1.浙江浙能臺(tái)州第二發(fā)電有限責(zé)任公司，浙江 臺(tái)州 317109；2.浙江浙能電力股份有限公司，浙江 杭州 310063；3.山東魯軟數(shù)字科技有限公司，山東 濟(jì)南 250001）

0 引言

汽輪機(jī)作為發(fā)電機(jī)組將蒸汽熱能轉(zhuǎn)換為機(jī)械能的核心設(shè)備，其設(shè)備效率直接影響機(jī)組的能耗水平［1-3］。而隨著其設(shè)備運(yùn)行時(shí)間的增長(zhǎng)，其各項(xiàng)性能指標(biāo)必將逐步劣化。通過(guò)分析歷年運(yùn)行數(shù)據(jù)，掌握汽輪機(jī)缸效率的劣化趨勢(shì)，對(duì)電廠的安全運(yùn)行和節(jié)能降耗有著重要意義。傳統(tǒng)方法通常需要進(jìn)行大量試驗(yàn)獲取高壓缸性能數(shù)據(jù)分析劣化趨勢(shì)，試驗(yàn)周期較長(zhǎng)，成本代價(jià)大，且試驗(yàn)數(shù)據(jù)往往不能滿足全工況分析要求。

20 世紀(jì)90 年代以來(lái)，通過(guò)分析大數(shù)據(jù)逐步成為工業(yè)優(yōu)化領(lǐng)域方向重要的方案［4-5］。文獻(xiàn)［6］使用ARMA 模型預(yù)測(cè)方法，用歷史數(shù)據(jù)預(yù)測(cè)得到的預(yù)測(cè)值與實(shí)際值做對(duì)比得到劣化數(shù)值，以此方法得出的劣化數(shù)值與預(yù)測(cè)模型的精確度有較大關(guān)系，難以保證其結(jié)果準(zhǔn)確性，文獻(xiàn)［7］提出了利用經(jīng)驗(yàn)?zāi)B(tài)分解方法與特征流通面積去計(jì)算劣化數(shù)值，而這種方法無(wú)法對(duì)后期微小劣化算得結(jié)果。

從電廠實(shí)際運(yùn)行數(shù)據(jù)入手，采用一種完全基于歷史數(shù)據(jù)的數(shù)值分析方案，通過(guò)對(duì)歷史數(shù)據(jù)進(jìn)行處理，以理論經(jīng)驗(yàn)為基礎(chǔ)設(shè)置工況判別條件，結(jié)合大數(shù)據(jù)分析手段挖掘并處理，在大量雜亂的歷史數(shù)據(jù)中篩選出真實(shí)有效且符合理論依據(jù)的缸效劣化結(jié)論。

1 算法原理

1.1 最小二乘法矩陣運(yùn)算

最小二乘法是一種數(shù)學(xué)優(yōu)化技術(shù)。它通過(guò)最小化誤差的平方和尋找數(shù)據(jù)的最佳函數(shù)匹配，由Adrien-Marie Legendre 在19世紀(jì)提出［8］，E(w)=∑(y-wx)2，其中y表示多組樣本的觀測(cè)值，wx為假設(shè)擬合函數(shù)的理論值，目標(biāo)函數(shù)E(w)即為損失函數(shù)。

假設(shè)函數(shù)y=anxn+a(n-1)x(n-1)+...+a0，其中y為因變量，x為自變量，目標(biāo)為求系數(shù)an、a(n-1)、a0將其構(gòu)建矩陣

目標(biāo)方程為Aw=Y，w=[a0，a1，…，an]T，Y=[y1，y2，…，ym]T，在m=n的情況下有w=A-1Y，在m≠n的情況下有ATAw=ATY，令T=ATA，則 得到（2）式：

1.2 SG濾波法

SG 濾波法（Savitzky Golay Filter）也稱多項(xiàng)式濾波，核心思想為對(duì)窗口內(nèi)數(shù)據(jù)進(jìn)行加權(quán)濾波，本質(zhì)上同是一種基于最小二乘原理的多項(xiàng)式平滑算法的信號(hào)處理方法［9-11］。相比于滑動(dòng)平均法、指數(shù)滑動(dòng)平均法，SG 濾波法在濾波平滑過(guò)程中能夠有效地保留信號(hào)的變化信息。

設(shè)濾波窗口長(zhǎng)度為n=2m+1，使用k階多項(xiàng)式對(duì)窗口長(zhǎng)度內(nèi)數(shù)據(jù)點(diǎn)進(jìn)行擬合，信號(hào)x(t)可以表示為時(shí)間t的函數(shù)

令θ=[a0，a1，…，ak]T為模型的參數(shù)向量，r(t)=[t0，t1，…，tk]T為回歸變量，則得到

在已知n個(gè)這樣的方程后，應(yīng)用最小二乘法確定擬合參數(shù)θ的最優(yōu)估計(jì)

在計(jì)算出矩陣后，即可快速地將觀測(cè)值轉(zhuǎn)為濾波值，某時(shí)刻電廠負(fù)荷數(shù)據(jù)濾波前后效果如圖1。

圖1 SG濾波效果

1.3 MCD異常檢測(cè)算法

最小協(xié)方差行列式（Minimum Covariance Determinant，MCD），是一種常用的多元位置參數(shù)和尺度參數(shù)抗差估計(jì)，而MCD 異常檢測(cè)算法核心為利用均值向量與協(xié)方差陣，結(jié)合馬氏距離計(jì)算思想得到更準(zhǔn)確的離群點(diǎn)探測(cè)，屬于一種魯棒性很強(qiáng)的位置和分布估計(jì)算法［12-15］。

假設(shè)有一組數(shù)據(jù)有n個(gè)樣本，每個(gè)樣本包含p個(gè)元素，構(gòu)成矩陣

經(jīng)典的容忍橢圓定義為一組數(shù)據(jù)x，馬氏距離為

式中：W為權(quán)重函數(shù)；c1為一致性因子；為數(shù)據(jù)集的協(xié)方差矩陣；W最為快捷有效的選擇為

2 具體流程

流程分為10步，如圖2所示，首先對(duì)歷史數(shù)據(jù)進(jìn)行穩(wěn)態(tài)非穩(wěn)態(tài)識(shí)別與篩選，找到機(jī)組運(yùn)行穩(wěn)定階段的有效數(shù)據(jù)；通過(guò)等間距的劃分負(fù)荷與調(diào)閥參數(shù)的距離區(qū)間，將原始數(shù)據(jù)劃分成若干個(gè)區(qū)塊作為工況劃分；分別對(duì)每個(gè)區(qū)塊進(jìn)行包含異常值與小樣本工況剔除的數(shù)據(jù)篩選工作；對(duì)篩選后的數(shù)據(jù)集計(jì)算高壓缸效率；等間距劃分，將數(shù)據(jù)集中的z軸（即高壓缸效率）分為多層；對(duì)調(diào)閥與負(fù)荷兩個(gè)測(cè)點(diǎn)擴(kuò)維，通過(guò)最小二乘法擬合權(quán)重方法進(jìn)行曲面擬合，并對(duì)各層缺失數(shù)據(jù)的區(qū)塊進(jìn)行填補(bǔ)；對(duì)各年各區(qū)間數(shù)據(jù)進(jìn)行統(tǒng)計(jì)與分析歷年高壓缸效率值，根據(jù)歷年高壓缸效率值的變化情況即可得到汽輪機(jī)高壓缸的具體劣化情況。具體步驟如下。

圖2 高壓缸劣化研究流程

1）獲取某年歷史數(shù)據(jù)中與高壓缸效率計(jì)算相關(guān)測(cè)點(diǎn)。

2）通過(guò)對(duì)相關(guān)測(cè)點(diǎn)進(jìn)行穩(wěn)態(tài)識(shí)別，具體步驟包含：分別對(duì)各個(gè)參數(shù)的歷史數(shù)據(jù)進(jìn)行多項(xiàng)式濾波，得到新的數(shù)據(jù)集；對(duì)濾波前的數(shù)據(jù)指定時(shí)間窗口T1，計(jì)算其時(shí)間窗口內(nèi)數(shù)據(jù)的方差，滑動(dòng)時(shí)間窗口，得到方差數(shù)據(jù)集；對(duì)濾波后的數(shù)據(jù)集時(shí)間窗口T1內(nèi)數(shù)據(jù)進(jìn)行一階線性擬合，滑動(dòng)時(shí)間窗口得到的斜率數(shù)據(jù)集；通過(guò)設(shè)定合適的方差與斜率的閾值，在方差數(shù)據(jù)集與斜率數(shù)據(jù)集篩選同時(shí)處于方差與斜率閾值內(nèi)時(shí)間段的數(shù)據(jù)，判定其時(shí)間段內(nèi)的數(shù)據(jù)為穩(wěn)態(tài)工況數(shù)據(jù)。

3）通過(guò)不同的負(fù)荷與調(diào)閥開度對(duì)數(shù)據(jù)進(jìn)行工況劃分，將穩(wěn)態(tài)工況數(shù)據(jù)集分為若干個(gè)區(qū)塊。

4）對(duì)樣本量足夠多的區(qū)塊進(jìn)行保留，若該工況下樣本量低于N條，則舍棄該區(qū)塊下全部數(shù)據(jù)。隨后對(duì)不同區(qū)塊下的數(shù)據(jù)使用MCD算法進(jìn)行異常值剔除。

5）對(duì)步驟2）處理后得到的新數(shù)據(jù)集進(jìn)行逐條計(jì)算高壓缸效率，得到缸效數(shù)據(jù)集。

6）對(duì)不同的區(qū)塊獲取其歷史最大缸效值與歷史最小缸效，按照其差值將其等間距的劃分為10 層。再將不同工況下相同層的數(shù)據(jù)進(jìn)行合并得到同層數(shù)據(jù)集，共10個(gè)。

7）對(duì)同層數(shù)據(jù)表進(jìn)行負(fù)荷、調(diào)閥與缸效之間的進(jìn)行擴(kuò)維后再使用最小二乘法擬合各個(gè)參數(shù)與缸效之間的權(quán)重，共得到10 個(gè)不同的非線性回歸擬合模型。

8）將穩(wěn)態(tài)數(shù)據(jù)中每一條數(shù)據(jù)的負(fù)荷與調(diào)閥分別帶入10 個(gè)非線性回歸擬合模型中得到10 個(gè)推算缸效值，并與業(yè)務(wù)計(jì)算得到的計(jì)算缸效值做差，找到其距離最近（差值絕對(duì)值最?。┑耐扑愀仔е怠２⒕嚯x最近的推算缸效值所對(duì)應(yīng)的模型建立層級(jí)視為該歷史數(shù)據(jù)的層級(jí)。

9）對(duì)所有穩(wěn)態(tài)數(shù)據(jù)都按上一步驟獲取層級(jí)，統(tǒng)計(jì)不同層級(jí)下不同負(fù)荷調(diào)閥工況下的高壓缸效率均值與樣本量。

10）再代入下一年的歷史數(shù)據(jù)，重復(fù)上述步驟，比較兩年歷史數(shù)據(jù)10 層中的均值與樣本量即可得到其高壓缸性能變化數(shù)據(jù)。

3 算例測(cè)試

以華南某電廠百萬(wàn)機(jī)組高壓缸的數(shù)據(jù)對(duì)本文所提方法的有效性進(jìn)行驗(yàn)證，使用2016 年1 月1 日00：00 至2019 年12 月31 日24：00 的數(shù)據(jù)進(jìn)行逐年建模分析。

3.1 數(shù)據(jù)獲取

測(cè)點(diǎn)共選取主蒸汽壓力、大氣壓力、主蒸汽溫度、高壓缸排汽壓力、高壓排汽溫度、機(jī)組負(fù)荷、調(diào)閥開度共計(jì)7項(xiàng)指標(biāo)。

利用廠級(jí)信息監(jiān)控系統(tǒng)選取2016—2019 年共計(jì)4 年的運(yùn)行數(shù)據(jù)，為保證后續(xù)穩(wěn)態(tài)數(shù)據(jù)篩選精度，采用5 s 一條的采樣頻率進(jìn)行插值取數(shù)。每年樣本量均為620 萬(wàn)條。通過(guò)前5 項(xiàng)指標(biāo)主蒸汽壓力、大氣壓力、主蒸汽溫度、高壓缸排汽壓力、高壓排汽溫度通過(guò)業(yè)務(wù)計(jì)算得到每條數(shù)據(jù)的計(jì)算缸效值。

3.2 穩(wěn)態(tài)數(shù)據(jù)篩選

對(duì)4 年數(shù)據(jù)的機(jī)組負(fù)荷與主蒸汽壓力分別進(jìn)行多項(xiàng)式濾波算法，通過(guò)選取某點(diǎn)前后各30 條數(shù)據(jù)進(jìn)行該點(diǎn)數(shù)據(jù)的2階多項(xiàng)式擬合計(jì)算，最終得到4年每年620 萬(wàn)條數(shù)據(jù)的新擴(kuò)維參數(shù)機(jī)組負(fù)荷濾波值與主蒸汽壓力濾波值。

由于在負(fù)荷與主蒸汽溫度穩(wěn)定的情況下缸效在短時(shí)間內(nèi)即可以穩(wěn)定，不會(huì)有很強(qiáng)的滯后性，故這里分別對(duì)機(jī)組負(fù)荷濾波值與主蒸汽壓力濾波值通過(guò)滑動(dòng)時(shí)間窗口10 min，步長(zhǎng)為5 s 的方式逐段進(jìn)行一階線性擬合，將一階線性擬合出的斜率作為該段濾波后數(shù)據(jù)的斜率，得到4年每年620萬(wàn)條數(shù)據(jù)的新擴(kuò)維參數(shù)機(jī)組負(fù)荷斜率與主蒸汽壓力斜率。

分別對(duì)機(jī)組負(fù)荷與主蒸汽壓力進(jìn)行原始參數(shù)的方差計(jì)算，選取與步驟2）滑動(dòng)相同時(shí)間窗口10 min與步長(zhǎng)5 s 的方式逐段計(jì)算方差，得到4 年每年620萬(wàn)條數(shù)據(jù)的新擴(kuò)維參數(shù)機(jī)組負(fù)荷方差與主蒸汽壓力方差。

通過(guò)設(shè)定合適的斜率與方差閾值，這里通過(guò)分別計(jì)算兩個(gè)測(cè)點(diǎn)的斜率與方差頻次分布直方圖。在斜率頻次分布直方圖獲取其拐點(diǎn)附近合適的數(shù)值作為斜率閾值設(shè)定值；在方差頻次分布直方圖獲取其雙峰結(jié)構(gòu)中的波谷點(diǎn)附近合適的數(shù)值作為方差閾值設(shè)定值（前峰為穩(wěn)態(tài)數(shù)據(jù)集中分布峰，后峰為非穩(wěn)態(tài)數(shù)據(jù)分布峰）。2016 年斜率頻次分布如圖3 所示，方差頻次分布如圖4所示，最終斜率閾值選定0.035，方差閾值選定7.6。

圖3 斜率頻次分布

圖4 方差頻次分布

通過(guò)上一步中所得閾值對(duì)數(shù)據(jù)進(jìn)行篩選，閾值內(nèi)的作為其穩(wěn)態(tài)數(shù)據(jù)段，閾值外的作為非穩(wěn)態(tài)數(shù)據(jù)段。

對(duì)所得每段穩(wěn)態(tài)數(shù)據(jù)各剔除首尾各6 條數(shù)據(jù)（30 s）。

最終在穩(wěn)態(tài)篩選完成后，2016 年數(shù)據(jù)留存181萬(wàn)條，2017 年留存194 萬(wàn)條，2018 年留存187 萬(wàn)條，2019年留存202萬(wàn)條。

3.3 各工況下異常數(shù)據(jù)剔除

以2016 年相關(guān)數(shù)據(jù)為例，分別計(jì)算在穩(wěn)態(tài)識(shí)別前調(diào)閥開度、負(fù)荷與高壓缸效率的spearman 相關(guān)系數(shù)，得到負(fù)荷與高壓缸效率的相關(guān)系數(shù)為0.31，調(diào)閥開度與高壓缸效率的相關(guān)系數(shù)為0.76；穩(wěn)態(tài)識(shí)別后負(fù)荷與高壓缸效率的相關(guān)系數(shù)為0.14，調(diào)閥開度與高壓缸效率的相關(guān)系數(shù)為0.93。其余3年的相關(guān)系數(shù)數(shù)值也表明其調(diào)閥開度與高壓缸效率相關(guān)性較強(qiáng)，負(fù)荷與高壓缸效率相關(guān)性較弱。

由步驟1）的結(jié)果得出結(jié)論，穩(wěn)態(tài)條件下，缸效與調(diào)閥開度有絕對(duì)正相關(guān)關(guān)系，與負(fù)荷關(guān)系較弱，因此以調(diào)閥開度每1 個(gè)單位，負(fù)荷每10 個(gè)單位從業(yè)務(wù)邏輯上劃分工況，將數(shù)據(jù)分割區(qū)間。

由于機(jī)組運(yùn)行條件下負(fù)荷范圍為400～1 050 MW，調(diào)閥開度為20%～97%，所以在將穩(wěn)態(tài)數(shù)據(jù)按步驟2）條件分割區(qū)間后，共產(chǎn)生65×77 塊數(shù)據(jù)，分別對(duì)每一塊數(shù)據(jù)運(yùn)用MCD 算法進(jìn)行異常值剔除，以此保證異常值剔除的準(zhǔn)確性。

在異常值數(shù)據(jù)剔除完成后2016 年數(shù)據(jù)留存174萬(wàn)條，2017 年留存186 萬(wàn)條，2018 年留存180 萬(wàn)條，2019 年留存196 萬(wàn)條。剔除異常值后數(shù)據(jù)的三維散點(diǎn)圖如圖5所示。

圖5 剔除異常值后數(shù)據(jù)的三維散點(diǎn)圖

3.4 數(shù)據(jù)分層

分別統(tǒng)計(jì)上述每一個(gè)工況下數(shù)據(jù)量，若數(shù)據(jù)條數(shù)低于500，則將該區(qū)間下的數(shù)據(jù)全部舍棄，未舍棄的數(shù)據(jù)則統(tǒng)計(jì)其最大值與最小值，將其按照等間距方式在最大值與最小值中劃分成十層；再將該工況下的所有數(shù)據(jù)分別放入十層中，以此類推將全部工況下的數(shù)據(jù)放入十層中，以此達(dá)到模糊工況劃分的目的。

3.5 缸效擬合與曲面構(gòu)建

由于在機(jī)組真實(shí)運(yùn)行階段會(huì)出現(xiàn)某些工況的發(fā)生次數(shù)較少或完全不出現(xiàn)，故需要對(duì)這些未發(fā)生的工況或者發(fā)生次數(shù)較少因此在分層階段中舍棄掉的數(shù)據(jù)進(jìn)行推測(cè)補(bǔ)全。

同樣基于前述其2016 年兩個(gè)相關(guān)系數(shù)分別為0.93 與0.14 的原因，假設(shè)其調(diào)閥開度、負(fù)荷與高壓缸缸效有直接關(guān)系，故可以簡(jiǎn)單地用其前兩個(gè)參數(shù)推算高壓缸缸效。以缸效率作為因變量Yk，調(diào)閥開度、機(jī)組負(fù)荷作為自變量，分別計(jì)為a，b進(jìn)行3階擴(kuò)維，得到參數(shù)a，b，ab，a2，b2，a2b，ab2，a3，b3，加上常數(shù)項(xiàng)z的后的權(quán)重分別為w1，w2，w3，w4，w5，w6，w7，w8，w9，w10最終得到

將2016 年各層現(xiàn)有數(shù)據(jù)作為訓(xùn)練數(shù)據(jù)，通過(guò)最小二乘法擬合式中各參數(shù)權(quán)重w1，w2，w3，w4，w5，w6，w7，w8，w9，w10，最終分別得到十層數(shù)據(jù)擬合出的各層公式作為各層的非線性回歸擬合模型。

在填入全部數(shù)據(jù)后，得到2016 年十個(gè)非線性回歸擬合模型的曲面，各層曲面圖形極為相似，但其在高調(diào)閥開度上的，相同調(diào)閥開度與相同負(fù)荷對(duì)應(yīng)下的缸效值會(huì)有明顯差異，以其中第1 層為例曲面圖形如圖6所示。

圖6 2016年第一層數(shù)據(jù)回歸擬合曲面圖

3.6 數(shù)值統(tǒng)計(jì)與比較

將前述數(shù)據(jù)分層步驟前的2019 年全部數(shù)據(jù)逐條利用其機(jī)組負(fù)荷和調(diào)閥開度兩個(gè)參數(shù)代入至基于2016 年數(shù)據(jù)的十個(gè)非線性模型中，得到推測(cè)出的十個(gè)推算缸效值。

將該條數(shù)據(jù)的計(jì)算缸效值與推測(cè)得到的十個(gè)推算缸效值進(jìn)行差值計(jì)算，統(tǒng)計(jì)與其差值絕對(duì)值最小的推算缸效值隸屬于基于第幾層數(shù)據(jù)所建立的非線性回歸預(yù)測(cè)模型。則該條數(shù)據(jù)應(yīng)屬于第幾層。

在將剔除異常值的穩(wěn)態(tài)數(shù)據(jù)全部統(tǒng)計(jì)完后，分別對(duì)每一層不同工況數(shù)據(jù)計(jì)算其均值，則當(dāng)年處理年份的數(shù)據(jù)獲取完畢。

在用同樣的方法將2016—2019 共4 年的數(shù)據(jù)全部統(tǒng)計(jì)完畢后，通過(guò)選取其具有典型性與代表性的幾種運(yùn)行工況，分別得到該工況下每一年數(shù)據(jù)每一層數(shù)據(jù)的均值。

通過(guò)對(duì)4 年均值進(jìn)行逐條對(duì)比即可得到其高壓缸在未發(fā)生異常情況下逐年劣化情況。

展示算例中所選機(jī)組在2019 數(shù)據(jù)各層的均值如表1所示，各年各層平均缸效數(shù)值如表2所示。

表1 2019年各層非線性擬合函數(shù)各項(xiàng)數(shù)值

表2 2016—2019年各層擬合后缸效平均值 單位：%

4 結(jié)語(yǔ)

將大量歷史運(yùn)行數(shù)據(jù)進(jìn)行處理分析，得出了高壓缸性能逐年劣化趨勢(shì)的具體數(shù)據(jù)。機(jī)組在2016 年至2019 年中，基本展現(xiàn)出以每年0.03%～0.15%的劣化趨勢(shì)發(fā)展，如圖7 所示，其中2016 年至2017 年平均劣化0.143%，2017 年至2018 年平均劣化0.06%，2018 年至2019 年平均劣化0.03%。利用該方法不僅能夠做到較為準(zhǔn)確地劣化數(shù)值，用于分析設(shè)備劣化情況，對(duì)機(jī)組的運(yùn)行及優(yōu)化同樣具有著指導(dǎo)意義。以該方法得出的劣化趨勢(shì)數(shù)據(jù)為基礎(chǔ)，分析近期高壓缸性能數(shù)據(jù)亦可進(jìn)行短期設(shè)備風(fēng)險(xiǎn)預(yù)測(cè)，若設(shè)備參數(shù)偏離正常劣化區(qū)間則可判別高壓缸內(nèi)發(fā)生異常情況，為機(jī)組的安全評(píng)估與能效優(yōu)化方面提供幫助。

圖7 2016-2019年各層缸效