[发明专利]一种基于序贯蒙特卡洛原理的EMT图像重建方法

权利要求书

1.一种基于序贯蒙特卡洛原理的EMT图像重建方法，其特征在于，所述方法包括以下步骤：

步骤一，EMT图像重建过程建模：EMT图像重建的本质是通过检测电压和系统灵敏度场分布情况反向推演出物场电导率分布情况，三者的关系用下式表示：

V＝Sg (1)

其中，V是28×1的独立电压向量，由传感器阵列采集和电子控制单元处理后得到；S是28×1600的系统灵敏度矩阵，反映了物场中各个位置的电导率和电压的对应关系；g是1600×1的电导率向量，其分布情况可以代表物场中目标物体的位置与形状和管道周围磁性介质的分布情况；

由于EMT系统求解原始图像信号是一个病态和欠定的问题，将序贯蒙特卡洛方法引入EMT图像信号的求解过程；对于两相流参数检测问题，原始图像为二值图像，灰度值为0或1，假设管道内的截面被分割成d×d像素，其中d为整数，其取值由网格剖分数决定，二值图像表示为二维d×d矩阵或一维N×1向量g，其中N＝d×d；

首先利用蒙特卡洛思想将EMT系统建模为一个动态时变系统，将图像重建问题的求解建模为该动态时变系统对最优解的搜索，EMT系统的图像重建过程是利用实测电压数据V根据任意算法反向推演出电导率分布g的测量值，如果电导率分布g的求解是通过递归迭代的方法的，搜索过程的最优解可视为动态时变系统，将迭代次数作为离散的时间量，自变量g表示系统的状态变量，最优解表示动态系统的状态估计，假设一个动态时变系统可以由系统的运动模型和观测模型来描述其数学模型：状态变量g基于迭代次数的变化可视为系统的运动模型；包含观测噪声的被测电压矢量V与系统状态变量g之间的关系方程为系统的观测模型；

系统状态g的更新方程为：

g_k＝g_k-1+u_k-1 (2)

其中k为代表迭代次数的离散时间量；g_k和g_k-1分别为k时刻和k-1时刻对应的电导率分布向量，即反映电导率分布情况的图像信号；u_k-1为过程噪声；

系统的观测方程为：

V_k＝Sg_k+v_k (3)其中，S为灵敏度矩阵，V_k为k时刻的电压矢量，v_k为测量噪声；

在动态时变系统的运动模型中，假设在k时刻的状态估计g_k由上一时刻的状态g_k-1决定，k-1时刻的过程噪声u_k-1分布确定了系统搜索空间，假设过程噪声u_k-1属于高斯噪声，其均值为投影数据的修正量E(u_k-1)＝-γΔf(g_k-1)，其中γ为松弛因子，取一个较小的正数；Δf(g_k-1)是通过任意像素点的所有投影数据的误差值的统计，通过下式得到：

其中diag(SS^T)是由SS^T的对角线分量组成的向量，除法是指分子中向量的任意元素除以分母中向量的对应元素；

公式(5)根据上一次迭代中获得的图像信号g_k-1计算估计电压值V_k-1并将其与实际测量电压V进行差值运算，然后根据1/diag(SS^T)的权重对差值进行加权，最后通过反演得到g_k-1信号的修正值；没有先验知识的变量可以简单设置，所以过程噪声的概率密度函数为u_k-1～N(-γΔf(g_k-1),1)，状态空间的概率密度分布为g_k-1～N(g_k-1-γΔf(g_k-1),1)；

步骤二，初始化状态空间：由于LBP算法计算简单快速，可以用来生成初始状态空间；

g₀＝g_LBP＝S^TV (5)

以LBP算法得到的结果作为初始状态空间；

步骤三，状态空间采样及样本权值更新：得到初始样本空间后，对其进行重要性采样，根据样本的观测距离设置相应的权重；

步骤四，最优状态估计：基于获得的粒子及其权重，计算系统的状态空间估计，并通过迭代更新的方式搜索最优解，最后作为全局最优估计输出进行图像重建。