张量层析中的Fourier方法：从Radon变换到频域重建

最近在翻国际基础科学大会相关学术报告的资料列表时我注意到 David Omogbhe 的一个报告题目“A Fourier method in Tensor Tomography with Applications”中文直译是《张量层析中的一种 Fourier 方法及其应用》。这个题目初看并不友好Tensor Tomography 是什么意思Fourier method 在层析成像里不是已经很成熟了吗为什么还要单独拿出来讲这篇文章不打算逐句复述报告内容而是从标题中拆出三条线索Fourier method、Tensor Tomography、Applications。我会把“张量层析为什么比普通 CT 更难”“Fourier 方法在这里解决了什么问题”“这类方法可以落到哪些应用场景”讲清楚并用 Python 做几个可复现的数值实验。适合数学成像、数值计算、图像重建和反问题方向的读者收藏零基础的人也能从 Radon 变换开始跟读。1. 背景从“只看灰度”到“要看方向”1.1 普通 CT 为什么只关心标量场医院里的 X 射线 CT可以简化为如下的物理模型射线从光源发出穿过人体组织到达探测器。人体不同组织对 X 射线的衰减能力不同这个衰减系数通常用一个标量函数表示[ \mu(x), \quad x \in \mathbb{R}^2 ]也就是说空间中任意一点只需要一个数值来描述它对射线的衰减强弱。探测器收集到的是射线路径上所有点的衰减值之和数学上是一个线积分。图像重建要解决的问题是我知道大量方向上的线积分如何反推出每个位置的 (\mu(x))。这是一个经典的积分几何反问题。大家熟悉的“滤波反投影FBP”“代数重建ART”“直接 Fourier 重建”都是针对这种标量场模型设计的。1.2 张量层析是什么但现实世界有很多物理量并不是“一个点一个数”而是“一个点一个方向相关的量”。举几个例子地壳中的应力状态需要二阶对称张量描述透明材料内部的残余应力同样是一个张量场流体中的速度场是向量场某些各向异性介质中的波传播速度还可能与方向有关。当目标物理量是向量场或张量场时投影过程不再只是沿着射线做简单累加。投影结果会同时取决于射线方向目标张量场在该点的取值张量场与射线方向的耦合方式。也就是说同样的位置换一个方向照射得到的线积分结果不仅会因为路径不同而改变还会因为“张量分量在这个方向上的投影权重”不同而改变。这类问题统称为张量层析Tensor Tomography。1.3 为什么值得专门研究如果目标只是“标量函数”Radon 变换和 Fourier 切片定理已经给出了一套非常优美的理论。但向量场和张量场比标量函数复杂得多。以二阶对称张量场为例在二维空间中每个位置需要三个独立分量来刻画在三维空间中则需要六个独立分量[ A(x) \begin{bmatrix} A_{11}(x) A_{12}(x) \ A_{12}(x) A_{22}(x) \end{bmatrix} ]这就像是把普通 CT 中“每个像素一个未知数”的问题升级成“每个像素三个未知数”甚至“每个像素多个耦合未知数”的问题。可观测的投影数量却不一定成倍增加。因此张量层析的核心问题从一开始就是什么样的测量数据能够恢复张量场能够恢复哪一部分哪些部分是本质上无法恢复的Fourier 方法在这类问题中的价值不只是提供一种快速算法更重要的是在频域中把“可恢复信息”和“不可恢复信息”区分开。2. 数学地基Radon 变换与 Fourier 切片定理2.1 从 Radon 变换说起要理解张量层析中的 Fourier 方法最好先把标量 CT 的数学框架写清楚。二维空间中任何一个方向可以用单位向量表示[ \theta (\cos\varphi, \sin\varphi) ]垂直于 (\theta) 的方向记为[ \theta^{\perp} (-\sin\varphi, \cos\varphi) ]一条直线可以参数化为[ x(t) s\theta t\theta^{\perp} ]这里 (\varphi) 是直线的法向角度(s) 是原点到直线的带符号距离(t) 是沿直线的位置。标量函数 (f(x)) 的 Radon 变换定义为[ Rf(\varphi, s) \int_{-\infty}^{\infty} f(s\theta t\theta^{\perp}) dt ]也就是沿着这条直线的线积分。CT 采集到的正弦图sinogram本质上就是不同角度、不同偏移下的 (Rf(\varphi, s))。2.2 Fourier 切片定理为什么是核心Radon 变换的逆问题看起来非常复杂每个投影都是一条路径上的积分要把所有路径上的积分整合成二维函数。Fourier 切片定理给出了一个漂亮的等价描述对投影 (Rf(\varphi, s)) 关于 (s) 做一维 Fourier 变换等价于对 (f(x)) 做二维 Fourier 变换后在过原点的、方向为 (\theta) 的直线上取值。用公式表达为[ \mathcal{F}_1 \left Rf(\varphi, \cdot) \right\mathcal{F}_2 f(\sigma\theta) ]其中 (\mathcal{F}_1) 是一维 Fourier 变换(\mathcal{F}_2) 是二维 Fourier 变换(\sigma) 是空间频率。这个定理的意义非常重大Radon 变换把“实空间中的线积分”变成了“频域中的一条直线采样”。一条投影数据经过一维 FFT 后可以直接填到二维频谱中的一条直线上。如果采集了足够多角度的投影就能用这些一维谱“拼”出一张完整的二维频谱图。最后做一次二维逆 FFT就得到了重建图像。2.3 Fourier 方法给后续研究留下了什么从工程角度看Fourier 方法可以实现“快捷重建”但也会遇到插值、采样不均匀、角度稀疏等实际问题因此后来又发展出了非均匀快速傅里叶变换NUFFT、迭代重建方法等。从数学角度看Fourier 切片定理真正重要的地方在于它把一个积分几何问题转化成了频域抽样问题。频域里哪些位置有数据、哪些位置没有数据哪些频率分量可识别、哪些频率分量不可识别都变得一目了然。张量层析之所以强调“Fourier method”也正是因为这种方法论可以延续到张量场先推导张量投影在频域中的对应关系再分析可恢复的子空间最后做重建。3. 张量层析的问题为什么比标量 CT 难3.1 张量场的纵向投影先固定一个相对简单的场景二维空间中的二阶对称张量场。假设有一条方向为 (\theta) 的射线射线方程为[ \gamma(t) s\theta t\theta^{\perp} ]张量场与射线方向水平投影耦合之后可以得到一种最基本的张量 Radon 变换通常称为纵向变换longitudinal transform[ g(\varphi, s)\int_{-\infty}^{\infty} \theta^T A(\gamma(t)) \theta , dt ]把二阶对称张量展开[ \theta^T A \theta\cos^2\varphi \cdot A_{11} 2\sin\varphi\cos\varphi \cdot A_{12} \sin^2\varphi \cdot A_{22} ]所以一条纵向投影数据并不是某个独立分量 (A_{11}) 的投影而是三个分量按方向权重组合后的线积分。这是张量层析与标量 CT 的第一个重要差异标量 CT 中每次投影直接对应同一个标量函数张量 CT 中每次投影对应的是三个分量的线性组合。3.2 仅靠纵向投影为什么不够一个很自然的问题是如果我把所有方向的纵向投影都采集齐了能不能恢复出 (A_{11}, A_{12}, A_{22}) 三个分量答案通常是不能完全恢复。张量场中有一类特殊结构被称为“势张量场”potential tensor field。它由某个向量场 (v) 的对称梯度构成[ (\nabla^s v)_{ij}\frac{1}{2} \left( \frac{\partial v_i}{\partial x_j} \frac{\partial v_j}{\partial x_i} \right) ]在合适的边界条件下这类张量场沿任何直线的纵向投影结果可能恰好为零。也就是说它们位于纵向 Radon 变换的零空间中。这意味着如果你只测量纵向投影无论测量多少个角度都无法区分一个张量场和叠加了某个“势张量场”之后的张量场。这个结论最早由 Sharafutdinov 在张量积分几何理论中系统研究过。很多张量层析论文讨论“数据完备性”时都会回到这个问题上。3.3 我们能恢复什么既然完整张量场不可恢复那实际工作中能恢复的是哪一部分答案是张量场分解之后的可测部分。二阶张量场可以进行 Helmholtz 型分解。一个常用的做法是把张量场分解为“无源部分”solenoidal part和“势部分”potential part。纵向 Radon 变换丢失的正是势部分也就是对称梯度部分。因此通过纵向投影可以恢复无源部分。这就带来两个工程方向如果应用场景只需要无源分量那么纵向投影数据就够用如果应用场景要求完整张量场则需要引入额外的测量信息例如横向投影、混合型投影或者已知的物理边界条件。在后一种情况下Fourier 方法能起到至关重要的作用它可以把“缺什么数据”“哪些频率缺失”这类信息放在频域中显式表达出来。4. Fourier 方法在张量层析中的核心思路4.1 频域里的切片关系现在把 Fourier 切片定理推广到张量场纵向变换。设 (\theta (\cos\varphi, \sin\varphi))纵向投影为[ g(\varphi, s)\int_{-\infty}^{\infty} \left[ \cos^2\varphi A_{11} 2\sin\varphi\cos\varphi A_{12} \sin^2\varphi A_{22} \right] (s\theta t\theta^{\perp}) , dt ]对 (s) 做一维 Fourier 变换利用标量场的 Fourier 切片定理每个分量都会得到对应的频域切片。因此[ \mathcal{F}_1 g(\varphi, \cdot)\cos^2\varphi \cdot \widehat{A_{11}}(\sigma\theta) 2\sin\varphi\cos\varphi \cdot \widehat{A_{12}}(\sigma\theta) \sin^2\varphi \cdot \widehat{A_{22}}(\sigma\theta) ]这个公式看起来只是把标量切片定理逐项套用但它提供的信息量非常大。它告诉我们在频域中沿任何一条过原点的直线张量场三个分量的 Fourier 变换并不是独立采样而是只能采到它们的某种加权组合。如果只有纵向投影在每个频域采样中心处三个未知分量共享一个方程。这个欠定关系就解释了为什么完整张量场无法从纵向投影中恢复。4.2 如何设计 Fourier 重建方法如果只有以上一个方程问题是欠定的。要让问题可解通常需要增加约束。常见思路包括采集横向投影或混合型投影。不同方向的二次型耦合方式不同在频域中可以提供更多独立的线性方程利用张量场的物理约束。例如在应力场中张量分量之间可能满足平衡方程把未知场的结构限制在无源子空间内。这在很多应用中本身就是物理上合理的假设在频域中用最小二乘或正则化方法求解。David Omogbhe 报告标题中的 Fourier method正是这一类把测量数据变换到频域、建立代数方程、再反演目标张量场的方法。它相比直接在实空间推导积分方程更直观也更方便分析高频信息和稳定性。4.3 有限角度与频域空洞问题实际工程中测量角度往往有限。比如工业 CT 中旋转机构有角度限制地震层析中射线角度由震源位置决定无法自由布设。在标量 CT 中有限角度会导致频域中形成“空洞”重建图像会出现方向性伪影。在张量 CT 中有限角度的影响更加复杂某些方向上的高频信息完全缺失即使频域覆盖相对完整欠定部分仍然存在噪声对高频分量的放大效应可能比标量 CT 更严重若只采集纵向投影还需要额外处理势分量的不确定性。因此张量层析的 Fourier 方法不只是“快速重建”工具更多时候它是分析的起点先判断哪些频率能识别再决定用什么样的正则化策略。5. Python 数值实验动手验证核心思想前面的理论比较抽象下面用 Python 做三个实验离散验证标的 Fourier 切片定理生成一个二维二阶对称张量场计算该张量场的纵向投影正弦图。实验均在 Jupyter Notebook 或普通 Python 脚本中完成。核心依赖只有 numpy、scipy、matplotlib。5.1 环境准备版本不需要完全一致以下组合是比较常见的Python 3.8 及以上numpy 1.21 及以上scipy 1.7 及以上matplotlib 3.5 及以上。代码原理是通用的不同版本通常都能运行。5.2 实验 1验证 Fourier 切片定理先构造一个简单的二维测试图像模拟“一个圆形目标内含两个小结构”的灰度场。import numpy as np import matplotlib.pyplot as plt from scipy.ndimage import rotate def make_test_image(n256): 构造一个二维测试场边界为 0内部有圆盘和小目标。 x np.linspace(-1, 1, n) xx, yy np.meshgrid(x, x) f np.zeros((n, n)) # 大圆 f[xx**2 yy**2 0.55**2] 0.9 # 内部的小圆 f[(xx - 0.18)**2 (yy 0.08)**2 0.13**2] 0.4 # 内部的小方块 f[(xx 0.22)**2 (yy - 0.10)**2 0.10**2] -0.6 return f def approximate_radon(f, angles_deg): 近似 Radon 变换 将图像旋转 angle 度再沿某个轴求和。 这只是演示用近似实现不是精确解析投影。 projections [] for angle in angles_deg: rotated rotate(f, angle, reshapeFalse, order1, cval0.0) projections.append(rotated.sum(axis0)) return np.asarray(projections)这里我没有调用 skimage 的 radon 函数而是用旋转求和近似。原因有两个减少第三方依赖能直观看到“旋转图像后按行求和”和“沿某个方向做线积分”之间的关系。需要注意这种近似的数值精度不如精确射线追踪因此后续的频谱比较以趋势判断为主。接下来验证中心切片定理。def ifftshift_fft_1d(signal): 把信号中心移到索引 0再做一维 FFT最后把零频搬回中央。 return np.fft.fftshift(np.fft.fft(np.fft.ifftshift(signal))) def ifftshift_fft_2d(image): 二维版本的中心化 FFT。 return np.fft.fftshift(np.fft.fft2(np.fft.ifftshift(image))) # 生成测试图像 f make_test_image(256) n f.shape[0] # 二维频谱 F2 ifftshift_fft_2d(f) # 0 度方向投影的一维频谱 p0 approximate_radon(f, [0])[0] P0 ifftshift_fft_1d(p0) # 90 度方向投影的一维频谱 p90 approximate_radon(f, [90])[0] P90 ifftshift_fft_1d(p90) # 二维频谱的水平和垂直中心切片 horizontal_slice F2[n // 2, :] vertical_slice F2[:, n // 2] # 计算相关系数 corr0 np.corrcoef(np.abs(P0), np.abs(horizontal_slice))[0, 1] corr90 np.corrcoef(np.abs(P90), np.abs(vertical_slice))[0, 1] print(f0 度方向投影谱与二维频谱水平切片相关系数: {corr0:.4f}) print(f90 度方向投影谱与二维频谱垂直切片相关系数: {corr90:.4f}) freq_axis np.fft.fftshift(np.fft.fftfreq(n)) fig, axes plt.subplots(1, 2, figsize(12, 4)) axes[0].plot(freq_axis, np.abs(P0), label1D FFT of 0 degree projection) axes[0].plot(freq_axis, np.abs(horizontal_slice), --, label2D FFT horizontal slice) axes[0].set_title(0 degree slice) axes[0].legend() axes[1].plot(freq_axis, np.abs(P90), label1D FFT of 90 degree projection) axes[1].plot(freq_axis, np.abs(vertical_slice), --, label2D FFT vertical slice) axes[1].set_title(90 degree slice) axes[1].legend() plt.tight_layout() plt.show()运行后可以看到两条曲线基本重合。相关系数通常