颈动脉管壁搏动位移超声检测

• 实验目的

1、掌握血管壁搏动位移的超声检测方法；

2、初步掌握信号相似性的评价方法；

3、初步掌握信号滤波的方法。

• 实验内容

1、简述基于超声射频信号的管壁搏动位移及脉搏波速超声测量过程。

2、 基于数据sim_bmd6_45，使用归一化相关系数为管壁追踪的代价函数，测量13线空间位置的管壁搏动位移曲线。

3、 对第2题中获得的管壁搏动位移曲线进行滤波，消除台阶状的量化误差。

4、 使用第3题中已滤波的13线管壁搏动位移曲线估计脉搏波速（PWV）。

必要的实验参数：

        1）13线空间位置平均分布在38mm的颈动脉上。

        2）每线空间位置上，有1333个射频信号。

        3）1333个射频信号的采样时间为1秒。

        4）射频信号的采样频率为100MHz。

• 方法步骤

1、超声探头对颈部血管位置扫描时发射13个声束（如图1所示），这时我们取其中一个位置来分析。 

图 1

 图 2

        图2中显示的是第一个位置（声束）的成像，可以看到这个声束中不同扫描深度处的成像。由于扫描时对每一个声束进行了1333次超声扫描，所以可以将1333次扫描看作是时间轴。我们用第一次扫描结果中的血管壁位置（我选的是左边的血管壁）与其它 1332次扫描结果进行相关系数测量，然后找到相关系数最大的位置，用这个位置减起始位置得到像素位移（以像素个数为单位），再将像素位移转换为实际位移（以mm为单位），得到的结果即为这个位置的管壁搏动位移。同理，计算得到13个位置的管壁搏动位移，观察它们相同幅度在时间上的延时dt，然后再计算得到13个测量位置dx，计算dx/dt即可计算出脉搏波速pwv。

1. 管壁搏动位移曲线测量

1. 1参数选择

up=3101;
down=3320;
sampling_width=down-up;%3100~3320
number_k=1333;
pixed_width=1540/10^8*1000/2;

       这里我选取的测量位置是图2中左边血管壁的位置，在3100—3320处。sampling_width是采样宽度，即血管壁宽度，pixed_width是每个像素实际表示多少mm的宽度，这里除2是因为超声成像是回声成像。

1. 2数据读取

for f=1:13

url="line_noise"+num2str(f);

load(url,'-mat')

data=RFline_noise;

 这里用循环读入13组数据。

1. 3位移

for m=1:number_k

center=data(up:down,1);

search=data(up-300:down+300,m);

for i=1:300*2-sampling_width

        ncc(i)=corr2(center,search(i:i+sampling_width,:));

end

l=300-find(ncc==max(ncc))+1;

distanse(m,f)=l*pixed_width;

end

       这里外循环用于计算13个位置中每一个位置的位移，内循环用于在一个位置中计算第一次扫描结果与剩下的1332次扫描结果的相关系数。具体实现步骤是：以第一次扫描结果为中心块，中心块左右300的位置为搜索块，在搜索块中找到与中心块最匹配的位置，然后计算像素位移再转换为实际位移distanse，这里用300减最匹配位置是因为脉搏跳动的时候血管壁要向up一侧扩张。

1. 4数据保存

       由于上述过程使用了三个循环嵌套，所以shij2复杂度非常高，为了后续计算速度能更快，所以要保存上述的数据。

1. 5绘图

x_axis=(0:1332).*ones(13,1);

x_axis=x_axis';

y_range=(0:12)*38/12;

y_axis=y_range.*ones(1333,1);

figure

plot3(x_axis,y_axis,distanse,'Color',[0 0 0]),hold on,grid on

ylabel('Position(mm)'),xlabel('Time(ms)'),zlabel('Distension(mm)')

       在这一部分我计算了13个位置的x坐标，即时间坐标，y坐标。然后用黑色的线条绘制了13条管壁搏动位移曲线。

2. 滤波、消除台阶状的量化误差

wp=6/1333;

ws=10/1333;%设置阻带频率

Rp=0.5;%设置通带波纹系数

Rs=1.5;%设置阻带波纹系数

       这里我设计了低通滤波器对13个位置的脉搏位移曲线进行了滤波，这里通带频率为6Hz，阻带频率为10Hz，通带波纹系数为0.5，阻带衰减系数为1.5。其中对通带频率和阻带频率进行了归一化处理。

[N,Wn]=buttord(wp,ws,Rp,Rs);

[b,a]=butter(N,Wn);

filter_datas(:,f)=filter(b,a,distanse(:,f));

  使用巴特沃斯滤波器对每一条曲线进行滤波。

3. 估计脉搏波速（PWV）

s0=0.06;

for i=1:13

dt0=find(filter_datas(:,i)>s0);

   dts=dt0(1);

   dt(i)=dts;

end

dt(2:13)=dt(2:13)-dt(1);

dt(1)=0;

ds=0:38/12:38;

cftool(dt,ds)

        这里我找了大于0.06的第一个位置，然后在13条曲线中找到这个位置对应的时间即为dt，对dt和ds进行数据拟合，斜率即为脉搏波速PWV。

• 实验结果

       对第一个位置进行管壁搏动位移曲线测量得到的曲线如图3所示，可以看到搏动最大位移为0.4mm，曲线存在阶梯状的量化误差。

图 3

图 4

       通过计算得到13条管壁搏动位移曲线如图4所示，图中的曲线存在阶梯状的量化误差，为了消除这个误差我对曲线进行了滤波处理。

图 5

       通过滤波得到13条管壁搏动位移曲线如图5所示，对比曲线的趋势可以初步判断这是中年组的管壁搏动位移曲线，所以测量得到的PWV应该为6.0m/s左右。

图 6

       拟合ds和dt得到的结果如图6所示，观察拟合曲线的斜率可以得到PWV为6.378m/s，这个结论和上面的预期一致，是中年组的PWV。

完整代码：

%%
clear
clc
up=3101;
down=3320;
sampling_width=down-up;%3100~3320
number_k=1333;
wp=6/1333;
ws=10/1333;
Rp=0.5;
Rs=1.5;
pixed_width=1540/10^8*1000/2;
for f=1:13
    url="line_noise"+num2str(f);
    load(url,'-mat')
    data=RFline_noise;
    for m=1:number_k
        center=data(up:down,1);
        search=data(up-300:down+300,m);
        for i=1:300*2-sampling_width
            ncc(i)=corr2(center,search(i:i+sampling_width,:));
        end
        l=300-find(ncc==max(ncc))+1;
        distanse(m,f)=l*pixed_width;
    end
    [N,Wn]=buttord(wp,ws,Rp,Rs);
    [b,a]=butter(N,Wn);
    filter_datas(:,f)=filter(b,a,distanse(:,f));
end
save distanse
save filter_datas
save distanse
%%
close all
load filter_datas.mat
load distanse.mat
x_axis=(0:1332).*ones(13,1);
x_axis=x_axis';
y_range=(0:12)*38/12;
y_axis=y_range.*ones(1333,1);
figure
plot3(x_axis,y_axis,distanse,'Color',[0 0 0]),hold on,grid on
ylabel('Position(mm)'),xlabel('Time(ms)'),zlabel('Distension(mm)')

figure
plot3(x_axis,y_axis,filter_datas,'Color',[0 0 0]),hold on,grid on
ylabel('Position(mm)'),xlabel('Time(ms)'),zlabel('Distension(mm)')
s0=0.06;
for i=1:13 
    dt0=find(filter_datas(:,i)>s0);
    dts=dt0(1);
    dt(i)=dts;
end
dt(2:13)=dt(2:13)-dt(1);
dt(1)=0;
ds=0:38/12:38;
cftool(dt,ds)
%view([0 0])%正视图
%%