企业官网建设流程全解析

你有没有想过，CT 图像并不是传统意义上“拍”出来的，反而是“算”出来的？探测器真正接收到的只是 X 射线穿过物体后的强度衰减数据，至于物体内部哪里密度高、哪里密度低，全都要靠数学模型从这些投影数据里反推。CT 技术与图像重建的核心，就是把“重建”变成一个求解代数方程组 Ax=b 的问题，而在迭代求解过程中又天然地引入了差分方程递推的思想。这篇文章我结合一个纯 Python 的模拟实验，把从投影采集、滤波反投影到代数迭代重建的完整链路拆开讲清楚，适合数学建模竞赛选手、图像处理入门者和想理解 CT 重建底层逻辑的开发者参考。

1. 从射线到方程：CT 重建的第一块积木

1.1 投影数据到底是什么

先回到物理层。一束强度为 I₀ 的 X 射线穿过物体后，强度会衰减为 I，这个过程可以用朗伯-比尔定律近似描述：

I = I₀ · exp(-∫ μ(l) dl)

其中 μ(l) 是路径上每个位置的线衰减系数，积分沿射线路径进行。CT 探测器记录的是 I，而我们真正关心的其实是投影值：

p = ln(I₀ / I) = ∫ μ(l) dl

这个 p 是一个标量吗？不是。针对每一个入射角度，探测器上会有成百上千个单元，每个单元都测量一条独立射线的投影值。把所有角度的投影值按顺序拼在一起，就得到一张二维的“正弦图”，英文叫 sinogram。之所以叫这个名字，是因为任意一个内部点在不同角度下的投影位置会画出一条正弦曲线。

所以，CT 重建并不是“看图恢复”，而是要回答这样一个问题：给了我一系列角度下的一维投影值 p，如何反推出二维空间里每个位置的 μ 值？

1.2 离散化：把积分变成加权和

连续积分在计算机里没法直接求解，必须离散化。把待重建的二维区域划分成 n×n 的像素网格，每个像素里的 μ 值就是未知数，把它们按行展开成一个长度为 N = n² 的列向量 x。对应地，每条射线沿路径的积分可以近似为路径穿过的每个像素的 μ 值与穿行长度的加权和。于是整条射线对应一个方程：

∑ w_ij · x_j ≈ p_i

把所有角度的所有射线方程堆在一起，就得到标准的线性方程组：

Ax = b

这里 A 是投影矩阵，行数是“角度数 × 探测器单元数”，列数是像素数 N；b 是实测投影向量。

从纯数学模型角度看，图像重建问题到这里已经变成了一个普通的代数方程求解问题。A 的规模通常大得吓人，比如一个 128×128 的图像配上 180 个角度、185 个探测器单元，A 就有 33300 行、16384 列。如果行数小于列数，方程组欠定；反过来则是超定。实际中为了稳定性和抗噪声能力，通常会布置尽量多的测量数据，使系统接近“满定”或“超定”，然后通过最小二乘或迭代方法求近似解。

1.3 差分方程模型在这里扮演什么角色

标题里的“代数方程与差分方程模型”不是生硬拼凑。代数方程模型显然对应 Ax=b；差分方程模型则体现在迭代重建的过程中。无论是 ART 还是 SIRT，迭代重建都不是一步到位的求解，而是构造一个离散递推序列：

x_{k+1} = x_k + Δx_k

每一次修正量 Δx_k 都是基于当前残差 r_k = b - Ax_k 计算出来的。这个公式本质就是一个一阶差分方程。更广义地说，对连续 Radon 变换在角度方向和探测器方向上的采样，本身也是一种离散化，和数值解偏微分方程时构造差分格式的思路完全一致。

把图像重建同时理解为“代数的静态求解”和“差分的动态逼近”其实很有价值，后面讲迭代算法时你会发现，很多反直觉的现象都能从这层理解推出来。

2. 三类主流重建思路，原理才是根本

2.1 中心切片定理：所有线性重建方法的地基

如果只记一个 CT 理论，选中心切片定理。它的表述非常漂亮：物体在角度 θ 下的一维投影做傅里叶变换，结果等于物体二维傅里叶变换在过原点、方向角 θ 的直线上的切片。

翻译成人话就是：每一个角度的投影频率信息，都只是二维频谱里的一根线。想要完整还原物体，理论上需要无穷多个角度的投影，把频域平面铺满。实际只能采集有限角度，所以频域某些区域信息不足，这就为伪影埋下了伏笔。

这个定理同时解释了为什么重建问题本质上是困难问题：因为我们永远只能拿到有损采样后的频域信息。

2.2 直接反投影为什么图像是糊的

最直觉的重建想法是反投影：把某个角度的投影值沿着入射方向均匀地“抹回”图像平面，把所有角度叠加起来。

这个过程数学上相当于对每个像素把所有穿过它的射线投影值求和。听起来合理，结果却不忍直视：重建出来的图像被严重模糊，物体边缘会有放射状星形伪影，均匀区域出现杯状畸变。原因是反投影本身并不是 Radon 变换的逆运算，它其实是在频域里给频谱乘了一个 1/|ω| 的权重因子，也就是说高频信息被系统性削弱了。

用一个生活类比想这件事：反投影像把几张“月光投影”叠在一起，每张投影都只是模糊的大致轮廓，叠完了还是模糊，而且多叠几次还会产生杂乱的星芒。

2.3 滤波反投影：补偿被削掉的高频

既然反投影在频域乘了 1/|ω|，那重建前先把投影在频域乘回 |ω| 不就行了？这就是滤波反投影 FBP 的核心思路。

具体分三步：

1. 对每个角度的一维投影做傅里叶变换；
2. 在频域乘一个斜坡滤波器，也就是幅值响应为 |ω| 的滤波器；
3. 逆傅里叶变换回空间域，得到滤波后的一维投影，再做反投影。

斜坡滤波器是理论上的精确选择，但直接把 |ω| 用起来会把高频噪声放大到离谱。所以工程上几乎总是要加窗函数做截断，最常用的是 Hamming 窗和 Hann 窗。加了窗之后，高频部分被压下来，图像更光滑，代价是分辨率略有下降。FBP 的优势是快、稳定，只要有 sinogram 就能在线性时间内重建，所以传统 CT 设备里长期都是主力算法。

2.4 迭代重建：绕开频率域直接解方程

FBP 建立在频域分析之上，而迭代重建选择直接面对代数方程组 Ax=b，不搞傅里叶变换。

最经典的代数重建技术 ART 是 Kaczmarz 方法的直接应用。思路是每次只拿一条射线的方程做修正：

x ← x + λ · (p_i - aᵢᵀx) / ‖aᵢ‖² · aᵢ

其中 aᵢ 是 A 的第 i 行，λ 是松弛因子。翻译成直觉语言：看当前解对这条射线的预测值，和实际测量差多少，然后把误差按照该射线穿过像素的权重大小“分摊”回像素上去。

SIRT 则是 ART 的温和版本，每轮迭代把所有射线的修正量累积起来，取平均后再统一更新一次。因为不逐线突变，SIRT 的抗噪声能力更强，收敛也更平稳，但需要更多轮迭代。

这类方法最大的优势在于：即使投影角度很少、测量不完全，迭代法仍然能给出比 FBP 更合理的结果，因为频率域里缺失的信息在迭代过程中会被空间域约束补回来一些。这也是低剂量 CT 和稀疏角度 CT 里迭代方法大受欢迎的根本原因。

3. 用 Python 实现一个完整的最小 CT 重建闭环

纸上谈兵没意思，直接上代码。我用 numpy 和 scipy 造一个极简但完整的模拟实验：生成模拟体模、构造投影矩阵、计算正弦图、跑 FBP 和迭代重建、最后对比质量。整个流程大概一百来行。

3.1 准备模拟数据：叠加椭圆的体模

真正的 CT 验证常用某个公开的合成体模，里面用一堆椭圆模拟不同组织的密度。这里我手动生成一个简化版，直接叠三个区域。它不涉及任何现实病人数据，纯粹是数学实验。

import numpy as np from scipy.sparse import csr_matrix from scipy.sparse.linalg import lsqr def make_phantom(size=128): yy, xx = np.mgrid[0:size, 0:size] cx, cy = size / 2.0, size / 2.0 img = np.zeros((size, size), dtype=np.float32) # 背景 img[(((xx - cx) / (size * 0.45)) ** 2 + ((yy - cy) / (size * 0.45)) ** 2) <= 1.0] = 0.2 # 一个大椭圆 img[(((xx - cx) / (size * 0.22)) ** 2 + ((yy - cy - size * 0.14) / (size * 0.16)) ** 2) <= 1.0] = 0.9 # 一个小圆 img[(xx - cx - size * 0.16) ** 2 + (yy - cy + size * 0.12) ** 2 <= (size * 0.04) ** 2] = 0.6 return img

生成出来是这样一块“灰色背景 + 亮色椭圆 + 亮色圆点”的合成物。内部灰度差异足够大，重建质量一眼就能看出来。

3.2 构造投影矩阵：最简单的权重模型

要跑迭代重建，光有正弦图还不够，我需要 A 矩阵。正规工程里会用射线逐像素求交的方式精确计算权重，典型代表是 Siddon 算法。演示代码里我用一个更粗暴的近似：把一个像素中心点投影到探测器的一维坐标上，落在某个探测器 bin 的宽度范围内，就认为这条射线对该像素有权重，权重直接取 1。

def build_weight_matrix(size, angles_deg, ndet, det_len): rows = len(angles_deg) * ndet cols = size * size cy = cx = size / 2.0 bins = np.linspace(-det_len / 2.0, det_len / 2.0, ndet) bin_w = (det_len) / ndet ri, ci, val = [], [], [] for ai, ang in enumerate(angles_deg): theta = np.deg2rad(ang) for di, bpos in enumerate(bins): r = ai * ndet + di for j in range(size): for i in range(size): x = i + 0.5 - cx y = j + 0.5 - cy s = x * np.cos(theta) + y * np.sin(theta) if abs(s - bpos) <= bin_w / 2.0: ri.append(r) ci.append(j * size + i) val.append(1.0) return csr_matrix((val, (ri, ci)), shape=(rows, cols)), bins

注意两点。第一，这个三重循环在 128×128 下能跑，但已经偏慢，实际工程中一定会向量化或用射线追踪。第二，权重值统一取 1 会损失定量精度，但定性演示完全够用，后面你可以看到伪影特征和理论预期完全一致。

3.3 正投影与带噪正弦图

有了 A 矩阵，正投影就是一次稀疏矩阵乘向量。

size = 128 angles = np.arange(0, 180, 2) # 90 个角度 ndet = 185 det_len = size * 1.5 phantom = make_phantom(size) A, bins = build_weight_matrix(size, angles, ndet, det_len) x_true = phantom.ravel() b_clean = A @ x_true rng = np.random.default_rng(42) noise = rng.normal(0, 0.03 * b_clean.std(), size=b_clean.shape) b = b_clean + noise

这里我加了 3% 标准差的高斯噪声。不加噪声也能重建，但加了才真实，因为真实探测数据一定是有噪声的，而且 FBP 对噪声的放大特性很典型。

3.4 手写滤波反投影

直接用某个图像库的现成函数当然可以，但既然我们要聊原理，最好自己写。滤波过程很简单：

def ramp_filter(p, det_pitch): n = p.size freq = np.fft.rfftfreq(n, d=det_pitch) pf = np.fft.rfft(p) * np.abs(freq) return np.fft.irfft(pf, n=n)

反投影部分，我直接用 A.T 来实现。理解这张矩阵：A 的行是“射线”，列是“像素”；A.T 就是把一条射线的值加到它经过的那些像素上，正是反投影操作。把每个角度的投影滤波后，乘到对应行上，再由 A.T 聚合回图像空间。

def fbp_reconstruct(A, b_sino, angles, ndet, det_len): ncol = A.shape[1] recon = np.zeros(ncol) det_pitch = det_len / ndet for r in range(len(angles) * ndet): row_vec = A.getrow(r).toarray().ravel() if row_vec.sum() == 0: continue di = r % ndet filtered = ramp_filter(b_sino[r - di : r - di + ndet], det_pitch) recon += filtered[di] * row_vec return recon.reshape(size, size)

严格来说，这里每个角度内部逐线滤波再反投影是可行的，只是效率不高。跑出来的结果会看到典型问题：边缘比直接反投影清楚得多，但噪声被抬高了，圆点边缘出现亮暗交替的振铃。

3.5 ART 和 SIRT 迭代求解

ART 的核心循环：

def art_reconstruct(A, b, size, n_epoch=20, lam=0.1): x = np.zeros(A.shape[1]) m = A.shape[0] for _ in range(n_epoch): for i in range(m): ai = A.getrow(i).toarray().ravel() denom = ai @ ai if denom < 1e-12: continue resid = b[i] - ai @ x x = x + lam * resid / denom * ai return x.reshape(size, size)

逐线更新的好处是每一轮都在快速修正，前几轮误差下降极快。缺点是噪声大的时候容易过拟合，λ 调太大还会发散。我试过 0.3 以上时结果剧烈震荡，0.05 到 0.1 比较稳。

SIRT 则换成整轮更新：

def sirt_reconstruct(A, b, size, n_iter=200, lam=0.5): x = np.zeros(A.shape[1]) row_sum = np.asarray(A.sum(axis=1)).ravel() col_sum = np.asarray(A.sum(axis=0)).ravel() row_sum[row_sum < 1e-12] = 1.0 col_sum[col_sum < 1e-12] = 1.0 for _ in range(n_iter): r = b - A @ x update = A.T @ (r / row_sum) x = x + lam * update / col_sum return x.reshape(size, size)

SIRT 每轮的计算量比 ART 大，但每一轮都不是“过冲式”修正，而是把各路信息平均后再走一小步，所以曲线更平滑。在噪声比较大的情况下，SIRT 重建出来的图像往往最干净。

3.6 结果对比：不要只用眼睛看

肉眼判断太主观，我再算两个常用指标：相对误差和结构相似度 SSIM。用一个小函数：

def rel_err(pred, true): diff = pred - true return np.sqrt(np.sum(diff ** 2)) / np.sqrt(np.sum(true ** 2))

在我的模拟配置下，典型结果是这样：

实际跑的时候你会注意到，直接反投影虽然难看，但它其实是 FBP 的“未滤波版本”，两者的计算速度都极快。迭代算法在 CPU 上跑几百轮要好几秒到十几秒，换回来越强的抗噪声能力和稀疏角度适应性。这里没有绝对最优，只有取舍。

4. 从模拟到工程，几个必须注意的坑

4.1 角度采样不够：一切伪影的根源

模拟中我把角度步长从 2 度改到 6 度，FBP 图像立刻出现大面积的放射状条纹，看起来像从几个方向向四周散开的光线。这是因为频域里的采样不够密，高频信息沿某些方向完全缺失。实际项目里，角度数和探测器单元数是直接决定重建上限的硬约束，算法再好也补不出信息都没采集进来的频率内容。

所以拿到一个重建任务，先问一句：投影数据质量如何、角度覆盖够不够。忽略这个前提去调算法纯属缘木求鱼。

4.2 A 矩阵的内存爆炸与稀疏化

如果不注意数据结构，直接用稠密矩阵存储 A，128×128 图像、33300 条射线时矩阵尺寸是 33300×16384，每个 float64 占 8 字节，总内存约 4.4GB。这还只是一个小规模的演示。

但理论上 A 里绝大多数元素都是 0，因为一条射线最多穿过几十个像素，所以正确做法是压缩稀疏行格式。上面代码里用的稀疏矩阵就是这个原因。工程上，投影矩阵的构建还会继续优化：用射线追踪遍历、只在像素边界处求交、只生成有效射线对应的行，能再省一个数量级的时间和内存。

4.3 迭代发散：松弛因子不是越大越快

我最初跑 ART 时，为了“快速收敛”，把 λ 直接设成 0.5，结果前几轮误差确实降得很快，但后面突然开始震荡，最终图像出现黑白相间的棋盘状噪声。这是典型的发散现象。

松弛因子的作用相当于差分方程递推里的步长。步长太大，残差修正会跨过真解；步长太小，收敛慢到让人失去耐心。经验值上，ART 的 λ 在 0.05～0.3 之间比较安全，SIRT 的 λ 可以到 0.5～1.0，因为它是平均更新，天然更稳。判断是否发散，不要只看误差绝对值，要看每轮误差曲线的末端趋势是下降还是反弹。

4.4 滤波器选择：噪声放大与分辨率损失的博弈

FBP 里斜坡滤波器的窗函数选择是个经典权衡。用纯斜坡滤波，图像最锐利，但投影里有噪声时会出现密密麻麻的高频毛刺。用 Hann 窗，噪声被压住了，但小细节也跟着糊了。

工程上的做法通常是：先估算投影数据的噪声水平，噪声高就选更保守的窗，噪声低可以激进一些。还有一个技巧是同时显示滤波前后的一维投影曲线，如果投影曲线在物体边缘处看着已经有明显抖动，那 FBP 出来基本没法看，这时候应该转身考虑迭代法。

4.5 常见问题速查表

我自己的习惯是：小规模实验先跑全流程，确认正投影和反投影配对正确；正式处理带噪声数据时，优先上 SIRT，用 200 轮左右的迭代，对比 FBP 的结果再下结论。经典方程模型的价值在于它让整条链路保持透明，哪里效果好、哪里效果差，全都可以追溯到矩阵和频率的细节上。做完这一次完整模拟，我对“CT 重建是一个病态反问题”这句话算有了实打实的理解。

这一期的“代数方程与差分方程模型（一）”先聊到这儿，CT 图像重建里还有动态序列和时空递推的场景，那就是差分方程模型进一步延伸的方向了。