MATLAB：一种超声图像的运动分析方法（块匹配法）

0 引言

心血管疾病( cardiovascular diseases,CVD)危险因素对居民健康的影响愈加显著,CVD 的发病率仍持续增高。 CVD 给居民和社会带来的经济负担日渐加重,已成为重大的公共卫生问题，已成为危害各国人民健康的三大疾病之一，其患病率和死亡率居于全球之首，中国 CVD 患病率处于持续上升阶段。推算 CVD 现患人数 3. 3 亿,其中脑卒中 1300 万,冠心病1139 万,心力衰竭 890 万,肺源性心脏病 500 万,心房颤动 487 万,风湿性心脏病 250 万,先天性心脏病 200 万, 外周动脉疾病 4 530 万,高血压 2. 45 亿[1]。动脉硬化是CVD的主要原因之一，研究表明动脉硬化的早期表现是血管壁内中膜的变厚，导致动脉壁运动的弹性变差，由于血管壁的二维运动的检测是血管弹性检测的基础。所以检测血管壁的二维运动能为血管弹性的检测提供效信息，从而对CVD的早期诊断起到很大帮助[2]。

超声成像技术凭借其无创、实时、低成本、无电离辐射等独特优势，已广泛渗透于临床医学诊断、生物医学研究、工业无损检测等多个核心领域。在临床医学中，超声图像不仅可用于脏器形态结构的直观观察，更能通过对组织、血流等运动信息的精准分析，为心血管疾病、肝脏疾病、胎儿发育异常等病症的早期筛查、病情评估及疗效监测提供关键依据。因此，基于超声图像进行血管壁的二维运动（如图1）是更具优势。基于超声图像运动分析方法，根据核心原理主要基于光流场的方法及基于块匹配的方法等，其中，块匹配法（Block Matching Method, BMM）凭借其原理简洁、计算效率高、抗噪声能力较强、易于工程实现等突出优势，成为当前超声图像运动分析中应用最广泛、最成熟的方法之一。块匹配法的核心思想是将参考帧超声图像划分为若干互不重叠或部分重叠的矩形图像块，通过在目标帧的对应搜索区域内寻找与参考块灰度相似度最高的匹配块，进而确定各图像块的位移矢量，最终实现整个图像序列的运动分析。块匹配法无需复杂的预处理和迭代过程，能在保证一定分析精度的前提下，满足超声图像运动分析对实时性的要求，尤其适用于临床实时诊断、动态监测等对时间响应速度敏感的场景。因此，本文主要介绍块匹配算法的原理以及基于MATLAB如何运行。

注：这个B超图像来自于[3]。

0 引言

1 块匹配法的匹配准则

2 径向运动的估计（MATLAB）

2.1 NCC 匹配法则

2.2 MAD 匹配法则

2.3 MSD 匹配法则

3 结果分析

3.1 时间的差异

3.2 径向位移的差异

4 讨论

5 参考

1 块匹配法的匹配准则

块匹配算法将图像看作是由一个个块组成，且认为各个块中的所有像素的光流是一样的。在下一帧(comparison image)中的一定范围(search range)按照一定的匹配准则进行搜索，寻找当前帧(reference image)的某个块(模板块,template block)的最相近的对应块(候选块,candidate block)，前后两个对应的块的位移向量就作为该块的运动向量[4]，如下图2所示：

在块匹配过程中，用于衡量两图像块匹配程度的量化指标种类繁多，其中应用较为广泛的包括归一化相关系数（NCC, Normalized Cross Correlation）、平均均方误差（MSD, Mean Square of Differences）以及平均绝对差（MAD, Mean Absolute of Differences）等。在这三类指标中，归一化相关系数（NCC）与平均绝对差（MAD）因各自的特性优势，在实际应用中更为普遍。针对序列超声图像中待匹配的M×N尺寸图像块X与Y，上述三种匹配度衡量指标的定义分别如下[4]：

其中，。分别表示前后两帧超声图像中M×N尺寸（见图2）的第m行和第n列的像素值。

2 径向运动的估计（MATLAB）

现有研究大多数都是针对径向运动进行估计[5-8]，因此这里仅对径向运动进行分析，使用上述三种匹配法则（NCC、MAD和MSD）分别进行估计。

2.1 NCC 匹配法则

为实现上述 NCC 公式的计算（选择最大值），使用 MATLAB 编写的核心函数如下：

function rho = calc_corr(cw1, cw) % 二维互相关公式 % cw1 -匹配块 % cw - 参考块 %计算cw和cw1 的均值矩阵 cw1_mean = mean(cw1(:)) * ones(size(cw1)); cw_mean = mean(cw(:)) * ones(size(cw)); %计算cw和cw1与均值矩阵的差值 term1 = cw1 - cw1_mean; term2 = cw - cw_mean; % 对差值相乘求和 numerator = sum(term1(:) .* term2(:)); % 对差值先平方 再求和 然后相乘，再开方 denominator = sqrt( sum(term1(:).^2) * sum(term2(:).^2) ); % 计算相关系数 rho = numerator / denominator; end

本方法采用 10×10 像素的图像模块进行运动特征分析，固定第一帧为参考基准帧，将其余所有帧作为待匹配帧依次完成匹配运算，最终输出各帧相对参考帧的位移量，以及整个分析过程的运行耗时。结果如下：

2.2 MAD 匹配法则

为实现上述 MAD 公式的计算(选择最小值），使用 MATLAB 编写的核心函数如下：

function mad_val = calc_mad(cw1, cw) % 二维平均绝对差（MAD）公式 % cw1 - 匹配块 % cw - 参考块 % 计算两个块对应元素的绝对差 abs_diff = abs(cw1 - cw); % 计算所有绝对差的平均值 mad_val = -mean(abs_diff(:)); end

本方法同上采用 10×10 像素的图像模块进行运动特征分析，固定第一帧为参考基准帧，将其余所有帧作为待匹配帧依次完成匹配运算，最终输出各帧相对参考帧的位移量，以及整个分析过程的运行耗时。结果如下：

2.3 MSD 匹配法则

为实现上述 MSD 公式的计算(选择最小值），使用 MATLAB 编写的核心函数如下：

function msd_val = calc_msd(cw1, cw) % 二维均方差（MSD）公式 % cw1 - 匹配块 % cw - 参考块 % 计算两个块对应元素的差值平方 squared_diff = (cw1 - cw).^2; % 计算所有平方差的平均值 msd_val = -mean(squared_diff(:)); end

三种匹配法则的代码如下：

clc; clear; % 启动计时器 tic; load("B1.mat"); num_frames = size（data, 3); r = zeros(1,num_frames); c = zeros(1, num_frames); % 确定参考帧和匹配帧的尺寸大小 %参考帧：cw起始行a 终点行b 起始列d 终点列e %匹配帧：cw1起始行a1 终点行b1 起始列d1 终点列e p=10;%最大像素位移 a=321; b=330; d=41; e=50; a1=a-p;b1=b+p;d1=d-p;e1=e+p; for k = 1:num_frames cw = data(a:b, d:e, 1);%参考帧 cw_t =data(a1:b1, d1:e1, k);%匹配帧 R = zeros(2*p+1, 2*p+1); % 计算相关系数矩阵 for i = 0:2*p for j = 0:2*p %R(i+1,j+1) = calc_corr(cw_t(i+1:i+10,j+1:j+10), cw);#调用NCC %R(i+1,j+1) = calc_mad(cw_t(i+1:i+10,j+1:j+10), cw);#调用MAD R(i+1,j+1) = calc_msd(cw_t(i+1:i+10,j+1:j+10), cw);;#调用MSD end end [sorted_values, sorted_indices] = sort(R(:), 'descend'); [rows, cols] = ind2sub(size(R), sorted_indices); found = false; for idx = 1:length(sorted_values) current_row = rows(idx); current_col = cols(idx); row_indices = current_row; col_indices = current_col; end r(k) = row_indices; c(k) = col_indices; end pluse=-(r-(p+1)); plot(pluse); elapsed_time = toc; % 输出运行时间 fprintf('程序总运行时间：%.4f 秒\n', elapsed_time);

3 结果分析

3.1 时间的差异

从运行时间来看，NCC (0.4358s) > MSD (0.3781s) > MAD (0.3428s)，这一结果直观反映了三种相似性度量算法的计算效率差异：（1）计算复杂度决定耗时差异：NCC 耗时最长，核心原因是其需要先计算两个块的均值矩阵、差值矩阵，再完成乘法求和与开方运算，涉及均值计算、矩阵减法、逐元素乘法、平方求和、开方多步操作；而 MSD 仅需逐元素平方差 + 均值，MAD 仅需逐元素绝对差 + 均值，运算步骤更少、浮点操作数更低，因此耗时依次降低。（2）效率与精度的取舍：虽然 MAD 耗时最短，但需注意算法特性的适配性 ——NCC 因做了均值归一化，对亮度偏移、对比度变化的鲁棒性最强，适合光照不稳定的场景；MSD 对像素差值的平方放大效应使其对大误差更敏感，适合需要精准检测像素偏差的场景；MAD 计算最快且对异常值鲁棒性优于 MSD，适合实时性要求高的轻量化场景。

3.2 径向位移的差异

从图6中可见，基于 MAD 计算的径向位移在约 15–20 区间出现了明显的峰值，数值突破 6，显著高于 NCC（图中未出现明显峰值）和 MSD 的同期最高值（约 5）。尽管 MAD 的波动幅度更大，但三者的整体变化趋势基本保持同步：在初期阶段均呈快速上升态势，随后在中期进入平台期并伴随阶梯式下降，最终在末期共同趋于零值。这种趋势上的一致性，表明三种算法在反映径向位移的整体演化规律上具有良好的一致性，而 MAD 对局部突变的高敏感性。

代码如下：

clc; clear; NCC=load("pluse_NCC.mat");MAD=load("pluse_MAD.mat");MSD=load("pluse_MSD.mat"); ncc = NCC.pluse; mad = MAD.pluse; msd = MSD.pluse; figure('Position', [100, 100, 1400, 800], 'Color', 'white'); plot(ncc, 'r-', 'LineWidth', 1.2, 'DisplayName', 'NCC'); hold on; plot(mad, 'g-', 'LineWidth', 1.2, 'DisplayName', 'MAD'); plot(msd, 'b-', 'LineWidth', 1.2, 'DisplayName', 'MSD');

4 讨论

本文主要介绍块匹配算法的原理以及基于MATLAB如何运行。从算法原理层面分析，MAD作为一种基于像素灰度值绝对误差的匹配准则，其计算逻辑为逐像素求取参考块与待匹配块的灰度差值绝对值并取平均，这种计算方式对像素灰度的局部突变更为敏感，因此在径向位移峰值处表现出更高的数值；而 NCC通过对灰度值进行归一化处理，有效降低了光照变化等因素的干扰，MSD则通过平方放大了误差的权重，二者对局部极值的响应弱于 MAD，这也是三者峰值差异的原因。其中 MAD 函数的核心代码仅需一行即可实现，计算效率最高；NCC 函数则需先通过cov函数计算协方差，再结合标准差完成归一化，代码复杂度稍高；MSD 函数需通过完成均方差求解。实际测试中，在 MATLABR2023a 环境下，对斑块大小（10X10）进行块匹配计算时，MAD 算法的单次运行耗时约 0.34s，NCC 约 0.44s，MSD 约 0.38s，这与算法原理的计算复杂度完全契合。