一种基于自适应扩展卡尔曼滤波的ect图像重建方法
【专利摘要】本发明公开了一种基于自适应扩展卡尔曼滤波的ECT图像重建方法，包括以下步骤：建立系统方程和测量方程，扩展卡尔曼迭代滤波，直至达到预设的迭代终止条件。本发明采用抗噪性能较好的扩展卡尔曼滤波方法对电容层析成像进行图像重建，在建立系统方程时加入系统噪声，并在线修正系统噪声的协方差矩阵，使状态空间模型更能反映真实的多相流运动，提高图像重建精度。
【专利说明】
一种基于自适应扩展卡尔曼滤波的ECT图像重建方法
技术领域
[0001] 本发明涉及一种ECT图像重建方法，尤其是一种基于自适应扩展卡尔曼滤波的ECT 图像重建方法，属于两相流/多相流运动的可视化检测技术领域。
【背景技术】
[0002] 电容层析成像(Electrical Capacitance Tomography，ECT)是20世纪80年代发展 起来的一种检测技术，该技术具有非侵入性、无辐射、响应快、成本低等优点，多用于工业两 相流/多相流运动检测中。工业中常见的两相流/多相流运动包括:气/液两相流、气/固两相 流、液/固两相流、液/液两相流、气/液/固三相流和油/气/水三相流。两相流/多相流检测中 需要检测的主要参数包括:流型、分相含率、流量、流速、密度、温度和压强等，其中精确的流 型识别对分相含率等其他参数的检测具有重要意义。两相流/多相流的流型通过观察流体 流动的几何形态进行识别，目前主要的识别方法有三种:通过视觉观测识别流型、通过测量 流体中某些参数的变化识别流型、通过层析成像技术识别流型，电容层析成像即是一种层 析成像技术。电容层析成像技术进行多相流检测时主要分为两个步骤:正问题和逆问题。其 中，正问题通过介电常数分布确定电极间电容值;逆问题通过测量电容值确定介电常数分 布，用以识别两相流/多相流的流型。
[0003] 由于逆问题求解过程中的欠定性、病态性和"软场"效应，许多图像重建算法被提 出，如Landweber迭代算法、Tikhonov正则化算法等。这些算法虽在仿真实验中取得质量不 错的重建图像，但在现场应用中却因为抗噪性不佳，重建效果大打折扣。Kalman滤波方法利 用噪声的统计特性对图像灰度值进行估计，具有较好的抗噪性，但Kalman滤波能够得到最 优估计值的前提条件是建立准确的状态空间模型。目前建立的ECT状态空间模型大体分为 两种情况：
[0004] 1、未考虑系统噪声。
[0005] 2、虽然考虑了系统噪声，但将系统噪声的协方差矩阵设置为常数。
[0006] 但当多相流运动时，流型会随时间发生变化，这两种状态空间模型都无法准确的 描述多相流运动情况。因此，要解决这一问题，需要在系统模型中添加系统噪声，并实时地 在线修正系统噪声的协方差矩阵，以使状态空间模型与真实的多相流运动更加接近。
【发明内容】
[0007] 本发明要解决的技术问题是提供一种基于自适应扩展卡尔曼滤波的ECT图像重建 方法。
[0008] 为解决上述技术问题，本发明采用的技术方案是：
[0009] 一种基于自适应扩展卡尔曼滤波的ECT图像重建方法，包括以下步骤：
[0010] 步骤1:建立系统方程和测量方程，由以下具体步骤组成：
[0011] 步骤1-1:建立系统方程：
[0012] gk = gk-i+?k-i (1)
[0013] 式中，gk和分别代表k时刻和k_l时刻的图像灰度值；为k_l时刻的系统噪 声，所述系统噪声为白噪声，其均值为〇，协方差矩阵为Qk;
[0014] 步骤1-2:建立非线性方程：
[0015] Vk = Uk(gk)+vk (2)
[0016] 式中，Vk为k时刻的电容测量值;Uk(gk)表示图像灰度电容测量值与之间的非线性 关系;vk为k时刻的测量噪声，所述测量噪声为白噪声，其均值为0,协方差矩阵为R;
[0017] 步骤1-3:对所述步骤1-2建立的非线性方程进行泰勒级数展开：
[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向量，误差协方差矩阵初值Po为α · I， I为单位矩阵，α为比例系数。
[0026]步骤2-2:计算系统噪声的协方差矩阵Qk的更新方程：
[0027] Qk^Kk-t-K^(β)
[0028] 式中，么表示新息序列％ =?-為的方差矩阵，zk表示k时刻的归一化测量值； 表示从时刻k-Ι到时刻k的一步量测预测值;K k为滤波增益矩阵;Pk和Ph分别表示时刻k 和时刻k-1的滤波误差协方差矩阵；
[0029]步骤2-3:扩展卡尔曼滤波：
<7)：
[0031 ]式中，么和分别表示时刻k和时刻k-Ι的图像灰度的估计值;么/μ表示从时刻k 到时刻k-1图像灰度的一步预测值;Pk/k-i表示从k-1时刻到k时刻的一步预测协方差;〇!^表 不k_l时刻系统噪声协方差矩阵；
[0032]步骤3:判断是否达到预设的迭代终止条件为终止阈值，如果是，转 向步骤4;否则转向步骤2;
[0033]步骤4:输出图像灰度的最优估计值。
[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时刻的状态转移矩阵Λ%表示k-1时刻到k时刻的状态预测值。
[0043] 采用上述技术方案所产生的有益效果在于：
[0044] 本发明采用抗噪性能较好的扩展卡尔曼滤波方法对电容层析成像进行图像重建， 在建立系统方程时加入系统噪声，并在线修正系统噪声的协方差矩阵，使状态空间模型更 能反映真实的多相流运动，提高图像重建精度。
【附图说明】
[0045] 图1是本发明的流程图。
【具体实施方式】
[0046] 下面结合附图和【具体实施方式】对本发明作进一步详细的说明。
[0047] 一种基于自适应扩展卡尔曼滤波的ECT图像重建方法，包括以下步骤：
[0048] 步骤1:建立系统方程和测量方程，由以下具体步骤组成：
[0049] 步骤1-1:建立系统方程：
[0050] gk = gk-i+?k-i (1)
[0051] 式中，gk和gk-^别代表k时刻和k_l时刻的图像灰度值；为k_l时刻的系统噪 声，所述系统噪声为白噪声，其均值为〇，协方差矩阵为Qk;
[0052] 所述系统方程利用系统噪声反映多相流流型变化，在建立系统方程时，加入系统 噪声，以反映当流体流型发生变化时，系统模型的改变，其噪声方差阵未知，需要进行估计； [0053] 步骤1-2:建立非线性方程：
[0054] Vk = Uk(gk)+vk (2)
[0055] 式中，Vk为k时刻的电容测量值;Uk(gk)表示图像灰度电容测量值与之间的非线性 关系;vk为k时刻的测量噪声，所述测量噪声为白噪声，其均值为0,协方差矩阵为R;
[0056] 所述非线性方程用于描述电容测量值与图像灰度关系；电容测量值与图像灰度值 之间为非线性关系，利用测量噪声表示测量过程中产生的噪声，认为测量噪声的方差矩阵R 在测量过程中不变；
[0057]步骤1-3:对所述步骤1-2建立的非线性方程进行泰勒级数展开：
[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向量，误差协方差矩阵初值Ρο为α · I， I为单位矩阵，α为比例系数。
[0066]步骤2-2:计算系统噪声的协方差矩阵Qk的更新方程：
[0067] (6)
[0068] 式中，?表示新息序列％ =? 的方差矩阵，zk表示k时刻的归一化测量值； 表示从时刻k_l到时刻k的一步量测预测值;Kk为滤波增益矩阵;Pk和Ph分别表示时刻k 和时刻k-1的滤波误差协方差矩阵；
[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时刻的一步预测协方差;〇!^表 不k_l时刻系统噪声协方差矩阵；
[0072]本发明基于极大似然准则，对系统噪声的协方差矩阵Qk进行实时修正；
[0073] 步骤3 :判断是否达到预设的迭代终止条件Ih S，在本实施例中ε =〇 . 001， 如果是，转向步骤4;否则转向步骤2;
[0074] 步骤4:输出图像灰度的最优估计值。
[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时刻的状态转移矩阵，毛/w表示k-1时刻到k时刻的状态预测值。
[0084] 所述扩展卡尔曼滤波方程基于误差方差最小的准则得到。
【主权项】
1. 一种基于自适应扩展卡尔曼滤波的ECT图像重建方法，其特征在于： 包括W下步骤： 步骤1:建立系统方程和测量方程，由W下具体步骤组成： 步骤1-1:建立系统方程： 阱= gk-l+?k-l (1) 式中，gk和gk-1分别代表k时刻和k-1时刻的图像灰度值；ω k-1为k-1时刻的系统噪声，所 述系统噪声为白噪声，其均值为0,协方差矩阵为Qk; 步骤1-2:建立非线性方程： Vk = Uk(阱）+vk (2) 式中，Vk为k时刻的电容测量值;化(gk)表示图像灰度电容测量值与之间的非线性关系； vk为k时刻的测量噪声，所述测量噪声为白噪声，其均值为0，协方差矩阵为R; 步骤1-3:对所述步骤1-2建立的非线性方程进行泰勒级数展开： 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向量，误差协方差矩阵初值Ρο为α · 1，1为 单位矩阵，α为比例系数。 步骤2-2:计算系统噪声的协方差矩阵Qk的更新方程：式中，每表示新息序列呢=毎-％…的方差矩阵，Zk表示k时刻的归一化测量值A…表 示从时刻k-1到时刻k的一步量测预测值;Kk为滤波增益矩阵;Pk和Pk-i分别表示时刻k和时刻 k-1的滤波误差协方差矩阵； 步骤2-3:扩展卡尔曼滤波：(7) 式中，么和為_1分别表示时刻k和时刻k-1的图像灰度的估计值;么表示从时刻k到时 亥化-1图像灰度的一步预测值;Pk/k-i表示从k-1时刻到k时刻的一步预测协方差;Qk-i表示k-1 时刻系统噪声协方差矩阵； 步骤3:判断是否达到预设的迭代终止条件I私-耗为终止阔值，如果是，转向步 骤4;否则转向步骤2; 步骤4:输出图像灰度的最优估计值。 所述步骤2-2中序列新息序列成=? ~S.w_，的方差矩阵为：(8) 所述步骤2-2中序列新息序列％ =句-^;._1的方差矩阵的计算方法为：(9) 式中，电容测量值Zk表示为： Zk =化Hi+Vk (10) 一步测量预测值表示为： =巧A-/1--1 U1) 式中，化表示k时刻的状态转移矩阵，爲康示k-1时刻至化时刻的状态预测值。
【文档编号】G06T5/10GK106097285SQ201610575305
【公开日】2016年11月9日
【申请日】2016年7月21日
【发明人】张立峰
【申请人】华北电力大学（保定）