超声相控阵R-Θ线性插值扫描转换与双线性插值实现

简介：这份PPT面向医学影像、超声成像与信号处理方向的学生和研究人员，围绕R-Θ线性插值方法讲解超声波图像重建的完整思路，适合用于课题汇报、课程展示或实验复盘。压缩包内共1个ppt文件，体积约1.17MB，内容以原理图、公式推导和实验参数表格为主，便于直接引用到汇报材料中。已有208人学习下载。内容覆盖实验B超成像的波束形成、图像存储、坐标变换与DSC部件流程，并给出16倍数据抽取、探头35C50HA的3.5MHz标准频率、40MHz采样频率、50mm半径、128阵元与0.498mm阵元间距等关键参数，以及240条扫描线、68度扇扫角度下二进制裸数据的提取方法。RF信号部分涉及时域TGC补偿效果、频域中心频率识别，以及移频、滤波、抽取构成的数字下变频解调过程，读者可据此理解坐标变换与插值如何提升图像连续性和分辨率，为超声图像质量优化提供参考。

1. 从扇形扫描说起：为什么超声重建离不开 R-Θ 线性插值

超声相控阵采回来的原始数据，长得像一张按角度排列的表格：第 i 行是探头在偏转角 θ_i 上的一条 A 扫，第 j 列对应深度 r_j，整块数据躺在 R-Θ 极坐标里。人眼要的是屏幕上那张扇形切片图，把极坐标网格换成直角坐标网格的这一步，业内叫扫描转换，R-Θ 线性插值就是其中最常用的一把刀。反直觉的地方在于，重建图像里那些放射状的亮暗条纹，多半不是探头坏了，也不是噪声大，而是 θ 方向采样太稀、插值核选得太糙，甚至是在射频域而不是包络域做的插值。这份内容适合三类人：刚拿到极坐标数据不知道怎么显示的入门者、要把超声波重建图象塞进实时流水线的工程师，以及要做课题汇报、需要把整个重建流程讲清楚的人。下面从采集链路一直讲到映射表固化，每一步都给可跑的命令和参数。

2. R-Θ 极坐标数据是怎么来的：采集链路与重建网格的对应关系

2.1 相控阵与环阵的极坐标采样，以及 (nθ, nR) 数据矩阵

不管是线阵做虚拟阵元偏转，还是凸阵、环阵本身带弧度，采集端形成的都是一组「角度 + 深度」的样本。每一次发射接收事件对应一个偏转角 θ_i，沿声轴方向按采样率 fs 采一条 A 扫，得到 nR 个点；扫完 nθ 个角度，就得到形状为 (nθ, nR) 的矩阵。这个矩阵的 dtype 值得讲究：原始射频信号通常是 int16，做包络检波（Hilbert 变换取模）和对数压缩之后，用 float32 参与插值最省心，一旦中途落到 uint8，双线性插值的四个权重相乘会把量化误差放大成可见的台阶。

这里有个决定成败的顺序问题：插值必须在包络检波之后做，不能在射频域做。射频信号是双极性的，过零点附近数值剧烈翻转，双线性插值本质上是四个点的加权平均，在正负交替的过零点上平均会得到接近零的值，重建出来就是一排排伪条纹，看起来像干扰，实际是数学造成的。

2.2 为什么用反向映射，而不是把极坐标样本打到直角网格上

正向映射是把每个极坐标样本按坐标变换打到直角像素上，近场角分辨率富余，一个像素被反复覆盖，远场则出现大量空洞，除非再做一轮散点填充或者形态学补洞，工程上极少这么干。反向映射反过来：遍历输出直角网格上的每个像素，反算它在极坐标里的位置，再做插值。好处是每个输出像素一定有一组确定的 (i, j, a, b)，不会出现空洞，代码可以用纯向量化写出来，没有循环分支。

代价是极坐标数据里有些区域永远不会被访问到，比如扇形之外的部分，这属于可接受的浪费。另一个常见误用是先把极坐标数据沿 r 方向重采样到和输出像素一样密，再直接取整索引——这等于放弃了插值，径向会出现明显的阶梯状条纹。

2.3 重建网格参数表与最小可跑的参数计算

import numpy as np c = 1.54 # 声速，单位 mm/μs，配合 MHz 采样率使用 fs = 40.0 # 采样率，单位 MHz，即 samples/μs dr = c / (2 * fs) # 径向步长，约 0.0193 mm r0 = 0.0 n_r = 1024 n_th = 128 th0, th1 = np.deg2rad(-45.0), np.deg2rad(45.0) dth = (th1 - th0) / (n_th - 1) # 复现极坐标采集网格 r_axis = r0 + np.arange(n_r) * dr # (n_r,) 深度轴 th_axis = th0 + np.arange(n_th) * dth # (n_th,) 角度轴 # 输出直角网格 dxy = 0.3 # mm/pixel Rmax = r0 + (n_r - 1) * dr # 视场半径取数据满量程 W = int(2 * Rmax / dxy) + 1 H = int(Rmax / dxy) + 2 apex = (H - 1, W // 2) # 顶点放底边中点 print(f"dr={dr:.4f} mm, W={W}, H={H}, apex={apex}")

这段代码解决的是「两套网格怎么对齐」的问题。dr由声速和采样率决定，是不能随手改的物理量，想提高径向分辨率只能提高采样率或者做插值。dxy是显示参数，可以自由选，选得比dr小意味着输出像素比原始样本还密，插值会变得平滑但不会凭空多出信息。apex必须和真实探头阵元的等效相位中心对上，凸阵和线阵偏转扫描的顶点位置差别很大，拿不准的时候用点目标标定一下。

3. 用 NumPy 写出第一版 R-Θ 双线性插值扫描转换

3.1 反向映射的坐标变换与索引公式

设输出像素坐标为 (u, v)，u 向右，v 向下，物理尺寸换算到毫米：

x = (u - cx) · Δx y = (cy - v) · Δx r = √(x² + y²) θ = atan2(x, y)

角度用 atan2(x, y) 而不是 atan2(y, x)，是把声轴（深度方向）当成 0 角，横向偏移决定正负。索引换算成：

f_i = (θ - θ0) / Δθ f_j = (r - r0) / Δr

再取 i = ⌊f_i⌋、j = ⌊f_j⌋，小数部分 a = f_i - i、b = f_j - j，双线性输出为 (1-a)(1-b)·P[i,j] + a(1-b)·P[i+1,j] + (1-a)b·P[i,j+1] + ab·P[i+1,j+1]。

3.2 向量化双线性插值的完整实现

import numpy as np def scan_convert(polar, r0, dr, th0, dth, out_hw, dxy, apex=None): """ polar : (n_th, n_r) 已做包络检波与对数压缩的极坐标数据 r0/dr : 径向起点与步长（mm） th0/dth: 角度起点与步长（rad） out_hw : 输出图像 (高, 宽) dxy : 输出像素间距（mm/pixel） apex : 扇形顶点在图像中的 (行, 列)，默认底边中点 """ n_th, n_r = polar.shape H, W = out_hw cy, cx = apex if apex is not None else (H - 1, W // 2) # 1) 输出像素 -> 物理坐标（mm） u = np.arange(W, dtype=np.float32) v = np.arange(H, dtype=np.float32) xx, yy = np.meshgrid((u - cx) * dxy, (cy - v) * dxy) # 2) 直角坐标 -> 极坐标 r = np.hypot(xx, yy) th = np.arctan2(xx, yy) # 3) 物理坐标 -> 极坐标索引（浮点） fi = (th - th0) / dth fj = (r - r0) / dr # 4) 取整、取小数部分 i0 = np.floor(fi).astype(np.int32) j0 = np.floor(fj).astype(np.int32) a = (fi - i0).astype(np.float32) b = (fj - j0).astype(np.float32) # 5) 越界掩膜，同时把索引夹到合法范围，避免读越界 valid = (i0 >= 0) & (i0 < n_th - 1) & (j0 >= 0) & (j0 < n_r - 1) i0c = np.clip(i0, 0, n_th - 2) j0c = np.clip(j0, 0, n_r - 2) f00 = polar[i0c, j0c] f10 = polar[i0c + 1, j0c] f01 = polar[i0c, j0c + 1] f11 = polar[i0c + 1, j0c + 1] out = (f00 * (1 - a) * (1 - b) + f10 * a * (1 - b) + f01 * (1 - a) * b + f11 * a * b) out[~valid] = 0.0 return out.astype(np.float32)

四个数组f00/f10/f01/f11用花式索引一次性取出形状为 (H, W) 的四张图，再做逐元素加权，整段没有 Python 循环，1280×800 的图在普通笔记本上几十毫秒能跑完。参数上的关键点有三个：np.clip必须在索引之前做，否则i0c + 1会越过数组边界；valid掩膜用逻辑与串起来，一定要包含+1之后的边界，只判i0 >= 0是不够的；输出的 0 值代表扇形之外，如果后面要做对数显示，记得先加一个小量再取对数，不然 log(0) 会给出 -inf。

3.3 同一件事交给 cv2.remap 和 LUT 做

import cv2 # 复用 3.2 中的 fi/fj，注意 remap 的 map1 是列索引、map2 是行索引 map_x = fj.astype(np.float32) # 对应 polar 的列，即 r 方向 map_y = fi.astype(np.float32) # 对应 polar 的行，即 θ 方向 img = cv2.remap(polar.astype(np.float32), map_x, map_y, interpolation=cv2.INTER_LINEAR, borderMode=cv2.BORDER_CONSTANT, borderValue=0.0)

cv2.remap的参数顺序容易踩坑：map_x是极坐标矩阵的列索引（半径方向），map_y是行索引（角度方向），和我们平时说的 x/y 是两回事，写反了图像会被转置。borderMode选常量填充 0，越界像素自动变黑，等于自带扇形掩膜，不用再手写valid。interpolation换成cv2.INTER_CUBIC就是三次插值，代价大约是双线性的四倍访存。

3.4 扇形掩膜与径向截止

除了越界，还有两个区域需要显式处理：Rmax 之外没有数据，即使索引在数组范围内也要置零；扇形开角之外的像素同理。最省事的做法是在映射阶段就把r > r_max的像素写进掩膜。

mask = (r <= r0 + (n_r - 1) * dr) & (np.abs(th - (th0 + (n_th - 1) * dth / 2)) <= (n_th - 1) * dth / 2) out = np.where(mask, out, 0.0)

如果后续要做图像拼接或者定量测量，建议把掩膜单独存一份布尔数组，而不是依赖 0 值来判断背景，因为真实的超声回波在小深度处也可能接近 0。

4. 参数调不好就出伪影：R-Θ 插值的误差来源与排错清单

4.1 远场横向欠采样，是扇形条纹的唯一主因

相邻两条扫描线在深度 R 处的弧长是 R·Δθ，这就是该深度处的横向采样间隔。R = 100 mm、Δθ = 0.7° = 0.0122 rad 时，弧长约 1.22 mm；而输出像素间距 Δx 取 0.3 mm，Nyquist 要求采样间隔不大于 0.6 mm。差了两倍，混叠就出现了，表现为从顶点向外发散的放射状明暗条纹，越远越明显。解决方向只有两个：把 Δθ 采小，或者在 θ 方向先插值加密。前者受脉冲重复频率和帧率约束，后者是纯计算，实际项目里多数选后者。

4.2 插值核怎么挑：一张对照表

选择逻辑很直白：如果下游是给医生看的实时画面，双线性足够，糊一点比振铃好看；如果是量点目标的主瓣宽度、做分辨率标定，就得用三次，否则量出来的宽度被双线性抹平，会偏大百分之十几。

4.3 先在 R-Θ 域上采样，再做扫描转换

这是抗混叠的正确顺序。先把 (nθ, nR) 沿 θ 方向加密，再做反向映射，等效于把 Δθ 降下来。

from scipy.ndimage import zoom # 沿角度方向加密 4 倍，径向保持原样 polar_up = zoom(polar, (4, 1), order=3, mode='nearest', prefilter=True) dth_up = dth / 4 # 后续 scan_convert 用这个新步长 # 如果担心三次插值的振铃，可以先在 θ 方向做一次短窗低通 kernel = np.array([0.25, 0.5, 0.25], dtype=np.float32) polar_lp = np.apply_along_axis(lambda m: np.convolve(m, kernel, mode='same'), 1, polar)

zoom的第二个参数是各维缩放因子，(4, 1)表示只加密角度维；order=3用三次样条；prefilter=True是样条预滤波，关掉它结果会明显偏软；mode='nearest'处理边界，用零边界会把扇形最外侧两条线拉黑。代价是内存涨 4 倍、插值耗时涨若干倍，但相比把硬件角度步长降到四分之一（帧率直接掉到四分之一），这笔账怎么算都划算。

4.4 重建出问题时的排查顺序

1. 先用最近邻跑一遍。如果最近邻下扇形轮廓正确、只是有马赛克，说明坐标变换对了，问题在插值；如果轮廓就是歪的，先回去查顶点和角度定义。
2. 检查图像是否左右镜像。把数据矩阵的行序倒过来跑一次，形状对上了就是行序约定不符。
3. 检查是否有斜向的直线状伪影。这通常是在射频域插值造成的，把插值挪到包络检波之后。
4. 检查是否有大片 NaN 或 -inf。对数压缩在 0 值上取对数，先加 1e-6 的底噪再压缩。
5. 检查远场条纹。按 4.1 的公式算一下 R·Δθ 和 2Δx，超了就是欠采样，不是 bug。

5. 从点目标验证到课题汇报：R-Θ 重建结果的自检与演示技巧

5.1 用点目标 PSF 定量验证重建质量

仿真一个点散射子，按 2.3 的网格生成极坐标数据，走完包络检波和扫描转换，然后过点取剖面量主瓣宽度。

import numpy as np def profile_width(img, row, col, thr=0.5, axis=1): """量 -6dB 主瓣宽度（像素）。img 为线性幅度图，thr=0.5""" line = img[row, :] if axis == 1 else img[:, col] peak = line.max() idx = np.where(line >= peak * thr)[0] return (idx[-1] - idx[0] + 1) if len(idx) else 0

阈值取 0.5 是因为线性幅度图的 -6dB 点正好等于峰值的一半，20·log10(0.5) ≈ -6 dB。如果图已经做过对数压缩，阈值要改成peak - 6再反算回线性，直接套 0.5 量出来的宽度会明显偏大。除了主瓣宽度，还可以沿 θ 方向对图像做一次一维 FFT，如果谱上出现与 Δθ 对应的尖峰，就是角度欠采样的直接证据，这一招在汇报里比嘴上说「感觉有伪影」有说服力得多。

5.2 映射表固化与逐帧复用

映射表只跟网格参数有关，跟数据无关，所以只要探头开角、深度、输出尺寸不变，算一次就够。

# 把 float32 映射表压缩成 16 位整数 + 小数部分，访存减半 map1, map2 = cv2.convertMaps(map_x, map_y, cv2.CV_16SC2) for frame in stream: # 逐帧只做一次 remap img = cv2.remap(frame, map1, map2, cv2.INTER_LINEAR, borderMode=cv2.BORDER_CONSTANT, borderValue=0.0) yield img

convertMaps把映射量化到 1/32 像素，对 8 位灰度显示来说完全看不出来，但 remap 的访存量减半，实测通常能快 1.5 到 2 倍。注意映射表一旦固化成 16 位，换开角或改深度就必须重算，否则整幅图会错位到无法解释。

一个实用的小技巧：把convertMaps之后的映射表连同 Δθ、Δr、顶点坐标、开角一起序列化存盘（np.save或cv2.FileStorage都行），换机器、换语言重写重建模块的时候直接读表对齐，能省掉一整轮坐标约定的反复对拍。

本文还有配套的精品资源，点击获取