从‘弹道’到‘散射’：手把手教你用Python模拟光子在不同散射介质中的传输路径

从‘弹道’到‘散射’：手把手教你用Python模拟光子在不同散射介质中的传输路径

当一束光穿过牛奶或生物组织这类浑浊介质时，光子们会经历怎样的冒险旅程？有的光子像子弹般直线穿出，有的像蛇一样蜿蜒前行，还有的则像醉汉般跌跌撞撞。理解这些不同行为的光子，对于医学成像、大气遥感和材料检测等领域至关重要。本文将带你用Python代码亲手构建一个光子传输模拟器，直观展示三种典型光子的运动轨迹。

1. 光子传输的物理基础与建模思路

在浑浊介质中，光子的命运由**散射平均自由程（SMFP）**决定——这是光子两次散射之间平均走过的距离。根据穿透深度与SMFP的比值，我们可以将光子分为三类：

• 弹道光子：穿透深度 < 1 SMFP，保持原始方向
• 蛇形光子：5 SMFP < 穿透深度 < 10 SMFP，方向轻微偏转
• 散射光子：穿透深度 > 10 SMFP，方向完全随机

为了模拟这个过程，我们将采用蒙特卡洛方法——通过大量随机采样来近似物理过程。每个光子将被视为一个独立的"探针"，在其运动过程中根据概率决定是否发生散射。

蒙特卡洛模拟的核心在于：用随机数生成器代替物理定律，通过统计规律再现真实现象。

2. 构建Python模拟环境

2.1 初始化设置

首先导入必要的库并设置参数：

import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # 介质参数 SMFP = 1.0 # 散射平均自由程 medium_thickness = 15 # 介质厚度(SMFP单位) n_photons = 1000 # 模拟光子数量 # 散射参数 g = 0.9 # 各向异性因子，控制前向散射概率

2.2 光子步进模拟

我们定义一个函数来模拟单个光子的运动：

def simulate_photon(max_steps=1000): positions = [(0, 0, 0)] # 起始位置 direction = (0, 0, 1) # 初始方向(z轴正向) for _ in range(max_steps): # 计算下一步移动(固定步长=0.1*SMFP) step_length = 0.1 * SMFP new_pos = ( positions[-1][0] + direction[0]*step_length, positions[-1][1] + direction[1]*step_length, positions[-1][2] + direction[2]*step_length ) positions.append(new_pos) # 检查是否穿出介质 if new_pos[2] >= medium_thickness: break # 决定是否散射(指数分布) if np.random.exponential(scale=SMFP) < step_length: direction = scatter(direction, g) return np.array(positions)

2.3 散射模型实现

采用Henyey-Greenstein相位函数模拟散射方向：

def scatter(old_dir, g): """根据Henyey-Greenstein相位函数计算新方向""" # 计算散射角 cos_theta = (1 + g**2 - ((1 - g**2)/(1 - g + 2*g*np.random.random()))**2)/(2*g) theta = np.arccos(cos_theta) phi = 2 * np.pi * np.random.random() # 建立旋转矩阵 sin_theta = np.sin(theta) cos_phi = np.cos(phi) sin_phi = np.sin(phi) # 计算新方向 new_dir = ( sin_theta * cos_phi, sin_theta * sin_phi, cos_theta ) # 确保归一化 norm = np.sqrt(sum(x**2 for x in new_dir)) return tuple(x/norm for x in new_dir)

3. 可视化光子轨迹

3.1 单光子轨迹模拟

让我们先观察几个典型光子的运动路径：

# 模拟三个典型光子 ballistic = simulate_photon() # 弹道光子 snake = simulate_photon() # 蛇形光子 diffusive = simulate_photon() # 散射光子 # 绘制3D轨迹 fig = plt.figure(figsize=(12, 8)) ax = fig.add_subplot(111, projection='3d') # 弹道光子(红色) ax.plot(ballistic[:,0], ballistic[:,1], ballistic[:,2], 'r-', linewidth=2, label='Ballistic') # 蛇形光子(绿色) ax.plot(snake[:,0], snake[:,1], snake[:,2], 'g-', linewidth=1, label='Snake') # 散射光子(蓝色) ax.plot(diffusive[:,0], diffusive[:,1], diffusive[:,2], 'b-', alpha=0.5, label='Diffusive') ax.set_xlabel('X (SMFP)') ax.set_ylabel('Y (SMFP)') ax.set_zlabel('Z (SMFP)') ax.legend() plt.title('Photon Trajectories in Scattering Medium') plt.show()

3.2 群体统计特性

为了获得统计规律，我们需要模拟大量光子并分类：

def classify_photon(trajectory): """根据最终深度分类光子""" final_depth = trajectory[-1,2] if final_depth < 1.0 * SMFP: return "ballistic" elif final_depth < 10.0 * SMFP: return "snake" else: return "diffusive" # 模拟并分类1000个光子 results = {"ballistic": 0, "snake": 0, "diffusive": 0} trajectories = [] for _ in range(n_photons): traj = simulate_photon() category = classify_photon(traj) results[category] += 1 trajectories.append((traj, category)) # 打印统计结果 print("Photon distribution:") for cat, count in results.items(): print(f"{cat.capitalize()}: {count/n_photons*100:.1f}%")

典型输出可能如下：

4. 点扩散函数(PSF)分析

不同光子对成像的影响可以通过点扩散函数来量化。我们计算光子在出射面的分布：

def calculate_psf(trajectories, category): """计算指定类型光子的PSF""" exit_positions = [] for traj, cat in trajectories: if cat == category: exit_pos = traj[-1][:2] # 取x,y坐标 exit_positions.append(exit_pos) if not exit_positions: return None exit_positions = np.array(exit_positions) # 创建2D直方图 hist, xedges, yedges = np.histogram2d( exit_positions[:,0], exit_positions[:,1], bins=50, range=[[-5,5],[-5,5]] ) return hist # 计算各类PSF psf_ballistic = calculate_psf(trajectories, "ballistic") psf_snake = calculate_psf(trajectories, "snake") psf_diffusive = calculate_psf(trajectories, "diffusive") # 可视化 fig, axes = plt.subplots(1, 3, figsize=(15, 5)) titles = ["Ballistic PSF", "Snake PSF", "Diffusive PSF"] for ax, psf, title in zip(axes, [psf_ballistic, psf_snake, psf_diffusive], titles): if psf is not None: im = ax.imshow(psf.T, origin='lower', extent=[-5,5,-5,5], cmap='hot') ax.set_title(title) plt.colorbar(im, ax=ax)

从PSF图像可以明显看出：

• 弹道光子集中在中心点，保持原始信息
• 蛇形光子形成模糊光斑，但仍保留部分结构
• 散射光子完全随机分布，信息几乎丢失

5. 参数影响与优化思路

通过调整模拟参数，我们可以研究不同条件下的光子行为：

5.1 各向异性因子g的影响

g值决定了散射的前向性（g=1为完全前向，g=0为各向同性）：

g_values = [0.1, 0.5, 0.9] results_g = {} for g_val in g_values: # 临时修改全局g值 global g g = g_val # 重新模拟并分类 temp_results = {"ballistic": 0, "snake": 0, "diffusive": 0} for _ in range(500): traj = simulate_photon() cat = classify_photon(traj) temp_results[cat] += 1 results_g[g_val] = temp_results

数据可以整理为下表：

5.2 介质厚度的影响

固定其他参数，改变介质厚度：

thickness_values = [5, 10, 15, 20] results_thickness = {} for thick in thickness_values: global medium_thickness medium_thickness = thick temp_results = {"ballistic": 0, "snake": 0, "diffusive": 0} for _ in range(500): traj = simulate_photon() cat = classify_photon(traj) temp_results[cat] += 1 results_thickness[thick] = temp_results

可视化结果更直观：

# 提取数据准备绘图 thicknesses = sorted(results_thickness.keys()) ballistic_rates = [results_thickness[t]["ballistic"]/500*100 for t in thicknesses] snake_rates = [results_thickness[t]["snake"]/500*100 for t in thicknesses] diffusive_rates = [results_thickness[t]["diffusive"]/500*100 for t in thicknesses] plt.figure(figsize=(10,6)) plt.plot(thicknesses, ballistic_rates, 'ro-', label='Ballistic') plt.plot(thicknesses, snake_rates, 'gs-', label='Snake') plt.plot(thicknesses, diffusive_rates, 'b^-', label='Diffusive') plt.xlabel('Medium Thickness (SMFP)') plt.ylabel('Percentage (%)') plt.legend() plt.grid(True) plt.title('Photon Distribution vs Medium Thickness')

从曲线可以看出，随着介质增厚，弹道光子比例急剧下降，散射光子比例显著上升。

6. 实际应用与扩展方向

基于这个基础模拟框架，我们可以进一步开发更实用的功能：

6.1 时间分辨模拟

记录光子传输时间，模拟超短脉冲激光的时域响应：

def simulate_photon_with_time(max_steps=1000, c=1.0): """模拟光子并记录传输时间""" positions = [(0, 0, 0)] direction = (0, 0, 1) times = [0.0] # 时间记录(单位：SMFP/c) for _ in range(max_steps): step_length = 0.1 * SMFP new_pos = ( positions[-1][0] + direction[0]*step_length, positions[-1][1] + direction[1]*step_length, positions[-1][2] + direction[2]*step_length ) positions.append(new_pos) times.append(times[-1] + step_length/c) # 时间累加 if new_pos[2] >= medium_thickness: break if np.random.exponential(scale=SMFP) < step_length: direction = scatter(direction, g) return np.array(positions), np.array(times)

6.2 多波长模拟

不同波长的光子可能具有不同的散射特性：

def wavelength_dependent_smfp(wavelength, base_smfp=1.0, exponent=-4): """根据波长计算SMFP(假设遵循幂律关系)""" return base_smfp * (wavelength/500)**exponent # 500nm为参考波长 # 模拟红光(650nm)和蓝光(450nm)的差异 smfp_red = wavelength_dependent_smfp(650) smfp_blue = wavelength_dependent_smfp(450) print(f"Red light SMFP: {smfp_red:.2f}") print(f"Blue light SMFP: {smfp_blue:.2f}")

6.3 GPU加速

对于大规模模拟，可以使用CUDA加速：

# 示例：使用numba进行GPU加速 from numba import cuda @cuda.jit def gpu_photon_simulation(positions, directions, g, n_photons, max_steps): """在GPU上并行模拟多个光子""" i = cuda.grid(1) if i < n_photons: # 每个线程处理一个光子 pass # 实现略

在实际项目中，我发现将光子模拟并行化可以轻松获得百倍加速，特别是当需要模拟数百万光子时，GPU的优势更加明显。另一个实用技巧是在模拟前预先计算好随机数序列，可以进一步优化性能。