CT 技术已经相当成熟的时候为什么还会有人把“层析成像”和“张量场”放在一起做研究今年国际基础科学大会ICBS上David Omogbhe 有一场报告标题是 “A Fourier method in Tensor Tomography with Applications”。很多做工程或者算法落地的同学看到这个题目第一反应往往是这和我熟悉的 X 射线 CT 是一回事吗如果不是它又多难Fourier 方法在这里到底起了什么作用这篇文章想帮你把这三个问题一次讲清楚。我的判断很明确张量层析成像不是“把一张图变成三张图再各自做 CT”那么简单它的难点在于物理量本身带有方向性而 Fourier 方法之所以重要是因为它能把复杂的耦合测量在频域里解开一部分让“哪些信息能恢复、哪些信息注定恢复不了”变得清晰。读完这篇文章你可以理解张量层析成像的基本测量模型知道 Fourier 方法在其中解决什么问题并且跟着一个最小 Python 实验直观体会频域重建流程是如何跑通的。1. 这篇报告到底在讲一个什么问题先看标题的层次。“A Fourier method in Tensor Tomography with Applications”里有三个关键词Fourier method、Tensor Tomography、Applications。Tensor Tomography 通常翻译为张量层析成像有些人也叫它张量断层重建。传统 CT 里的待重建对象是密度这类标量场。也就是说空间中每个位置只需要一个数来描述它。你拍一张投影得到的是一个数沿着射线方向的积分收集多个角度之后再用各种重建算法还原出空间分布。张量层析成像要处理的问题则是空间中每个位置不是一个数而是一个张量场。最常见的场景是二阶对称张量场可以理解为一个 2×2 或者 3×3 的矩阵矩阵里有多个分量。你沿一条射线测到的数据往往是“张量场在这条射线方向下的某种组合”而不只是某一分量的简单积分。这就直接带来了两个新问题。第一个问题是测量模型变复杂了。不同的入射方向会激发不同的分量组合。如果设计实验时不充分考虑方向覆盖很容易出现一种情况手头数据看似很多但某些待求分量根本没有被激励到也就是所谓的不可观测信息。第二个问题是重建算法的适用性问题。普通 CT 中有 Fourier 切片定理做基础把空间域里的射线积分转换到频域会得到一个非常漂亮的切片关系。到了张量情形每一个张量分量都可能有对应的投影关系同时分量之间又存在耦合算法设计不能直接套用标量重建。从公开数学脉络来看这类报告通常不会只停留在“换个名词”的层面而是会回答给定什么样的测量数据什么时候解唯一在不适定情况下哪些部分可以稳定重建哪些部分需要正则化以及 Fourier 方法如何把反演化成频域乘子或者更容易处理的积分算子。因此这篇报告对三类人特别有价值。第一类是做成像算法的研究者尤其是研究医学成像、无损检测、地震成像中向量或张量场恢复的人。第二类是做数值计算和反问题的人他们想了解 Fourier 工具在积分几何里如何自然出现。第三类是工程同学如果你的项目涉及应力光学测量、弹性波、各向异性介质理解这套数学语言的边界比记住几个公式更重要。2. Tensor Tomography 的核心概念与测量模型要给这篇报告建立直觉先要知道最经典的射线积分。对于一个标量场 (f(x))一条射线 (\gamma(t)) 上的 CT 测量可以写成[ P f(\theta,s) \int f(s\theta t\theta^{\perp}) dt ]从物理上看这就是平行束投影。不同的旋转角度 (\theta) 和偏移量 (s) 会组成一张正弦图 sinogram。重建目标就是从这张正弦图恢复 (f)。张量层析成像的不同之处在于被测对象 (F(x)) 是一个张量场。以二维二阶对称张量为例它由三个独立分量描述[ F(x) \begin{bmatrix} f_{11}(x) f_{12}(x) \ f_{12}(x) f_{22}(x) \end{bmatrix} ]如果只把它当成三个标量场分别用三次 Radon 变换重建就会忽略一个关键事实这三个分量在物理上并不是独立出现的它们常常满足某种偏微分约束。比如应力场要满足平衡方程应变场要满足协调性条件。如果强行拆开重建重建结果会违背这些约束产生无物理意义的伪影。张量层析的常见测量还多出一个自由度射线方向上的“方向性”。一种典型测量是对射线方向向量 (v) 做二次型投影[ I F(v,\xi) \int \langle v, F(\xitv)v \rangle dt ]这里 (v) 是射线前进方向数据不只来自某个分量而是 F 在 v 方向上的“有效权重”。这也解释了为什么张量层析不是标量 CT 的简单扩展标量 CT数据是一条数轴上的函数重建目标也是单值函数。张量 CT数据可能是多条数轴、多个极化方向的函数重建目标是一组分量。标量 CT 的 Fourier 切片关系是单个张量 CT 则要处理一组耦合关系。如果要用一个类比标量 CT 像你用一支笔扫描一页纸上的墨迹浓度黑白深浅是唯一信息张量 CT 则更像你推一块材料的不同方向观察它在各个方向的响应最后反推出材料内部“哪个方向强、哪个方向软、会不会回弹”这些更丰富的信息。理解这一点后就不会被“Tensor”这个词吓倒也不会低估它的难度张量场本身不是难点难点是数据里包含的方向混合。3. Fourier 方法在这一类问题里的角色为什么 Fourier 变换会在张量层析成像里反复出现这是整篇报告最关键的技术内核。在标量 CT 中Fourier 切片定理说的是对一个二维函数做某方向的一维投影其一维 Fourier 变换等于该函数二维 Fourier 变换在对应方向中心线上的取值。这意味着如果从 0 到 180 度都有投影你实际上是在频域里“画满”了原函数的频域支撑。这个定理的核心作用是解耦。一条射线投影在空间里像一个积分但是到频域之后它变成了一个“采样”。采样本身仍然是局部的每一个角度贡献频域的一条线但反演不再需要在空间域里逐个角度去解复杂的积分方程。张量层析中Fourier 方法的价值主要体现在三个层面。第一分析零空间和唯一性。当被测量对象从标量变成张量时会出现一些在空间域看不出原因的不可恢复信息。通过 Fourier 变换可以把这些不可恢复部分和对测量不敏感的部分映射到频域里的某些频率组合。如果某个方向的分量在频域里根本没被激励到无论优化算法多强都无法凭空恢复它。第二构造反演公式。对某些张量射线变换其法算子可以在频域里表示为乘子形式。于是反演问题从复杂的积分方程化简为对一个或多个频域乘子的处理。这也就是标题中 “A Fourier method” 的含义之一不是“用 FFT 加速一下”而是从算子谱结构出发设计反演策略。第三确定稳定恢复条件。实际采样不可能做到所有角度连续覆盖。有限角度、截断边界、噪声干扰都会让频域支撑产生空洞。Fourier 视角能让你判断哪些分量对角度覆盖更敏感哪些分量受噪声影响更大从而设计更合理的采集方案。不过需要注意Fourier 方法不是“银弹”。当张量场有内部约束时频域关系会变得更加复杂当数据离散且噪声明显时还要借助正则化。报告标题里的 Fourier method 通常解决的是“理想测量条件下如何理解和反演算子”工程上落地仍需要处理大量的数值细节。4. 环境准备与前置条件接下来做一个最小数值实验帮助你建立直观手感。这个实验不会完整复现论文里的张量反演算法但会覆盖三件事生成一个满足平衡条件的二维对称张量场而不是随便填一个矩阵。模拟射线方向的二次型投影得到多角度正弦图。用 Fourier 频域滤波加反投影的方式验证标量重建的最小流程。这样设计的目的是把一个大的数学问题拆开先理解张量场本身的约束再理解投影数据的生成最后理解频域重建中的“投影-滤波-反投影”套路。如果你能跑通这个最小实验再去看论文里的推导会轻松很多。实验环境并不复杂Python 3.8 以上。numpy 1.21 以上。scipy 1.7 以上用于插值采样。matplotlib 3.5 以上用于画图。对硬件没有要求普通笔记本即可。版本不需要刻意追新以上都是常见版本。如果安装依赖报错可以新建一个虚拟环境然后安装python -m venv tensor-tomo source tensor-tomo/bin/activate pip install numpy scipy matplotlibWindows 用户把source tensor-tomo/bin/activate换成tensor-tomo\Scripts\activate即可。5. 完整示例代码实现代码会一步一步展开。你可以把下面三段代码合到同一个脚本demo_tensor_tomography.py中也可以分段在 Jupyter Notebook 里运行。5.1 生成满足散度约束的张量场在二维无体力平衡应力场中一个常见的技巧是 Airy 应力函数。如果定义[ \sigma_{11} \partial_{yy} \phi, \quad \sigma_{22} \partial_{xx} \phi, \quad \sigma_{12} -\partial_{xy} \phi ]那么张量场自动满足平衡方程[ \partial_x \sigma_{11} \partial_y \sigma_{12} 0 ][ \partial_x \sigma_{12} \partial_y \sigma_{22} 0 ]这里选择 (\phi(x,y)\cos(\pi x)\cos(\pi y))得到一组既光滑、又满足约束的非平凡张量场。通过数值散度检查可以确认场生成的代码没有犯低级错误。# demo_tensor_tomography.py import numpy as np import matplotlib.pyplot as plt from scipy.ndimage import map_coordinates nx, ny 256, 256 xs np.linspace(-1.0, 1.0, nx) ys np.linspace(-1.0, 1.0, ny) X, Y np.meshgrid(xs, ys, indexingij) dx xs[1] - xs[0] dy ys[1] - ys[0] def d_dx(f): return np.gradient(f, axis0) / dx def d_dy(f): return np.gradient(f, axis1) / dy def stress_from_airy(X, Y): 由 Airy 应力函数 phicos(pi x)cos(pi y) 生成平衡张量场。 pi2 np.pi ** 2 S11 -pi2 * np.cos(np.pi * X) * np.cos(np.pi * Y) S22 -pi2 * np.cos(np.pi * X) * np.cos(np.pi * Y) # S12 -d_xy(phi)推导后得到 -pi^2 sin(pi x) sin(pi y) S12 -pi2 * np.sin(np.pi * X) * np.sin(np.pi * Y) return S11, S12, S22 S11, S12, S22 stress_from_airy(X, Y) # 验证张量场满足散度约束 div1 d_dx(S11) d_dy(S12) div2 d_dx(S12) d_dy(S22) print(div1 max abs , np.abs(div1).max()) print(div2 max abs , np.abs(div2).max()) fig, axes plt.subplots(1, 3, figsize(14, 4)) axes[0].imshow(S11.T, originlower, extent[-1, 1, -1, 1], cmapRdBu) axes[0].set_title(S11) axes[1].imshow(S12.T, originlower, extent[-1, 1, -1, 1], cmapRdBu) axes[1].set_title(S12) axes[2].imshow(S22.T, originlower, extent[-1, 1, -1, 1], cmapRdBu) axes[2].set_title(S22) plt.savefig(stress_components.png, dpi120, bbox_inchestight) print(saved stress_components.png)这段代码的关键在于用indexingij生成网格使得数组第 0 个轴对应物理 x 方向第 1 个轴对应物理 y 方向后面求散度时才不会方向错乱。运行后如果散度最大值在 (10^{-13}) 量级说明张量场满足约束可以用来生成测量数据。5.2 模拟射线方向的二次型投影有了张量场之后下一步是模拟一条射线。在二维平面中一条平行束射线可以这样描述单位方向向量为 (d)偏移向量为 (s) 沿单位法向量 (n)射线从侧边进入区域沿线进行积分。测量值取射线方向 (d) 对张量场的二次型[ g \int_{-R}^{R} d^T \sigma(x(t)) d , dt ]因为 S12 是对称分量展开后的被积函数是[ d_x^2 S11 2 d_x d_y S12 d_y^2 S22 ]采用map_coordinates做双线性插值获取射线上每个采样点处的分量值。这样避免了写复杂的判断同时保证代码能处理射线不完全落在网格点上的情况。def field_at(field, xs, ys, xq, yq): 在物理坐标 (xq, yq) 上对二维场做双线性插值。 row (xq - xs[0]) / (xs[1] - xs[0]) col (yq - ys[0]) / (ys[1] - ys[0]) return map_coordinates( field, [row, col], order1, modeconstant, cval0.0, ) def ray_second_form_integral(theta_deg, offset, R1.8, N256): 计算一条平行束射线上的 d^T sigma d 线积分。 theta_deg: 射线方向角单位度。 offset: 射线离原点的偏移量单位与网格坐标相同。 theta np.deg2rad(theta_deg) d np.array([np.cos(theta), np.sin(theta)]) n np.array([-np.sin(theta), np.cos(theta)]) center offset * n t_grid np.linspace(-R, R, N) xq center[0] t_grid * d[0] yq center[1] t_grid * d[1] v11 field_at(S11, xs, ys, xq, yq) v22 field_at(S22, xs, ys, xq, yq) v12 field_at(S12, xs, ys, xq, yq) measured d[0] ** 2 * v11 2.0 * d[0] * d[1] * v12 d[1] ** 2 * v22 return np.sum(measured) * (2.0 * R / (N - 1)) # 生成多角度正弦图 angles np.linspace(0.0, 179.0, 36) offsets np.linspace(-1.35, 1.35, 64) sinogram np.zeros((len(angles), len(offsets))) for ia
