基于優(yōu)化GM（1，N）模型的油氣管道腐蝕速率預(yù)測(cè)

趙年峰，夏艷權(quán)，崔凱強(qiáng)，孔淑穎，解 靜，梁昌晶

（1.中國(guó)石油集團(tuán)渤海石油裝備制造有限公司石油機(jī)械廠，河北 任丘 062552；2.河北華北石油工程建設(shè)有限公司，河北 任丘 062552；3.國(guó)家管網(wǎng)集團(tuán)西南管道有限責(zé)任公司 南寧輸油氣分公司，南寧 530200；4.中國(guó)石油長(zhǎng)慶油田分公司伴生氣綜合利用項(xiàng)目部，西安 710016）

腐蝕是油氣管道失效的主要原因之一，管道一旦發(fā)生穿孔泄漏事故，將對(duì)油氣田企業(yè)和周圍環(huán)境造成重大影響。引發(fā)腐蝕的因素眾多，因素之間具有多重非線性特征，科學(xué)有效地預(yù)測(cè)腐蝕發(fā)展趨勢(shì)，及時(shí)采取適當(dāng)?shù)姆栏胧?，?duì)于一定程度上降低管道失效帶來(lái)的風(fēng)險(xiǎn)具有重要意義[1-2]。目前，有關(guān)管道腐蝕預(yù)測(cè)的模型有回歸分析模型、神經(jīng)網(wǎng)絡(luò)模型和灰色模型（GM）等[3-4]，其中回歸分析模型的精度較低，對(duì)于非線性關(guān)系擬合效果較差；神經(jīng)網(wǎng)絡(luò)模型需要大量的數(shù)據(jù)支撐，對(duì)于小樣本數(shù)據(jù)回歸精度不高，經(jīng)常出現(xiàn)欠擬合現(xiàn)象；灰色模型可用于少數(shù)據(jù)、貧信息的數(shù)據(jù)預(yù)測(cè)，適用于油氣管道腐蝕這類樣本規(guī)律不明顯但系統(tǒng)內(nèi)部發(fā)展明確的預(yù)測(cè)。當(dāng)前灰色模型多以GM（1，1）模型為主，例如，李昊燃等[5]從無(wú)偏灰色、新陳代謝等方面對(duì)GM（1，1）模型進(jìn)行了改進(jìn)；王慶鋒[6]將GM（1，1）模型與馬爾科夫鏈相結(jié)合，通過(guò)狀態(tài)轉(zhuǎn)移矩陣，預(yù)測(cè)下一時(shí)間的腐蝕狀態(tài)；張新生等[7]對(duì)GM（1，1）模型的尾端殘差進(jìn)行了修正，改變了預(yù)測(cè)趨勢(shì)的慣性作用。但以上研究未能反映腐蝕因素對(duì)腐蝕速率的影響，模型具有一定的局限性。

基于此，本研究采用斯皮爾曼相關(guān)系數(shù)尋找腐蝕因素中相關(guān)性較高的變量，通過(guò)隨機(jī)森林（RF）對(duì)變量重要性進(jìn)行排序，去除相關(guān)性較高但重要性較低的變量，將優(yōu)化后的變量輸入GM（1，N）模型，采用人工蜂群算法（ABC）對(duì)GM（1，N）模型的背景值進(jìn)行動(dòng)態(tài)優(yōu)化，最后對(duì)預(yù)測(cè)值和實(shí)際值的差異進(jìn)行檢驗(yàn)分析。研究結(jié)果可為管道防腐措施的制定提供技術(shù)指導(dǎo)。

1 數(shù)學(xué)模型

1.1 斯皮爾曼相關(guān)系數(shù)

斯皮爾曼相關(guān)系數(shù)通過(guò)兩個(gè)變量間的等級(jí)關(guān)系確定相關(guān)關(guān)系，不受樣本數(shù)據(jù)分布和樣本容量影響，且對(duì)于缺失值和極值分布的存在具有較好的適應(yīng)性，斯皮爾曼相關(guān)系數(shù)公式為

式中：rs——斯皮爾曼相關(guān)系數(shù)；

n——樣本容量；

di——兩變量間的等級(jí)差值。

1.2 RF重要度評(píng)價(jià)

RF 在進(jìn)行數(shù)據(jù)選取訓(xùn)練決策樹(shù)時(shí)，對(duì)輸入樣本進(jìn)行了多次有放回的抽樣，但總有一些樣本無(wú)法參與決策樹(shù)訓(xùn)練，這部分樣本稱為袋外數(shù)據(jù)[8]。RF 進(jìn)行重要度評(píng)價(jià)的原理是計(jì)算袋外數(shù)據(jù)在每棵決策樹(shù)上的貢獻(xiàn)度，然后取平均值，即

式中：I——重要度指標(biāo)；

M——決策樹(shù)個(gè)數(shù)；

Error2i——第i個(gè)決策樹(shù)的干擾信號(hào)；

Error1i——第i個(gè)決策樹(shù)的袋外數(shù)據(jù)誤差，當(dāng)加入干擾信號(hào)后，I值下降幅度較大，則說(shuō)明該變量的重要程度較高。

1.3 GM（1，N）模型

GM（1，N）模型是多維灰色模型，其中“1”為1 階，“N”為變量個(gè)數(shù)，包括1 個(gè)自變量和n-1個(gè)因變量[9-10]。

設(shè)x1(0)為自變量序列，即x1(0)=(x1(0)(1)，x1(0)(2)，…，x1(0)(n))，因變量序列為

對(duì)原始序列xi(0)進(jìn)行一階累加形成xi(1)，即

其中，i=1，2，…，n；l=1，2，…，k。

隨后生成x1(1)緊鄰均值序列z1(1)，即

建立微分方程

該式即為GM（1，N）模型，a為發(fā)展系數(shù)，bi為驅(qū)動(dòng)系數(shù)。通過(guò)最小二乘法估計(jì)a和bi的參數(shù)值，并將數(shù)據(jù)重新進(jìn)行一階累減還原得到序列的預(yù)測(cè)值。

1.4 模型評(píng)價(jià)

為了確定模型預(yù)測(cè)結(jié)果的可靠性，采用相對(duì)誤差檢驗(yàn)法評(píng)價(jià)模型精度。計(jì)算殘差、相對(duì)誤差和平均相對(duì)誤差等三個(gè)指標(biāo)，評(píng)價(jià)標(biāo)準(zhǔn)見(jiàn)表1。

表1 模型精度檢驗(yàn)評(píng)價(jià)標(biāo)準(zhǔn)

2 GM（1，N）模型優(yōu)化

經(jīng)研究表明，GM（1，N）模型的誤差關(guān)鍵在于公式（4）中背景值的構(gòu)造存在缺陷，公式（4）是用相鄰空間步長(zhǎng)對(duì)應(yīng)的梯形面積（黑色實(shí)線與底邊圍成的面積）近似代替積分區(qū)域面積（紅色實(shí)線與底邊圍成的面積），如圖1所示。當(dāng)處理某些劇烈變化的數(shù)據(jù)序列時(shí)，誤差較大。基于此，采用動(dòng)態(tài)背景系數(shù)λ（i）替代公式（4）中的“0.5”，通過(guò)動(dòng)態(tài)調(diào)整各區(qū)間的背景系數(shù)，最大程度降低背景值誤差，即

圖1 背景值誤差來(lái)源

對(duì)于λ（i）的選取原則是要使自變量預(yù)測(cè)值與實(shí)際值的誤差最小，等同于多峰多維函數(shù)尋優(yōu)問(wèn)題。在此，采用人工蜂群算法（ABC）進(jìn)行求解[11-12]，算法步驟如下：

（1）初始化種群。根據(jù)公式（7）產(chǎn)生蜜源的初始位置，即目標(biāo)函數(shù)的可行域，根據(jù)公式（8）計(jì)算各蜜源的初始適應(yīng)度值。

式中：xij——第j維空間的第i個(gè)解；

xminj、xmaxj——第j維空間解的下限和上限；

fit（i）——第i個(gè)解的適應(yīng)度值；

f（i）——第i個(gè)解的目標(biāo)函數(shù)值。

（2）采蜜蜂根據(jù)公式（9）在蜜源附近進(jìn)行局部搜索，形成新的位置vij，并根據(jù)公式（8）更新各蜜源的適應(yīng)度值，將更新后的適應(yīng)度值與之前的進(jìn)行比較，并根據(jù)貪婪選擇策略保留較好的解。

式中：Φij——（-1，1）的正態(tài)隨機(jī)數(shù)。

（3）觀察蜂在蜂房中等待采蜜蜂的蜜源信息，然后根據(jù)概率篩選蜜源。

（4）當(dāng)采蜜蜂和觀察蜂對(duì)搜索空間內(nèi)的解均搜索完畢后，如一個(gè)蜜源的適應(yīng)度值還沒(méi)有提高，則放棄該蜜源，采蜜蜂變?yōu)閭刹榉?，重新尋找新蜜源?/p>

（5）迭代次數(shù)達(dá)到設(shè)定的最大迭代次數(shù)時(shí)，算法停止，否則轉(zhuǎn)移至步驟（2）循環(huán)執(zhí)行。

3 實(shí)例分析

以某Φ355.6 mm×6.3 mm 埋地輸氣管道為例，設(shè)計(jì)壓力3 MPa，運(yùn)行壓力2.5 MPa。對(duì)管道沿線某處采用實(shí)地掛片的方式監(jiān)測(cè)腐蝕速率，通過(guò)掛片與管道之間進(jìn)行電性連接模擬防腐層破損點(diǎn)，保證該點(diǎn)的陰極保護(hù)電位和管道一致。在相同埋深處埋置12 組掛片，每組3 個(gè)，每30 天取一組掛片數(shù)據(jù)，序號(hào)1為第1天監(jiān)測(cè)數(shù)據(jù)，序號(hào)2 為第30 天監(jiān)測(cè)數(shù)據(jù)，以此類推，計(jì)算每組掛片的平均腐蝕速率，結(jié)果見(jiàn)表2。

表2 管道腐蝕速率和相關(guān)影響因素

3.1 數(shù)據(jù)預(yù)處理

采用公式（1）計(jì)算表2 中影響因素的斯皮爾曼相關(guān)系數(shù)，結(jié)果見(jiàn)表3。篩選相關(guān)系數(shù)大于0.55 的變量，含鹽量與氯離子質(zhì)量分?jǐn)?shù)的相關(guān)系數(shù)、pH 值和氧化還原電位的相關(guān)系數(shù)均超過(guò)了0.55。

表3 斯皮爾曼相關(guān)系數(shù)矩陣

采用極差法對(duì)表2數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理，在最大決策樹(shù)100的條件下，考察不同葉節(jié)點(diǎn)個(gè)數(shù)對(duì)結(jié)果均方誤差的影響，結(jié)果如圖2所示。當(dāng)葉節(jié)點(diǎn)個(gè)數(shù)為5時(shí)，模型的均方誤差最低，當(dāng)決策樹(shù)個(gè)數(shù)為80時(shí)，均方誤差基本不再變化。因此，確定模型中葉節(jié)點(diǎn)個(gè)數(shù)5，決策樹(shù)個(gè)數(shù)80。

圖2 RF模型參數(shù)優(yōu)選

按照上述模型對(duì)變量進(jìn)行重要度評(píng)價(jià)，結(jié)果如圖3 所示。由圖3 可知，變量重要度排序從大到小依次為土壤電阻率、含鹽量、氧化還原電位、含水量、pH 值、氯離子質(zhì)量分?jǐn)?shù)、硫化物質(zhì)量分?jǐn)?shù)、陰極保護(hù)率和雜散電流，根據(jù)之前的相關(guān)系數(shù)分析，刪除相關(guān)性較高但重要度較低的變量，故剔除氯離子質(zhì)量分?jǐn)?shù)和pH 值，這兩個(gè)變量可以間接的通過(guò)含鹽量和氧化還原電位反映。此外，陰極保護(hù)率和雜散電流的重要性較小，且表1中的陰極保護(hù)率均超過(guò)了90%，雜散電流均小于2.5 mV/m，說(shuō)明管道在土壤中的最小保護(hù)電位基本滿足規(guī)范要求，陰極保護(hù)失效對(duì)腐蝕的影響作用不大，管道周邊不存在直流電氣化鐵路、有軌電車、直流電焊機(jī)等直流接地裝置。最終確定GM（1，N）模型的因變量為土壤電阻率、含鹽量、氧化還原電位、含水量和硫化物質(zhì)量分?jǐn)?shù)。

圖3 變量重要度排序

3.2 訓(xùn)練精度檢驗(yàn)

根據(jù)前文分析結(jié)果，N 取6，建立GM（1，6）模型，為了驗(yàn)證變量篩選的必要性，同時(shí)建立GM（1，10）模型（即未進(jìn)行變量篩選建立的GM模型）。取前8組數(shù)據(jù)建立預(yù)測(cè)模型，通過(guò)后4組數(shù)據(jù)進(jìn)行模型評(píng)價(jià)，參數(shù)列的計(jì)算結(jié)果見(jiàn)表4。

表4 參數(shù)列計(jì)算結(jié)果

在計(jì)算動(dòng)態(tài)背景系數(shù)λ（i）時(shí)，蜂群數(shù)量為20，蜜源數(shù)量為10，搜索限制為20 次，最大迭代次數(shù)為100次，此時(shí)待優(yōu)化問(wèn)題的目標(biāo)函數(shù)為10 維。為保證算法準(zhǔn)確性，反復(fù)運(yùn)行10 次，取誤差最小的一次作為最終結(jié)果，算法在迭代至50次左右收斂，迭代過(guò)程如圖4所示。最終得到優(yōu)化GM（1，6）的λ（i）為0.214 7、0.014 5、0.058 9、0.075 2、0.165 4、0.987 5、0.069 52。

GM（1，10）GM（1，6）優(yōu)化GM（1，6）序號(hào)實(shí)際值1 2 3 4 5 6 7 8 0.021 0.035 0.042 0.058 0.067 0.076 0.085 0.088模擬值0.021 0 0.037 4 0.031 8 0.041 5 0.069 2 0.089 1 0.092 5 0.093 4殘差0-0.002 4 0.010 2 0.016 5-0.002 2-0.013 1-0.007 5-0.005 4相對(duì)誤差/%0 6.857 14 24.285 71 28.448 28 3.283 58 17.236 84 8.823 53 6.136 36模擬值0.021 0 0.036 2 0.041 1 0.059 6 0.068 5 0.074 2 0.081 3 0.090 5殘差0-0.001 2 0.000 9-0.001 6-0.001 5 0.001 8 0.003 7-0.002 5相對(duì)誤差/%0 3.428 57 2.142 86 2.758 62 2.238 81 2.368 42 4.352 94 2.840 91模擬值0.021 0 0.035 8 0.042 9 0.057 4 0.065 7 0.074 2 0.083 5 0.087 4殘差0-0.000 8-0.000 9 0.000 6 0.001 3 0.001 8 0.001 5 0.000 6相對(duì)誤差/%0 2.285 71 2.142 86 1.034 48 1.940 30 2.368 42 1.764 71 0.681 82

圖4 ABC算法迭代過(guò)程

三種模型腐蝕速率訓(xùn)練精度對(duì)比結(jié)果見(jiàn)表5。GM（1，10）模型除了第5 個(gè)樣本相對(duì)誤差較小，其余均較大，平均相對(duì)誤差為11.88%，對(duì)照表1為“勉強(qiáng)”等級(jí)。GM（1，6）和優(yōu)化GM（1，6）模型訓(xùn)練階段的平均相對(duì)誤差分別為2.51%和1.52%，兩者均達(dá)到了表1中“好”等級(jí)。

表5 腐蝕速率訓(xùn)練值與誤差分析

同理，將GM模型中N按照?qǐng)D3的排序依次增加氯離子質(zhì)量分?jǐn)?shù)和pH 值，構(gòu)成GM（1，7）和GM（1，8）模型，統(tǒng)計(jì)訓(xùn)練樣本的平均相對(duì)誤差，結(jié)果如圖5所示。由圖5可見(jiàn)，增加變量后，平均相對(duì)誤差增大，說(shuō)明了變量間的冗余性對(duì)預(yù)測(cè)精度有較大影響，采用斯皮爾曼相關(guān)系數(shù)和隨機(jī)森林對(duì)變量進(jìn)行篩選具有必要性。優(yōu)化GM（1，6）模型的訓(xùn)練精度最高，可以進(jìn)行下一步的預(yù)測(cè)。

圖5 不同模型的平均相對(duì)誤差對(duì)比

3.3 預(yù)測(cè)精度檢驗(yàn)

根據(jù)上述訓(xùn)練模型，對(duì)后4組數(shù)據(jù)進(jìn)行預(yù)測(cè)和模型評(píng)價(jià)，結(jié)果見(jiàn)表6和圖6。其中，GM（1，10）模型預(yù)測(cè)值的相對(duì)誤差均在10%以上，且圍繞實(shí)際值上下波動(dòng)，說(shuō)明影響腐蝕速率的因變量中存在較多的相關(guān)冗余信息，這些信息使模型擬合困難，導(dǎo)致誤差較大；GM（1，6）模型在預(yù)測(cè)階段的相對(duì)誤差和數(shù)據(jù)波動(dòng)逐漸增大，說(shuō)明該模型不適合長(zhǎng)期預(yù)測(cè)，當(dāng)腐蝕速率發(fā)生較大變化時(shí)，下一時(shí)刻預(yù)測(cè)結(jié)果往往是前一時(shí)刻結(jié)果的延續(xù)；優(yōu)化GM（1，6）模型的相對(duì)誤差大幅降低，預(yù)測(cè)值與實(shí)際值吻合較好，可以根據(jù)訓(xùn)練時(shí)建立的非線性關(guān)系及時(shí)對(duì)預(yù)測(cè)結(jié)果進(jìn)行調(diào)整。GM（1，10）模型、GM（1，6）和優(yōu)化GM（1，6）模型預(yù)測(cè)階段的平均相對(duì)誤差分別為15.61%、9.74%和2.03%，分別對(duì)應(yīng)“勉強(qiáng)”、“合格”和“好”等級(jí)，說(shuō)明優(yōu)化后的GM（1，6）模型在訓(xùn)練樣本和預(yù)測(cè)樣本上的表現(xiàn)均較好，不存在欠擬合的現(xiàn)象。

圖6 不同模型預(yù)測(cè)結(jié)果對(duì)比

表6 腐蝕速率預(yù)測(cè)值與誤差分析

4 結(jié) 論

（1）針對(duì)單變量模型預(yù)測(cè)管道腐蝕速率方面的局限性，建立了GM（1，N）模型，通過(guò)斯皮爾曼相關(guān)系數(shù)和RF 算法對(duì)因變量進(jìn)行篩選，降低因變量相關(guān)性對(duì)自變量預(yù)測(cè)結(jié)果的影響，最終選取5個(gè)因素作為因變量，變量重要度排序從大到小依次為土壤電阻率、含鹽量、氧化還原電位、含水量和硫化物質(zhì)量分?jǐn)?shù)。

（2）針對(duì)GM（1，N）模型背景值構(gòu)造存在的缺陷，引入動(dòng)態(tài)背景系數(shù)，采用人工蜂群算法對(duì)最佳背景系數(shù)進(jìn)行尋優(yōu)，實(shí)現(xiàn)對(duì)腐蝕速率的準(zhǔn)確預(yù)測(cè)。

（3）相比GM（1，10）模型和GM（1，6）模型，優(yōu)化GM（1，6）模型在訓(xùn)練階段和預(yù)測(cè)階段的平均相對(duì)誤差大幅降低，分別為1.52%和2.03%，預(yù)測(cè)結(jié)果更接近實(shí)際值，可為油氣管道重點(diǎn)部位的腐蝕防護(hù)提供實(shí)際參考。