张量层析成像与傅里叶方法：从标量CT到方向场重建

假设你是一个做图像或反问题方向的技术工程师平时对 CT 重建、偏微分方程反演、医学影像都多少有些接触。那你大概率已经习惯了这样一种设定在三维空间里我们想重建的目标是一个“每个位置只有一个数值”的标量场例如 X 射线衰减系数、声速扰动、温度场。很多教科书和开源库也都是按这个设定来组织算法的滤波反投影也好、代数重建也罢第一步都是把测量值投影到像素或体素网格上。但真实世界并不总是这么友好。复合材料在受力之后内部的应力状态不仅有大小还有方向地震波穿过地壳时传播速度会随偏振方向和传播方向变化光通过透明晶体或纤维材料时延迟和相位改变也依赖入射方向。这些场景中需要重建的对象不再是某个点上“一个数”而是“一个带着方向约束的物理量”。用数学语言说目标从标量场变成了张量场。这种问题在成像领域有一个专门的名字张量层析成像Tensor Tomography。在近期可看到的学术会议信息中国际基础科学大会ICBS的系列学术报告里有一场报告标题为“A Fourier method in Tensor Tomography with Applications”分享人是 David Omogbhe。只看题目就知道讨论的是如何用傅里叶方法处理张量层析成像并给出实际应用。这件事之所以值得单独写一篇技术解读是因为“张量 CT 傅里叶”三者组合之后很多经典直觉会失效。如果你只是把标量 CT 的软件包套到一个二阶张量场上大概率会得到错误甚至无意义的结果。这篇文章会把问题拆开讲清楚张量层析成像到底要解决什么、为什么傅里叶方法在这里能起作用、经典中心切片定理如何才能被推广到张量场合以及如果你想动手实验代码层面应该从哪里开始。文末还会给出一批常见误区和工程建议方便你把这个方向接入自己的项目。1. 这篇文章真正要解决的问题1.1 为什么标量 CT 不够用先回忆经典的二维 CT 模型。一束 X 射线穿过物体物体对射线的衰减可以表示为一个标量函数 (f(x,y))。当射线沿某个方向穿过物体时探测器得到的是这条路径上衰减系数的线积分。如果收集足够多角度下的线积分就可以通过中心切片定理、滤波反投影或者迭代重建得到 (f(x,y))。这个模型有一个隐含假设物体对射线的响应在各个方向上是“可累加”的而且衰减系数本身不依赖探测方向。在很多问题中这个假设是成立的。但应力、应变、扩散张量、各向异性波速这类物理量天然带有方向信息。举个例子如果你想知道一块材料内部每一点的应力主轴方向你需要的不是一个灰度图而是每个点的 2×2 应力矩阵或 3×3 应力张量。这种每个位置都有一个矩阵或更高阶数组的场就叫张量场。1.2 张量层析成像要做的“升级”张量层析成像的研究目标是从一系列沿直线的积分测量中重建一个张量场而不是标量场。听起来像是一个简单的“多通道 CT”既然每个张量场有多个分量那就把每个分量当成一个独立的标量场分别做 CT 不就行了问题恰恰出在这里。实际测量值通常是张量某个方向和某个张量分量线性组合后的积分而不是张量分量本身。比如在偏振光应力测量中一束光沿方向 (v) 穿过样品时测量到的相位延迟可能正比于该点应力张量 (\boldsymbol\sigma) 在方向 (v) 上的二次型 (\sum_{i,j}\sigma_{ij}v_i v_j)。你测量到的不是某个分量的直接投影而是全部分量被方向 (v) 加权后的总和。因此单条射线的测量值是多个未知分量共同贡献的结果。这种“混叠”给重建带来了本质困难也是张量层析成像区别于普通 CT 的关键。1.3 谁最需要读这篇文章如果你属于以下任何一类读者这篇文章对你会有帮助你平时做 CT、MRI、超声或地震成像现在遇到了“各向异性”“方向相关”“张量约束”等新需求你阅读论文时看到 Fourier method、Tensor Tomography 这类词希望快速知道它们之间的关系你需要把张量场重建从理论公式变成可运行的 Python 代码你正在参加相关领域的数学报告想提前补足背景知识。这篇文章不会把全部篇幅放在纯数学推导上而是先用直观模型建立物理图像然后给出能落地的数值实验思路。2. 张量层析成像的核心概念与基础模型2.1 什么是张量场张量可以理解为对标量、向量、矩阵的推广。零阶张量是标量一阶张量是向量二阶张量通常用矩阵表示。在数学物理中应力、应变、介电常数、扩散系数等都是二阶张量。张量场就是空间中每个位置都定义了一个张量的函数。例如二维平面中的一个二阶对称张量场[ \mathbf f(x,y) \begin{bmatrix} f_{11}(x,y) f_{12}(x,y)\ f_{12}(x,y) f_{22}(x,y) \end{bmatrix} ]因为这里假设对称所以独立分量是三个标量场(f_{11})、(f_{12})、(f_{22})。更高阶情况类似只是分量数量更多。需要注意的是张量场中的每个分量在数学上往往是相互关联的不是因为图像处理上的“多通道”而是因为物理定律决定了它们之间要满足某些微分约束。比如在无旋条件下一个弹性应变场可能是某个位移场的对称梯度。这种内在约束是后续消歧问题的根源。2.2 张量场的方向敏感性为什么普通 CT 不能直接用来重建张量场因为张量场在不同方向上的“表现”不同。想象你在二维平面上观察一个应力张量场。选择一条单位方向向量 (v)张量场在 (v) 方向的“表现”通常是[ q(x,v)\sum_{i,j1}^{2} f_{ij}(x),v_i,v_j ]如果 (f) 是应力张量这个量就是应力在方向 (v) 上的法向分量如果 (f) 是扩散张量这个量就与方向 (v) 上的扩散系数有关。所以当你沿着某条射线做线积分时测量到的不是 (f_{11}) 独立积分也不是 (f_{22}) 独立积分而是它们被方向权重 (v_i v_j) 混合后的积分。为了重建完整张量场往往需要从多个方向获取数据并利用张量本身的结构把混叠解开。2.3 张量场射线变换的定义射线变换是一种把“函数沿直线积分”的数学算子。对于标量场 (f)在二维平面中沿方向 (v) 做线积分可以写成[ P f(x_0, v)\int_{-\infty}^{\infty} f(x_0 t v),dt ]但对于张量场这个定义要改成带权重的线积分。以一阶张量场向量场为例可以定义[ I_1 \mathbf u(x_0,v)\int_{-\infty}^{\infty} \mathbf u(x_0t v)\cdot v,dt ]也就是把向量 (\mathbf u) 先投影到方向 (v) 上再沿直线积分。对 (m) 阶张量场 (\mathbf f)可以定义更一般的权重[ I_m\mathbf f(x_0,v)\int_{-\infty}^{\infty} \mathbf f(x_0t v)\cdot v^{\otimes m},dt ]其中 (\cdot) 表示张量全收缩(v^{\otimes m}) 表示把 (v) 自己乘 (m) 次得到的高阶张量。对二阶对称张量 (\mathbf f(f_{ij}))上式具体写出来就是[ I_2\mathbf f(x_0,v) \int_{-\infty}^{\infty}\left( f_{11},v_1^2 2f_{12},v_1v_2 f_{22},v_2^2 \right)dt ]这里的分量都是在路径点 (x_0t v) 处取值的。为什么要定义得这么“绕”因为它直接对应物理测量。在很多探测手段中一个沿着方向 (v) 传播的波或粒子对张量场的响应天然就是这种二次型加权后的线积分。如果把加权去掉反而无法描述真实传感器。2.4 标量 CT 与张量 CT 的对比对比维度标量 CT张量层析成像待重建对象标量场 (f(x))张量场 (\mathbf f(x))每个位置未知数1 个多个分量典型线积分(\int fdt)(\int \mathbf f\cdot v^{\otimes m}dt)测量值与方向关系只依赖路径不依赖张量方向测量值随方向 (v) 明显变化经典重建方法FBP、代数迭代、傅里叶切片傅里叶域分解、消歧、反投影修正主要难点噪声、有限角、伪影分量混叠、核空间、结构消歧这张表最值得记住的是最后一行张量层析成像的主要难点不只是“噪声更多”或“计算量更大”而是存在一类数据无法分辨的分量。这是标量 CT 中没有的新问题。3. 从标量 CT 到张量 CT傅里叶方法的切入点3.1 重新理解中心切片定理在经典 CT 中傅里叶方法的核心是中心切片定理。假设一个二维标量场 (f(x,y))沿方向 (v) 做一系列平行射线投影并把投影结果关于射线位置 (t) 做一维傅里叶变换。这个定理说一维