[发明专利]基于马尔科夫链蒙特卡洛多目标追踪的经颅超声脑血管造影超分辨率成像方法

权利要求书

1.一种基于马尔科夫链蒙特卡洛多目标追踪的经颅超声脑血管造影超分辨率成像方法，其特征在于：包括以下步骤：

1)利用低频超声线阵换能器向颅内发射超声平面波，根据回波信号获得连续N帧脑血管超声造影图像B_n(n＝1,...,N)，成像帧数N为80～200；

2)对得到的N帧图像B_n(n＝1,...,N)，依次选取其中连续的N_med帧图像对用时域中值滤波函数进行处理，得到背景图像G_n：

其中，Distance(·)表示像素间距离，i表示像素点横坐标，j表示像素点纵坐标；

3)从图像B_n中减去背景图像G_n得到前景图像F_n；

4)对经过伪迹去除的前景图像F′_n设定阈值f_thr进一步滤波，然后进行局部最大值处理，将得到的μ_n个像素点作为造影微泡的位置，并记录造影微泡位置坐标及对应的速度信息所述阈值f_thr为-18dB～-23dB；

5)将第n帧图像中所有μ_n个像素点的坐标组合，得到将所有F′_n处理得到的Y_n组合，得到

6)在Y中设定距离范围为≤1mm，以各个造影微泡位置为中心，将小于距离范围的其他造影微泡位置预设为中心造影微泡位置可能产生轨迹的链接造影微泡位置；

7)利用马尔科夫链蒙特卡洛多目标追踪对造影微泡轨迹进行推演，所述推演的迭代次数L为1500～4000；

8)对推演得到的轨迹集进行成像；

所述步骤7)具体包括以下步骤：

7.1)根据抽取的推演状态，改变Y中的造影微泡位置，并推演得到轨迹集ω＝[ω₁ω₂…ω_M]，其中M为轨迹总数；

7.2)构造时域迭代状态方程：

其中A为过程矩阵：

Δt＝1/Frt，Frt为成像帧率，成像帧率Frt为500Hz～2kHz，为服从高斯分布的过程方差矩阵；

其中表示轨迹集中造影微泡的坐标信息，表示轨迹集中造影微泡的速度信息，k＝1,2,…K，K为轨迹集中造影微泡个数，n＝1,2,…N，N为成像帧数，角标x和z分别表示横向和纵向坐标；

7.3)根据卡尔曼滤波，预测的造影微泡位置为：

其中，为服从高斯分布的观测方差矩阵，为均匀分布；

7.4)统计计算轨迹终点数目z_t、当前轨迹总数n_t、新生轨迹起点个数a_t、被选中为轨迹的造影微泡位置数d_t，以及虚警目标总数f_t，并计算从上次迭代中衍生的新轨迹数目c_t及未被选中为轨迹的造影微泡位置数g_t：

c_t＝n_t-z_t-a_t

g_t＝n_t-d_t

7.5)计算当前推演状态下的先验概率值：

其中，p_z表示目标终止或消失概率，p_d表示目标被选中为轨迹点的概率，λ_b表示理想轨迹中造影微泡位置新生的概率值，λ_f表示理想轨迹中造影微泡位置新生的虚警概率值；

7.6)假设虚警数目和推演路径均服从均匀分布，计算基于当前推演轨迹集和各轨迹造影微泡位置的似然值：

其中，M为推演轨迹总数，τ_k为当前推演轨迹中造影微泡位置总数，为造影微泡位置的似然值，表示卡尔曼滤波中观测值预测向量，表示卡尔曼滤波中观测方差预测矩阵；

7.7)计算当前推演状态下的后验概率值：

其中，P(Y)表示探测造影微泡位置的先验概率；

7.8)以当前推演轨迹集ω为马尔科夫链的某状态，依据最大后验概率方法计算接受当前推演的概率：

其中，ω′为前次推演轨迹集，ω为本次推演轨迹集；

7.9)产生服从均匀分布的随机数ξ，0ξ1，若ξA(ω′,ω)，则接受本次推演，并用轨迹集ω替换轨迹集ω′；否则，拒绝本次推演，轨迹集ω′保持不变；

7.10)重复步骤7.1)-7.9)，直至完成路径推演迭代次数L。