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

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

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