一種基于自適應擴展卡爾曼濾波的ect圖像重建方法
【專利摘要】本發(fā)明公開了一種基于自適應擴展卡爾曼濾波的ECT圖像重建方法，包括以下步驟：建立系統(tǒng)方程和測量方程，擴展卡爾曼迭代濾波，直至達到預設的迭代終止條件。本發(fā)明采用抗噪性能較好的擴展卡爾曼濾波方法對電容層析成像進行圖像重建，在建立系統(tǒng)方程時加入系統(tǒng)噪聲，并在線修正系統(tǒng)噪聲的協(xié)方差矩陣，使狀態(tài)空間模型更能反映真實的多相流運動，提高圖像重建精度。
【專利說明】
一種基于自適應擴展卡爾曼濾波的ECT圖像重建方法
技術(shù)領(lǐng)域
[0001] 本發(fā)明涉及一種ECT圖像重建方法，尤其是一種基于自適應擴展卡爾曼濾波的ECT 圖像重建方法，屬于兩相流/多相流運動的可視化檢測技術(shù)領(lǐng)域。
【背景技術(shù)】
[0002] 電容層析成像(Electrical Capacitance Tomography，ECT)是20世紀80年代發(fā)展 起來的一種檢測技術(shù)，該技術(shù)具有非侵入性、無輻射、響應快、成本低等優(yōu)點，多用于工業(yè)兩 相流/多相流運動檢測中。工業(yè)中常見的兩相流/多相流運動包括:氣/液兩相流、氣/固兩相 流、液/固兩相流、液/液兩相流、氣/液/固三相流和油/氣/水三相流。兩相流/多相流檢測中 需要檢測的主要參數(shù)包括:流型、分相含率、流量、流速、密度、溫度和壓強等，其中精確的流 型識別對分相含率等其他參數(shù)的檢測具有重要意義。兩相流/多相流的流型通過觀察流體 流動的幾何形態(tài)進行識別，目前主要的識別方法有三種:通過視覺觀測識別流型、通過測量 流體中某些參數(shù)的變化識別流型、通過層析成像技術(shù)識別流型，電容層析成像即是一種層 析成像技術(shù)。電容層析成像技術(shù)進行多相流檢測時主要分為兩個步驟:正問題和逆問題。其 中，正問題通過介電常數(shù)分布確定電極間電容值;逆問題通過測量電容值確定介電常數(shù)分 布，用以識別兩相流/多相流的流型。
[0003] 由于逆問題求解過程中的欠定性、病態(tài)性和"軟場"效應，許多圖像重建算法被提 出，如Landweber迭代算法、Tikhonov正則化算法等。這些算法雖在仿真實驗中取得質(zhì)量不 錯的重建圖像，但在現(xiàn)場應用中卻因為抗噪性不佳，重建效果大打折扣。Kalman濾波方法利 用噪聲的統(tǒng)計特性對圖像灰度值進行估計，具有較好的抗噪性，但Kalman濾波能夠得到最 優(yōu)估計值的前提條件是建立準確的狀態(tài)空間模型。目前建立的ECT狀態(tài)空間模型大體分為 兩種情況：
[0004] 1、未考慮系統(tǒng)噪聲。
[0005] 2、雖然考慮了系統(tǒng)噪聲，但將系統(tǒng)噪聲的協(xié)方差矩陣設置為常數(shù)。
[0006] 但當多相流運動時，流型會隨時間發(fā)生變化，這兩種狀態(tài)空間模型都無法準確的 描述多相流運動情況。因此，要解決這一問題，需要在系統(tǒng)模型中添加系統(tǒng)噪聲，并實時地 在線修正系統(tǒng)噪聲的協(xié)方差矩陣，以使狀態(tài)空間模型與真實的多相流運動更加接近。

【發(fā)明內(nèi)容】

[0007] 本發(fā)明要解決的技術(shù)問題是提供一種基于自適應擴展卡爾曼濾波的ECT圖像重建 方法。
[0008] 為解決上述技術(shù)問題，本發(fā)明采用的技術(shù)方案是：
[0009] 一種基于自適應擴展卡爾曼濾波的ECT圖像重建方法，包括以下步驟：
[0010] 步驟1:建立系統(tǒng)方程和測量方程，由以下具體步驟組成：
[0011] 步驟1-1:建立系統(tǒng)方程：
[0012] gk = gk-i+?k-i (1)
[0013] 式中，gk和分別代表k時刻和k_l時刻的圖像灰度值；為k_l時刻的系統(tǒng)噪 聲，所述系統(tǒng)噪聲為白噪聲，其均值為〇，協(xié)方差矩陣為Qk;
[0014] 步驟1-2:建立非線性方程：
[0015] Vk = Uk(gk)+vk (2)
[0016] 式中，Vk為k時刻的電容測量值;Uk(gk)表示圖像灰度電容測量值與之間的非線性 關(guān)系;vk為k時刻的測量噪聲，所述測量噪聲為白噪聲，其均值為0,協(xié)方差矩陣為R;
[0017] 步驟1-3:對所述步驟1-2建立的非線性方程進行泰勒級數(shù)展開：
[0018] Vk = Uk(go)+Jk(go) · (gk-go)+H.O.Ts+vk (3)
[0019]式中，H.O.Ts表示高階項，為高斯白噪聲;Jk(g〇)為Jacobian矩陣，定義為：
(4)
[0021] 步驟1-4:建立測量方程：
[0022] z:k=Sgk+vk (5)
[0023] 式中，Zk為歸一化測量值；S為歸一化靈敏度矩陣；gk為歸一化圖像灰度值； 5 + ，為線性測量噪聲，其均值為〇,方差為歹；：
[0024] 步驟2:擴展卡爾曼迭代濾波，由以下具體步驟組成：
[0025]步驟2-1:設置濾波初值：圖像灰度初值go為0向量，誤差協(xié)方差矩陣初值Po為α · I， I為單位矩陣，α為比例系數(shù)。
[0026]步驟2-2:計算系統(tǒng)噪聲的協(xié)方差矩陣Qk的更新方程：
[0027] Qk^Kk-t-K^(β)
[0028] 式中，么表示新息序列％ =?-為的方差矩陣，zk表示k時刻的歸一化測量值； 表示從時刻k-Ι到時刻k的一步量測預測值;K k為濾波增益矩陣;Pk和Ph分別表示時刻k 和時刻k-1的濾波誤差協(xié)方差矩陣；
[0029]步驟2-3:擴展卡爾曼濾波：
<7)：
[0031 ]式中，么和分別表示時刻k和時刻k-Ι的圖像灰度的估計值;么/μ表示從時刻k 到時刻k-1圖像灰度的一步預測值;Pk/k-i表示從k-1時刻到k時刻的一步預測協(xié)方差;〇!^表 不k_l時刻系統(tǒng)噪聲協(xié)方差矩陣；
[0032]步驟3:判斷是否達到預設的迭代終止條件為終止閾值，如果是，轉(zhuǎn) 向步驟4;否則轉(zhuǎn)向步驟2;
[0033]步驟4:輸出圖像灰度的最優(yōu)估計值。
[0034] 所述步驟2-2中序列新息序的方差矩陣為：
[0035] Pvk=SPklkAST +R (8)
[0036] 所述步驟2-2中序列新息序列％ =? -#i/M的方差矩陣的計算方法為：
[0037] PvM = E[{zk - zk/k_^ ){zk - zk,k_A) ] (9)
[0038] 式中，電容測量值Zk表示為：
[0039] zk = HkXk+vk (10)
[0040] 一步測量預測值i表示為：
[0041] = ^k^k'k-i (.1.1.)
[0042] 式中，Hk表示k時刻的狀態(tài)轉(zhuǎn)移矩陣Λ%表示k-1時刻到k時刻的狀態(tài)預測值。
[0043] 采用上述技術(shù)方案所產(chǎn)生的有益效果在于：
[0044] 本發(fā)明采用抗噪性能較好的擴展卡爾曼濾波方法對電容層析成像進行圖像重建， 在建立系統(tǒng)方程時加入系統(tǒng)噪聲，并在線修正系統(tǒng)噪聲的協(xié)方差矩陣，使狀態(tài)空間模型更 能反映真實的多相流運動，提高圖像重建精度。
【附圖說明】
[0045] 圖1是本發(fā)明的流程圖。
【具體實施方式】
[0046] 下面結(jié)合附圖和【具體實施方式】對本發(fā)明作進一步詳細的說明。
[0047] 一種基于自適應擴展卡爾曼濾波的ECT圖像重建方法，包括以下步驟：
[0048] 步驟1:建立系統(tǒng)方程和測量方程，由以下具體步驟組成：
[0049] 步驟1-1:建立系統(tǒng)方程：
[0050] gk = gk-i+?k-i (1)
[0051] 式中，gk和gk-^別代表k時刻和k_l時刻的圖像灰度值；為k_l時刻的系統(tǒng)噪 聲，所述系統(tǒng)噪聲為白噪聲，其均值為〇，協(xié)方差矩陣為Qk;
[0052] 所述系統(tǒng)方程利用系統(tǒng)噪聲反映多相流流型變化，在建立系統(tǒng)方程時，加入系統(tǒng) 噪聲，以反映當流體流型發(fā)生變化時，系統(tǒng)模型的改變，其噪聲方差陣未知，需要進行估計； [0053] 步驟1-2:建立非線性方程：
[0054] Vk = Uk(gk)+vk (2)
[0055] 式中，Vk為k時刻的電容測量值;Uk(gk)表示圖像灰度電容測量值與之間的非線性 關(guān)系;vk為k時刻的測量噪聲，所述測量噪聲為白噪聲，其均值為0,協(xié)方差矩陣為R;
[0056] 所述非線性方程用于描述電容測量值與圖像灰度關(guān)系；電容測量值與圖像灰度值 之間為非線性關(guān)系，利用測量噪聲表示測量過程中產(chǎn)生的噪聲，認為測量噪聲的方差矩陣R 在測量過程中不變；
[0057]步驟1-3:對所述步驟1-2建立的非線性方程進行泰勒級數(shù)展開：
[0058] Vk = Uk(go)+Jk(go) · (gk-go)+H.O.Ts+vk (3)
[0059] 式中，H.O.Ts表示高階項，為高斯白噪聲;Jk(go)為Jacobian矩陣，定義為：
(4)
[0061 ] 步驟1-4:建立測量方程：
[0062] zk = Sgk + vk (, 5 )
[0063] 式中，Zk為歸一化測量值；S為歸一化靈敏度矩陣；gk為歸一化圖像灰度值； R ^//λΤΙν + ν)，為線性測量噪聲，其均值為〇,方差為瓦；
[0064] 步驟2:擴展卡爾曼迭代濾波，由以下具體步驟組成：
[0065] 步驟2-1:設置濾波初值:圖像灰度初值go為0向量，誤差協(xié)方差矩陣初值Ρο為α · I， I為單位矩陣，α為比例系數(shù)。
[0066]步驟2-2:計算系統(tǒng)噪聲的協(xié)方差矩陣Qk的更新方程：
[0067] (6)
[0068] 式中，?表示新息序列％ =? 的方差矩陣，zk表示k時刻的歸一化測量值； 表示從時刻k_l到時刻k的一步量測預測值;Kk為濾波增益矩陣;Pk和Ph分別表示時刻k 和時刻k-1的濾波誤差協(xié)方差矩陣；
[0069]步驟2-3:擴展卡爾曼濾波： Sk/t-ι - Sk-i Λ A· my r- Λ'* A. -| Sk = OA-/i--l + Kkizt
[0070] '、人廠 L fpu7 十歹]! ( 7 ) Pk.tk-l ~ Pk-1 + Qt-\ ，=[/-A'd 丨
[0071] 式中，分別表示時刻k和時刻k-1的圖像灰度的估計值;表示從時刻k 到時刻k-1圖像灰度的一步預測值;Pk/k-i表示從k-1時刻到k時刻的一步預測協(xié)方差;〇!^表 不k_l時刻系統(tǒng)噪聲協(xié)方差矩陣；
[0072]本發(fā)明基于極大似然準則，對系統(tǒng)噪聲的協(xié)方差矩陣Qk進行實時修正；
[0073] 步驟3 :判斷是否達到預設的迭代終止條件Ih S，在本實施例中ε =〇 . 001， 如果是，轉(zhuǎn)向步驟4;否則轉(zhuǎn)向步驟2;
[0074] 步驟4:輸出圖像灰度的最優(yōu)估計值。
[0075] 所述步驟2-2中序列新息序列％ =? 的方差矩陣為：
[0076] Pvk=SPk/k^ST+R (8)
[0077] 所述步驟2-2中序列新息序列％ ，的方差矩陣的計算方法為：
[0078] Pvk=E\izk-zki kA){zk~zk,kl) ] (9)
[0079] 式中，電容測量值Zk表示為：
[0080] zk = HkXk+vk (10)
[0081] 一步測量預測值心Η表示為：
[0082] %^=Hkxktk_l (115
[0083] 式中，Hk表示k時刻的狀態(tài)轉(zhuǎn)移矩陣，毛/w表示k-1時刻到k時刻的狀態(tài)預測值。
[0084] 所述擴展卡爾曼濾波方程基于誤差方差最小的準則得到。
【主權(quán)項】
1. 一種基于自適應擴展卡爾曼濾波的ECT圖像重建方法，其特征在于： 包括W下步驟： 步驟1:建立系統(tǒng)方程和測量方程，由W下具體步驟組成： 步驟1-1:建立系統(tǒng)方程： 阱= gk-l+?k-l (1) 式中，gk和gk-1分別代表k時刻和k-1時刻的圖像灰度值；ω k-1為k-1時刻的系統(tǒng)噪聲，所 述系統(tǒng)噪聲為白噪聲，其均值為0,協(xié)方差矩陣為Qk; 步驟1-2:建立非線性方程： Vk = Uk(阱）+vk (2) 式中，Vk為k時刻的電容測量值;化(gk)表示圖像灰度電容測量值與之間的非線性關(guān)系； vk為k時刻的測量噪聲，所述測量噪聲為白噪聲，其均值為0，協(xié)方差矩陣為R; 步驟1-3:對所述步驟1-2建立的非線性方程進行泰勒級數(shù)展開： Vk = Uk(go)+Jk(go) ·(阱-go)+H.O.Ts+vk (3) 式中，H.O.Ts表示高階項，為高斯白噪聲;Jk(go)為化cobian矩陣，定義為：式中，zk為歸一化測量值；S為歸一化靈敏度矩陣；gk為歸一化圖像灰度值； 町=//·ατ：ν + ν',,，為線性測量噪聲，其均值為0,方差為玄； 步驟2:擴展卡爾曼迭代濾波：由W下具體步驟組成： 步驟2-1:設置濾波初值：圖像灰度初值go為0向量，誤差協(xié)方差矩陣初值Ρο為α · 1，1為 單位矩陣，α為比例系數(shù)。 步驟2-2:計算系統(tǒng)噪聲的協(xié)方差矩陣Qk的更新方程：式中，每表示新息序列呢=毎-％…的方差矩陣，Zk表示k時刻的歸一化測量值A…表 示從時刻k-1到時刻k的一步量測預測值;Kk為濾波增益矩陣;Pk和Pk-i分別表示時刻k和時刻 k-1的濾波誤差協(xié)方差矩陣； 步驟2-3:擴展卡爾曼濾波：(7) 式中，么和為_1分別表示時刻k和時刻k-1的圖像灰度的估計值;么表示從時刻k到時 亥化-1圖像灰度的一步預測值;Pk/k-i表示從k-1時刻到k時刻的一步預測協(xié)方差;Qk-i表示k-1 時刻系統(tǒng)噪聲協(xié)方差矩陣； 步驟3:判斷是否達到預設的迭代終止條件I私-耗為終止闊值，如果是，轉(zhuǎn)向步 驟4;否則轉(zhuǎn)向步驟2; 步驟4:輸出圖像灰度的最優(yōu)估計值。 所述步驟2-2中序列新息序列成=? ~S.w_，的方差矩陣為：(8) 所述步驟2-2中序列新息序列％ =句-^;._1的方差矩陣的計算方法為：(9) 式中，電容測量值Zk表示為： Zk =化Hi+Vk (10) 一步測量預測值表示為： =巧A-/1--1 U1) 式中，化表示k時刻的狀態(tài)轉(zhuǎn)移矩陣，爲康示k-1時刻至化時刻的狀態(tài)預測值。
【文檔編號】G06T5/10GK106097285SQ201610575305
【公開日】2016年11月9日
【申請日】2016年7月21日
【發(fā)明人】張立峰
【申請人】華北電力大學（保定）