基于UKF的海上天然氣井數(shù)據(jù)驅(qū)動軟測量方法

王 丹 康 琦 楊居衡 宮 敬 張 奇

1. 中國石油大學(xué)（北京）經(jīng)濟管理學(xué)院 2. 中國石油經(jīng)濟技術(shù)研究院 3. 深圳清華大學(xué)研究院

4. 中國石油國際事業(yè)有限公司 5.油氣管道輸送安全國家工程實驗室·中國石油大學(xué)（北京）

0 引言

隨著海洋天然氣資源的開發(fā)，海上天然氣生產(chǎn)系統(tǒng)的運行監(jiān)測面臨著嚴峻挑戰(zhàn)：多相流量計安裝和維護成本高，用于系統(tǒng)運行分析的單井流量往往無法直接測量，只能得到與之有函數(shù)關(guān)系的溫度或壓力等變量的觀測數(shù)據(jù)，而觀測數(shù)據(jù)也受到各種誤差的干擾。以上因素導(dǎo)致生產(chǎn)系統(tǒng)運行數(shù)據(jù)獲取難度大、成本高、可靠性低，制約著海上氣田的安全生產(chǎn)[1-3]。為了解決工程中的此類問題，研究者們利用辨識建模法開發(fā)了多種計算速度快、在線適用性強的數(shù)據(jù)驅(qū)動模型，包括多項式模型和機器學(xué)習(xí)模型等[4-12]，表征系統(tǒng)狀態(tài)變量間的函數(shù)關(guān)系，用于工業(yè)系統(tǒng)關(guān)鍵狀態(tài)變量的軟測量，作為實體儀表的參考或備用。由于數(shù)據(jù)驅(qū)動模型的估計準確度有限，部分研究者基于最優(yōu)估計數(shù)據(jù)融合理論[13-14]，將濾波器引入軟測量，利用觀測數(shù)據(jù)實時校正模型模擬結(jié)果，充分利用工業(yè)系統(tǒng)多方面信息，得到更加貼近實際工況的估計結(jié)果。

濾波器是普遍應(yīng)用的、以最優(yōu)估計為基礎(chǔ)的動態(tài)數(shù)據(jù)融合方法，計算速度快，適用于動態(tài)系統(tǒng)的實時估計。濾波器自1960年產(chǎn)生以來，從僅適用于線性系統(tǒng)的卡爾曼濾波器[15]發(fā)展出了擴展卡爾曼濾波器[16]、集合卡爾曼濾波器[17]、無跡卡爾曼濾波器（Unscented Kalman Filter, UKF）[18-20]等一系列適用于非線性系統(tǒng)的濾波器[13,21]。多年來，濾波應(yīng)用范圍也從航天導(dǎo)航系統(tǒng)發(fā)展到電子、油氣生產(chǎn)及運輸?shù)裙I(yè)系統(tǒng)。Aguirre[22]等將濾波器與電子振蕩器的多項式數(shù)據(jù)驅(qū)動模型結(jié)合，用于改進電子振蕩器的模擬信號。才建[23]、王永紅[24]、Lorentzen[25-27]等將濾波器應(yīng)用于輸油和氣—液兩相流管道流量和壓力的估計，以及采油過程中井筒壓力和流量的估計。Teixeira[4]等將濾波器與多項式模型、多層感知器模型組成的數(shù)據(jù)驅(qū)動模型庫相結(jié)合，提升井底壓力軟測量效果；并利用模型冗余建立融合濾波器，進一步改善井底壓力估計結(jié)果。在油氣生產(chǎn)的數(shù)據(jù)驅(qū)動軟測量領(lǐng)域，目前研究者大多僅依賴數(shù)據(jù)驅(qū)動模型進行軟測量，對濾波器與數(shù)據(jù)驅(qū)動模型的結(jié)合，以及融合濾波器的應(yīng)用仍然較少，可見，軟測量準確度研究待進一步提高。

針對天然氣生產(chǎn)系統(tǒng)，引入UKF創(chuàng)新建立了數(shù)據(jù)驅(qū)動的軟測量方案，用于改善單井流量和井口壓力的估計準確度。①建立外源輸入非線性自回歸（Nonlinear Auto-Regressive with Exogenous Inputs，NARX）模型庫；②利用每個NARX模型構(gòu)建單個UKF；③利用NARX模型的冗余，創(chuàng)新建立兩種融合UKF。通過實際生產(chǎn)數(shù)據(jù)，驗證了濾波技術(shù)與數(shù)據(jù)驅(qū)動模型結(jié)合的優(yōu)勢，比較了單個UKF和融合UKF的準確性與計算成本，推薦適合工程應(yīng)用的最佳軟測量方案，從而提高流動管理系統(tǒng)的軟測量功能模塊，提高其獲取可靠運行數(shù)據(jù)的能力，以期推動智慧油氣田建設(shè)。

1 基于融合濾波的數(shù)據(jù)驅(qū)動估計算法建立

1.1 系統(tǒng)介紹和算法框架

以某海上天然氣兩井生產(chǎn)系統(tǒng)（圖1，兩井編號分別為Ⅰ、Ⅱ）為例建立基于融合濾波的數(shù)據(jù)驅(qū)動軟測量算法。該兩井生產(chǎn)系統(tǒng)位于某水下氣田，于2013年12月底投產(chǎn)。兩支路生產(chǎn)的天然氣—凝析液經(jīng)管匯進入海管，連入配有分離器的處理平臺。系統(tǒng)沿線分布多個配有壓力和溫度傳感器的觀測點，包括井底、井口、油嘴下游和平臺分離器，用數(shù)字1、2、3、4分別代表以上位置。單井未設(shè)置多相流量計，因此，單井流量不可觀測。

圖1 某海上天然氣兩井生產(chǎn)系統(tǒng)拓撲結(jié)構(gòu)圖

結(jié)合系統(tǒng)數(shù)據(jù)驅(qū)動模型和儀表測量值兩方面信息，兼顧模型和測量值的不確定度，基于數(shù)據(jù)融合理論，引入濾波技術(shù)建立“三步”數(shù)據(jù)驅(qū)動動態(tài)軟測量算法，用于單井流量和井口壓力的在線校正估算。所建軟測量算法是過程辨識和狀態(tài)估計技術(shù)的結(jié)合（圖2），其中，重構(gòu)井筒狀態(tài)空間方程和構(gòu)造合適的濾波器是兩個關(guān)鍵點。對于前者，對各井筒建立NARX模型庫用以重構(gòu)井筒狀態(tài)空間方程。對于后者，引入基于點集傳播的非線性濾波器以適應(yīng)非線性天然氣生產(chǎn)系統(tǒng)。

圖2 基于濾波器的數(shù)據(jù)驅(qū)動動態(tài)估計算法框架圖

1.2 NARX模型庫建立

建立NARX模型用以捕捉天然氣—凝析液在井筒中流動的非線性動態(tài)特性，描述當期輸出與先行輸出、先行輸入之間的動態(tài)約束關(guān)系[28-29]，為構(gòu)建狀態(tài)空間方程和濾波器奠定基礎(chǔ)。離散時間下NARX模型的通式表示為：

式中g(shù)表示NARX模型；R表示實數(shù)空間；RN×ζ表示N行ζ列的實數(shù)矩陣，其他同理；z(k)、u(k)、e(k)分別表示某時步下的輸出向量、輸入向量和誤差向量；z、u、e分別表示各時步下輸入向量、輸出向量、誤差向量組成的矩陣；ζ表示輸出變量的個數(shù)；μ表示輸入變量的個數(shù)；ρ=[ρ1, …,ρζ]表示輸出時滯向量；τ=[τ1, …,τμ]表示輸入時滯向量，每個輸出（輸入）變量的最大時滯不同Θ表示參數(shù)向量；k表示離散時步（k=1, 2, …,N，其中N表示采樣點個數(shù)）。

針對每個井筒，以單井流量和井口壓力為輸出變量（ζ=2），以系統(tǒng)沿線的可觀測變量為輸入變量，基于時間序列分析[30-31]，綜合動態(tài)和穩(wěn)態(tài)訓(xùn)練樣本（圖1），采用正交最小二乘回歸[28-29]和深度學(xué)習(xí)[9,10]等黑箱辨識技術(shù)，以及雙目標最小二乘[5,32]和遷移學(xué)習(xí)等參數(shù)校正技術(shù)，建立了一系列相互獨立的黑箱和灰箱NARX模型（gj，j=1, 2, …,m，其中m為模型個數(shù)），組成NARX模型庫。模型庫中包括：以多項式（Polynomial，POLY）為基礎(chǔ)的POLY-NARX黑箱動態(tài)模型（g1）、POLY-NARX灰箱動態(tài)模型（g2），以深度前饋網(wǎng)絡(luò)（Dense Neural Network，DNN）為基礎(chǔ)的DNN-NARX黑箱動態(tài)模型（g3、g4），DNNNARX灰箱動態(tài)模型（g5、g6）。其中，g3和g5以支路和平臺可觀測變量為輸入變量，g4和g6僅以支路可觀測變量為輸入變量。以上NARX模型相關(guān)建模過程和模型參數(shù)設(shè)置采用Wang等[33]研究成果。

1.3 UKF理論基礎(chǔ)

對于非線性的天然氣生產(chǎn)系統(tǒng)，經(jīng)典卡爾曼濾波器（Kalman Filter，KF）很多條件和性質(zhì)不再成立[13,21]。1990年后出現(xiàn)的Sigma點卡爾曼濾波器（Sigma Point Kalman Filter，SPKF）[21]適用于非線性系統(tǒng)。SPKF通過確定性采樣法得到Sigma點集，用Sigma點集非線性映射后的期望和方差近似狀態(tài)向量（X）概率分布。SPKF具有以下優(yōu)勢：①逼近概率分布比逼近任意非線性函數(shù)更容易實現(xiàn)，避免對非線性函數(shù)進行解析求導(dǎo)；②對均值和方差的逼近精度提高到二階及以上；③計算量較小，濾波穩(wěn)定性較高[13]。

Julier和Uhlmann于1994年提出的UKF是典型的SPKF方法，利用歷史狀態(tài)向量的后驗分布，采用無跡變換法（Unscented Transform，UT）生成規(guī)模為2nX+1（nX為X維度）、對稱分布的Sigma點集[20-21]，并以此進行非線性傳播。

1.3.1 狀態(tài)空間方程和濾波步驟

系統(tǒng)的非線性狀態(tài)空間方程是構(gòu)建UKF的基礎(chǔ)，由非線性狀態(tài)轉(zhuǎn)移方程和觀測方程組成[式（2）]，前者表示不同時步狀態(tài)向量間的遞推關(guān)系，后者表示狀態(tài)向量和觀測向量間的關(guān)系。

式中X(k)、U(k－1)、Y(k)表示某時步下的狀態(tài)向量、輸入向量、觀測向量；nX、nU、nY分別表示狀態(tài)變量、輸入變量和觀測變量的個數(shù)；C(k)表示某時步下的觀測矩陣；W(k－1)、V(k)分別表示某時步下的系統(tǒng)模型噪聲向量、觀測噪聲向量，假設(shè)是均值為0的白噪。

濾波是一種基于最小化均方誤差[13]的“兩步”遞推估計算法。第一步是一步預(yù)測，也叫時間更新，根據(jù)狀態(tài)向量前一時步的分布信息Px(k－1) }得到其一步預(yù)測值，即當前時步的先驗估計。第二步是測量更新，以濾波增益為權(quán)重，融合觀測向量的測量新息與狀態(tài)向量的先驗估計，得到狀態(tài)向量當前時步的后驗估計。確定狀態(tài)向量的初始分布Px(0)}后，按以上兩步不斷進行遞推估算，可得到每個時步狀態(tài)向量的后驗估計。用式（3）概括UKF的濾波步驟，具體計算過程詳見文獻[21]。

式中Forecast表示一步預(yù)測階段；DataAssimilation表示數(shù)據(jù)融合，也就是測量校正階段；Q(k－1)、R(k)分別表示某時步下系統(tǒng)模型噪聲和觀測噪聲的協(xié)方差陣，是系統(tǒng)模型和觀測數(shù)據(jù)的不確定度；(k)表示某時步下觀測向量的測量值；分別表示Sigma點集、觀測向量點集的期望；Px(k)表示狀態(tài)向量的誤差協(xié)方差陣；Py(k)表示觀測向量的誤差協(xié)方差陣；Pxy(k)表示狀態(tài)向量和觀測向量的誤差協(xié)方差陣；S(k)表示濾波增益，S(k)=Pxy(k)P-y1(k)；某時步下Y(k)的新息為，即觀測向量的估計殘差；k|k－1代表k－1時步到k時步的一步預(yù)測值或先驗估計。

1.3.2 噪聲統(tǒng)計特性設(shè)置

濾波算法是利用狀態(tài)空間方程和白噪聲激勵的統(tǒng)計特性形成。在白噪激勵的假設(shè)下，Q(k)=diag[qi2(k)]qi2(k)表示系統(tǒng)模型對每個狀態(tài)變量一步預(yù)測的誤差方差，通過模型模擬結(jié)果得到。ro2(k)表示每個觀測變量的噪聲方差，通過儀表的不確定度得到。

1.4 基于融合UKF的數(shù)據(jù)驅(qū)動動態(tài)估計

基于所建NARX模型庫重構(gòu)井筒狀態(tài)空間方程，以此構(gòu)建單個UKF。進而利用NARX模型的冗余，建立融合UKF。

1.4.1 天然氣生產(chǎn)井狀態(tài)空間方程重構(gòu)

以天然氣生產(chǎn)井為估計單元，對每個井筒構(gòu)建基于NARX模型的狀態(tài)空間方程。

1.4.1.1 狀態(tài)轉(zhuǎn)移方程

狀態(tài)轉(zhuǎn)移方程的一步預(yù)測算子（f）通過時滯的NARX模型（gj）變形得到[4,21]，gj的輸入向量、輸出向量與f的輸入向量、狀態(tài)向量之間具有明確的對應(yīng)關(guān)系[式（4）]。k時步下，U(k)對應(yīng)u(k)張成的空間，X(k)對應(yīng)z(k)張成的空間。z=[G, p2]，則X=[G, p2]；不同狀態(tài)轉(zhuǎn)移方程對應(yīng)的gj不同，相應(yīng)的輸入向量也不同。

式中wi(i=1, 2, …, nX)表示系統(tǒng)模型對每個狀態(tài)變量一步預(yù)測的誤差。

1.4.1.2 觀測方程

為便于驅(qū)動濾波器，取狀態(tài)向量中的可觀測變量組成觀測向量[22,25]，因此Y=[z2]=[p2]。Y(k)=[z2(k), …, z2(k－ ρ2+1)]= [p2(k), …, p2(k－ ρ2+1)]，C(k)=[0，I]，其中0是全零列，I是單位陣。

1.4.1.3 初始分布與噪聲協(xié)方差陣

取各狀態(tài)變量樣本集中第1個時刻的值作為狀態(tài)向量的初始期望ˉ(0)。狀態(tài)向量的初始誤差協(xié)方差陣 Px(0)=diag[σi2(0) ]，σi2(0)=10－5。

Qk中，qi2(k)根據(jù)各模型對輸出變量模擬誤差的方差賦值[4]。R(k)中ro2(k)=0.003 6。

1.4.2 融合UKF建立

為發(fā)揮多模型優(yōu)勢，提高軟測量的估計性能，利用所建NARX模型庫中模型的冗余，構(gòu)造了兩種融合UKF。第一種融合濾波器“先融合后濾波”，簡稱為FF1（圖3）；第二種融合濾波器“先濾波后融合”[4]，簡稱為FF2（圖3）。

圖3 融合濾波器FF1和FF2的計算流程圖

1.4.2.1 建立FF1

針對每個井筒構(gòu)建1個UKF，該UKF的狀態(tài)轉(zhuǎn)移方程是模型庫中多個模型的組合。模型庫中擇優(yōu)選擇多個模型并行計算，對所得多個模型模擬值進行加權(quán)融合得到一步預(yù)測值，在此基礎(chǔ)上進行測量校正得到后驗估計值。各時步下每個模型的相應(yīng)權(quán)重根據(jù)各自的模擬誤差賦值，誤差較小的模型對應(yīng)較大的權(quán)重。

1.4.2.2 建立FF2

每個模型分別構(gòu)建1個UKF，通過多個模型得到多個狀態(tài)空間方程互不相同的獨立UKF，組成濾波庫再進行融合。擇優(yōu)選擇多個獨立UKF并行計算，對所得多個UKF的后驗估計結(jié)果進行加權(quán)融合，得到融合濾波值。各時步下每個濾波器的相應(yīng)權(quán)重根據(jù)各自的新息賦值，新息較小的濾波器對應(yīng)較大的權(quán)重。

2 算例分析和結(jié)果

以上述海上天然氣兩井生產(chǎn)系統(tǒng)（圖1）中Ⅱ井的生產(chǎn)數(shù)據(jù)為例，從可行性、準確度和計算速度等方面對所建數(shù)據(jù)驅(qū)動估計算法進行驗證。根據(jù)前期研究結(jié)果，灰箱模型的全局準確性和穩(wěn)定性高于相應(yīng)黑箱模型[33]，因此，直接基于灰箱模型g2、g5、g6構(gòu)建3個相互獨立的UKF（表1），分別命名為UKF2、UKF5、UKF6。對以上灰箱模型進一步擇優(yōu)，用于建立融合濾波器FF1和FF2。在此基礎(chǔ)上比較單個UKF和FF1、FF2的估計結(jié)果與計算時間。

表1 Ⅱ井灰箱模型g2、g5、g6和濾波器UKF2、UKF5、UKF6設(shè)置表

利用生產(chǎn)系統(tǒng)的常規(guī)工藝參數(shù)和生產(chǎn)數(shù)據(jù)進行算法測試。常規(guī)工藝參數(shù)包括井身結(jié)構(gòu)和流體組分，生產(chǎn)數(shù)據(jù)由集散控制系統(tǒng)采集的儀表觀測數(shù)據(jù)和兩相流模擬軟件OLGA的模擬數(shù)據(jù)組成，共涵蓋了2014年1—4月的5 000組動態(tài)數(shù)據(jù)，并從中提取出10組穩(wěn)態(tài)生產(chǎn)數(shù)據(jù)（圖4）。其中，儀表觀測數(shù)據(jù)包括T1、p1、T2、p2、T3、p3、T4、p4、V；OLGA 的 模擬數(shù)據(jù)主要包括G。以2014年1—4月p2的儀表測量值和G的OLGA模擬值作為算法估算結(jié)果的對照，為評判估算結(jié)果的準確性提供依據(jù)。受限于現(xiàn)場生產(chǎn)條件，各生產(chǎn)井不配備流量計，因此目前難以通過現(xiàn)場生產(chǎn)數(shù)據(jù)來驗證OLGA計算結(jié)果。OLGA是工業(yè)界公認的多相流模擬專業(yè)軟件，現(xiàn)有條件下相對可行的手段是利用OLGA的流量模擬結(jié)果來驗證算法的流量估算結(jié)果。

圖4 Ⅱ井2014年1—4月動態(tài)和穩(wěn)態(tài)生產(chǎn)數(shù)據(jù)圖

使用均方根誤差（Root Mean Square Error，RMSE）、平均絕對誤差（Mean Absolute Error，MAE）、均方誤差（Mean Square Error，MSE）及決定系數(shù)R2作為算法準確度的評價指標[34]。

2.1 不同NARX模型所建UKF對比

將UKF2、UKF5、UKF6對Ⅱ井流量（G）和井口壓力（p2）的估算結(jié)果（圖5）、估算準確度指標（表2）與相應(yīng)NARX模型g2、g5、g6的結(jié)果一一對比。與g2、g5、g6的模擬結(jié)果相比，具有測量校正功能的UKF2、UKF5、UKF6對G和p2的估算準確度明顯提高，R2較大且RMSE、MRE、MAE較小。其中，UKF5和UKF6的R2達到0.9以上，而g5、g6的R2約為0.7、0.5；與g6相比，UKF6的RMSE、MRE、MAE減小50%以上。以上結(jié)果說明，通過濾波器反饋井口壓力測量值能夠有效校正模型對單井流量和井口壓力的模擬結(jié)果。

表2 Ⅱ井模型g2、g5、g6和濾波器UKF2、UKF5、UKF6的計算準確度指標統(tǒng)計表

圖5 Ⅱ井模型g2、g5、g6和濾波器UKF2、UKF5、UKF6對單井流量和井口壓力的計算結(jié)果圖

對不同模型構(gòu)建的UKF2、UKF5、UKF6進行對比可見，通過DNN-NARX灰箱模型構(gòu)建的UKF5、UKF6的估算準確度高于POLY-NARX灰箱模型構(gòu)建的UKF2。進一步比較UKF5與UKF6發(fā)現(xiàn)，UKF5的輸入向量中包含支路變量和平臺變量，而UKF6的輸入向量僅包含支路變量。雖然UKF6的輸入變量個數(shù)較少，但與UKF5的估算誤差相差不大。因此，DNN-NARX模型構(gòu)造的UKF具有較高的準確性和一定的穩(wěn)定性，在工程實際中，如果部分觀測點的儀表缺失導(dǎo)致輸入變量個數(shù)減少，仍能保持較高的估算準確度。因此，選擇DNN-NARX灰箱模型作為融合濾波器的首選模型。

2.2 單個UKF與融合UKF對比

基于DNN-NARX灰箱模型g5、g6構(gòu)建融合濾波器FF1、FF2，比較單個濾波器UKF5、UKF6和融合濾波器FF1、FF2對G和p2的估算準確度（圖6、表3）與計算速度。

表3 單個濾波器UKF5、UKF6和融合濾波器FF1、FF2的計算準確度指標統(tǒng)計表

圖6 Ⅱ井單個濾波器UKF5、UKF6與融合濾波器FF1、FF2對單井流量和井口壓力的計算結(jié)果圖

根據(jù)結(jié)果可知，兩個融合濾波器中，F(xiàn)F1的R2較高且RMSE、MAE、MRE等估算誤差明顯較小。FF1與單個濾波器中估算效果較好的UKF5相比，前者R2較高，且MAE和MRE較小。而FF2未體現(xiàn)絕對優(yōu)勢，其估算準確度介于UKF5和UKF6之間。因此，融合濾波器FF1的準確度最高。

對單個UKF和融合UKF的計算時間進行對比。計算程序使用Java語言進行編寫，在Windows 7 64-bit操作系統(tǒng)平臺上搭建JDK開發(fā)環(huán)境，測試環(huán)境為4核PC，處理器性能為Intel（R）Core（TM）CPU i5-5200 @ 2.20 GHz。UKF5、UKF6、FF1、FF2 的單步計算時間分別為8 ms、8 ms、10 ms、15 ms。可以看出，F(xiàn)F1和FF2比單個UKF計算耗時，但均能滿足工程中在線估計的需求。同時發(fā)現(xiàn)，F(xiàn)F2比FF1耗時長，這是因為，F(xiàn)F1是單個濾波器運行，而FF2是多濾波器同時并行，且每個濾波器涉及Sigma點集的并行計算，計算量較大。可以推斷，隨著FF1和FF2所含模型個數(shù)增加，與單個UKF計算時間的差距會加大。

綜上所述，工程應(yīng)用中首先推薦FF1作為數(shù)據(jù)驅(qū)動軟測量方案。FF1能夠綜合不同模型的優(yōu)勢，整體估算準確度最高，計算量和耗時較少。其次推薦DNN-NARX灰箱模型構(gòu)造的單個UKF作為數(shù)據(jù)驅(qū)動軟測量方案，該方案也具有較高的估算準確度和較快的計算速度。

3 結(jié)論

1）基于UKF軟測量方案的準確度高于純模型軟測量方案的準確度，所建軟測量算法準確度較高，計算量小，拓展性和在線適用性強，能夠提高流動管理系統(tǒng)獲取可靠運行數(shù)據(jù)的能力。

2）所建包含POLY-NARX和DNN-NARX的模型庫中，DNN-NARX灰箱模型構(gòu)造的UKF準確性和穩(wěn)定性較高。

3）單個UKF與FF1、FF2相比，“先融合后濾波”的FF1估計準確度較高，計算速度與單個UKF相當，推薦作為首選軟測量方案。