做光学实验的同学应该都有过这种经历一束激光打到毛玻璃或者生物组织上后面的屏幕上不是干净的光斑而是一片密密麻麻的颗粒状图案这就是散斑。更让人头疼的是当你想透过这块散射介质去看后面的物体图案乱成一团常规成像手段直接失效。但你如果试着把入射光束稍稍平移或偏转一点会发现出射的散斑图样在一定范围内几乎没变只是整体挪了个位置偏转幅度再大一些图样才慢慢失相关。这个“还能记住原光路信息”的有效范围就是光学记忆效应。这篇文章我就用 Python NumPy 从零搭一个二维散射模拟把散射介质简化成随机相位屏用角谱法做衍射传播最后同时给出散斑图样和记忆效应相关曲线并附完整可运行的代码。整个实现用到的只有 numpy 和 matplotlib 两个库非常适合刚接触光学模拟、又想快速看到物理图像的人。你不需要有很强的物理背景但最好对二维数组和傅里叶变换有个基本概念没有也没关系我会把每一步在算什么都讲清楚。1. 模拟之前先搞清楚散射和记忆效应到底是怎么一回事1.1 散射过程在光学里为什么难算光在介质中传播如果介质折射率分布不均匀波前就会被分割成无数个子波源每个子波源向四面八方发射彼此干涉最终在观察面上形成复杂的强度分布。这个过程的严格数学描述需要求解 Maxwell 方程组而实际样品生物组织、雾、毛玻璃的折射率分布又极难精确测量所以直接“硬算”几乎不可能。在实际模拟里我们很少去追每个粒子的散射而是用一个叫“随机相位屏”的模型来做等效把散射介质对光波的影响近似看成在光通过它时给波前叠加了一个随机空间相位。这个相位扰动会让原本平滑的波前变得皱巴巴的之后再经过一段自由空间衍射传播就能形成统计特征上和真实散斑很接近的图案。这个模型当然有一定局限性比如它没有考虑介质内部的多次散射和偏振变化但在描述“记忆效应”、散斑相关这类统计现象时它已经足够准确而且计算效率极高。用电脑做光学实验的优势就在这儿参数随便改改完重跑一次就行不用到光学平台上重新调光路。1.2 记忆效应散射背后的有序性如果散射完全是“混乱”的那散斑图样应该对入射条件极其敏感稍微动一下就面目全非。但实验和理论都发现在相当多散射介质的透射或反射光路里入射光发生小角度偏转或小范围平移时出射的散斑图样会保持很强相关性。这种“记住了原有照射条件”的现象就叫光学记忆效应。通俗点说散射介质像一面打碎的镜子碎片的走向看起来是随机的但如果你把整面镜子轻轻转一个小角度碎片的分布虽然整体变了但碎片之间的相对位置关系还是被“记住”了因此反射图案会有保留地跟着变化。记忆效应之所以重要是因为它给我们指了一条路散射不是完全不可逆的散射中仍然携带着入射光的信息只是被编码得很复杂。后面要讲的透过散射介质成像、波前整形本质上都是在利用这个“有序性”。1.3 模拟模型的取舍我为什么用相位屏而不是真的随机粒子搭建这个模拟时我第一个考虑的问题是要不要真的在二维网格里扔一堆随机折射率小球然后用逐点散射去算答案是不需要也不划算。真实粒子模型的边界条件非常麻烦散射角度、粒子大小、密度分布每个都要调算一次要很久新手很难判断结果到底对没对。随机相位屏虽然看起来“物理味”不够浓但它抓住了散射的核心效应——波前随机化。对一个薄散射样品来说光在样品内部发生的相位积累可以近似写成折射率起伏沿光路的积分这个积分结果在统计特性上就等价于一个随机相位分布。所以在二维模拟里把介质面替换成一张“随机相位图”是完全站得住脚的简化也是散射光学模拟里最常见的做法。2. 环境准备和数值参数设计别让 numpy 装不上拖后腿2.1 安装 numpy 和 matplotlib从命令行到 PyCharm先用最简单的命令行方式。打开终端确保当前 Python 环境是你平时写代码的那个然后执行pip install numpy matplotlib如果你想指定国内镜像加速可以用pip install numpy matplotlib -i https://pypi.tuna.tsinghua.edu.cn/simple安装完以后在 Python 里执行import numpy as np和import matplotlib.pyplot as plt不报错就说明基础环境没问题。如果你用的是 PyCharm更推荐在项目设置里安装。因为 PyCharm 每个项目通常都会创建一个独立虚拟环境你在系统命令行里装的包不一定能被这个项目看到。正确做法是点击File - Settings - Project - Python Interpreter在列表右侧点搜索 numpy 和 matplotlib然后 Install Package。这样装出来的库一定属于当前项目不会出现“明明装了却找不到”的问题。VSCode 用户也容易踩同一个坑右下角选择的 Python 解释器和你在终端里默认用的不是同一个。建议先创建一个.venv虚拟环境然后在 VSCode 命令面板里执行Python: Select Interpreter选到这个虚拟环境再在集成终端里pip install。2.2 为什么 PyCharm 显示有 numpyimport 却报 ModuleNotFoundError这个问题的原因几乎都是解释器错位。你很可能在 PyCharm 的某个虚拟环境里看到 numpy 已经安装了但当前运行脚本用的解释器却指向了另一个环境那个环境里没有 numpy。排查方法很简单在 PyCharm 里打开File - Settings - Project - Python Interpreter看当前环境列表再在编辑区最右下角看当前脚本实际使用的解释器确保两者一致。还有一个非常隐蔽的情况你的项目里有多个虚拟环境比如一个在.venv目录一个在系统 Python 里PyCharm 可能把系统 Python 当作默认解释器导致你在项目环境里装的包看不到。如果上面都排查过了还是报错就在 PyCharm 的 Python Console 里执行一句import sys, numpy print(numpy.__file__) print(sys.executable)如果sys.executable指向的路径和你在 PyCharm 设置里看到的解释器不一致直接把 Interpreter 重新选一遍、重启项目就好。这类环境错位问题占了新手报错里至少三分之一值得一次性搞清楚。2.3 网格、波长与 FFT 的关系物理量不是凭空来的光学模拟里网格参数是物理正确性的基础。我的数组尺寸是 N×N每个网格点的物理尺寸是 pixel那么整个模拟区域边长就是 L N × pixel。这里有个关键约束网格必须足够细能够采样到光场的最小结构。拿波长 λ 532 nm 来说这是绿色激光的典型波长。网格间距取 pixel 2 μm也就是 2 微米远大于波长但又不是大得离谱。光斑半径取 100 μm 左右在 512 × 512 的网格里大约占 100 个像素能保证散斑的细节被充分采样。观察距离 z 我取 0.05 m也就是 5 厘米这比较接近桌面光学实验里的常见光路尺寸。另外还要记住 FFT 的两个特性频域坐标的分辨率是 1/L也就是 1/N×pixel频域能表示的最大频率是 1/(2×pixel)对应 Nyquist 频率。计算角谱传播因子时所有频率分量都要落在这个范围内否则传播结果会出现混叠伪影。本文参数下最大空间频率对应的 λf 约等于 0.133远小于 1所以角谱里的开根号始终有实数解不会出现倏逝波的麻烦。3. 核心代码逐段拆解从随机相位屏到散斑图样3.1 随机相位屏怎么写先固定随机种子保证每次跑出来的散斑图一样方便复现rng np.random.default_rng(2024) phase rng.normal(0, 1.5, (N, N))这里生成的是均值为 0、标准差为 1.5 弧度的高斯随机相位。(N, N)是二维数组每个位置代表散射屏上那个点的相位扰动。标准差 1.5 弧度意味着相位扰动幅度不小足以让波前发生明显畸变但又不会大到让相位出现过于剧烈的跳变这样散斑统计比较接近实验。如果你把标准差调小到 0.1 rad光波基本感觉不到散射出射图案还是接近干净的光斑调大到 10 rad相位在 0 到 2π 之间疯狂跳动也能出散斑但会产生大量高频细节对网格分辨率要求更高。1.5 是我反复试下来视觉和计算都比较舒服的取值。3.2 高斯光束入射场入射光用最简单的高斯光束形式x (np.arange(N) - N // 2) * pixel x2d, y2d np.meshgrid(x, x) w0 60 * pixel E_in np.exp(-(x2d ** 2 y2d ** 2) / w0 ** 2)这里w0是高斯光束的束腰半径取 60 个像素也就是 120 μm。用减半再乘 pixel 的方式是为了让坐标原点落在数组中心。E_in 是二维复场分布现在虚部为 0只有振幅后续做倾斜或平移时会乘上额外的相位因子。为什么要用高斯光束而不是一个圆形光斑因为高斯光束的边缘是平滑衰减的在频域里带宽小不容易在和随机相位屏相乘后产生强烈的边缘衍射。圆形硬边光斑会在边界处产生一圈明显的 Fresnel 环干扰散斑的统计分布。3.3 用角谱法实现衍射传播衍射传播是整个过程的核心。我的做法是把经过相位屏的光场做二维 FFT把它分解成一系列平面波然后让每个平面波各自传播一段距离 z再逆 FFT 合成新的光场。频域坐标要和无 FFT 时的排列顺序对应。直接用np.fft.fftfreq生成的是符合 FFT 输出顺序的频率序列不需要再做 fftshift这样最后乘以传播因子时才不会错位fx np.fft.fftfreq(N, dpixel) fx2d, fy2d np.meshgrid(fx, fx) k0 2 * np.pi / lam prop np.exp(1j * k0 * z * np.sqrt(1 - (lam * fx2d) ** 2 - (lam * fy2d) ** 2))传播因子的物理含义是频率为 (fx, fy) 的平面波沿 z 方向传播距离 z 后相位变化为 k0 * z * sqrt(1 - (λfx)^2 - (λfy)^2)。这里 λfx 和 λfy 代表平面波方向余弦的横向分量开根号后是纵向方向余弦。对近轴光来说这个因子非常接近 1但保留完整形式能让模型在更大角度下也成立。接下来定义单次“照明散射传播”的过程def propagate(E_in, phase, prop): E_screen E_in * np.exp(1j * phase) E_out np.fft.ifft2(np.fft.fft2(E_screen) * prop) return E_out先给入射场乘上随机相位屏等效于光场透过散射介质然后 FFT 到频域乘传播因子再逆 FFT得到观察面上的复场。强度分布就是复场模长的平方def intensity(E): return np.abs(E) ** 23.4 完整代码画出散斑图与记忆效应曲线合并上面的所有片段加上记忆效应计算的循环就得到下面这份可以直接运行的完整代码。我会在下一节详细解释记忆效应那段循环在做什么。import numpy as np import matplotlib.pyplot as plt # 物理参数 lam 532e-9 # 波长 532 nm pixel 2e-6 # 网格采样间距 2 μm N 512 # 网格大小 512x512 z 0.05 # 散射屏到观察面距离 5 cm # 实空间坐标 x (np.arange(N) - N // 2) * pixel x2d, y2d np.meshgrid(x, x) # 频域坐标注意顺序要和 fft2 输出对应 fx np.fft.fftfreq(N, dpixel) fx2d, fy2d np.meshgrid(fx, fx) # 角谱传播因子 k0 2 * np.pi / lam prop np.exp(1j * k0 * z * np.sqrt(1 - (lam * fx2d) ** 2 - (lam * fy2d) ** 2)) # 随机相位屏标准差 1.5 rad rng np.random.default_rng(2024) phase rng.normal(0, 1.5, (N, N)) # 入射高斯光束 w0 60 * pixel E_in np.exp(-(x2d ** 2 y2d ** 2) / w0 ** 2) # 传播函数 def propagate(E_in): E_screen E_in * np.exp(1j * phase) E_out np.fft.ifft2(np.fft.fft2(E_screen) * prop) return E_out # 参考散斑图 E_ref propagate(E_in) I_ref np.abs(E_ref) ** 2 # 计算记忆效应曲线横向平移入射光束观察散斑失相关速度 shifts np.arange(-30, 31, 2) # 平移量单位像素 corrs [] # 限定的相关计算区域光斑附近的圆形掩模避免边缘空区域干扰 mask (x2d ** 2 y2d ** 2) (100 * pixel) ** 2 for s in shifts: E_shift np.roll(E_in, s, axis1) # 横向移动入射光束 I_shift np.abs(propagate(E_shift)) ** 2 c np.corrcoef(I_ref[mask], I_shift[mask])[0, 1] corrs.append(c) # 绘制结果 plt.figure(figsize(12, 4)) plt.subplot(1, 3, 1) plt.imshow(I_ref, cmapinferno) plt.title(Speckle pattern) plt.colorbar() plt.subplot(1, 3, 2) plt.plot(shifts * pixel * 1e6, corrs, o-) plt.xlabel(Lateral shift (um)) plt.ylabel(Correlation) plt.title(Memory effect curve) plt.grid(True) plt.subplot(1, 3, 3) plt.imshow(np.sqrt(I_ref), cmapinferno) plt.title(Speckle amplitude) plt.colorbar() plt.tight_layout() plt.show()这份代码在我机器上跑一遍大约需要十几秒主要时间花在记忆效应循环里对每个平移量做一次完整的 FFT 传播。如果你想跑得更快可以把 N 从 512 降到 256曲线形状基本不变。4. 记忆效应曲线的完整计算如何把“相似程度”量化出来4.1 散斑图解读你看到的并不是噪声跑完代码第一张图是观察面上的散斑强度图。它看起来像随机噪声但其实每个亮斑都有自己的物理尺寸也就是“散斑颗粒”的大小。散斑平均尺寸大约正比于 λz/DD 是照明光斑的直径。在这个参数下λz/D 大约一百多微米换算成像素大概几十个像素所以你看到的亮斑不是单像素闪烁而是一小块一小块的。散斑图有一个重要统计特性强度分布服从负指数分布也就是大部分地方是暗的少数亮斑强度特别高。你可以用np.histogram(I_ref.ravel(), bins100)看一眼会发现暗像素占比很高这和真实激光散斑的统计规律是一致的。如果你做出来的散斑图亮得像块白板多半是相位屏方差太小散射不够强。4.2 相关函数计算为什么不能直接比整张图记忆效应曲线计算的核心是比较参考散斑图 I_ref 和平移后散斑图 I_shift 的相似度。直接在整个 512×512 数组上算相关系数有一个隐患光斑只照亮了画面中央一小块区域其余地方都是接近 0 的黑色背景。黑色背景参与计算会强行拉高相关系数因为无论怎么平移黑色背景都长一个样最后曲线看起来就像一条接近 1 的直线失去区分度。所以我加了圆形掩模 mask只取中心半径 100 μm 范围内的像素做相关。这个范围比照明光斑略大能保证覆盖到主要散斑区域又不会把大片背景卷进来。相关系数用np.corrcoef的计算公式是[ C \frac{\sum (I_{ref} - \bar{I}{ref})(I{shift} - \bar{I}{shift})} {\sqrt{\sum (I{ref} - \bar{I}{ref})^2 \sum (I{shift} - \bar{I}_{shift})^2}} ]这个公式衡量的是两组强度值的线性相关程度。当两张散斑图完全一样时C 等于 1完全无关时C 接近 0。注意散斑强度不服从高斯分布所以这是一种“基于强度的线性相关”不是严格的散斑对比度测量但对于教学演示和定性判断记忆效应范围足够直观。4.3 结果如何与记忆效应理论对应跑出来的记忆效应曲线典型形状是中间一个峰值接近 1两侧快速下降最后在 0 附近波动。峰值半高宽就是“记忆范围”的直观体现入射光束横向平移在这个范围以内散斑图样能保持较强相似性超出这个范围图样才真正“翻篇”。我之前用一组参数实测曲线半高全宽大约在几十微米量级。这个结果和散斑颗粒尺寸在一个量级符合直觉你把照明光斑移走一个散斑颗粒的距离原来那一小块散斑区域已经被完全不同的光场覆盖了相关性自然大幅下降。要注意的是真实散射介质的记忆效应范围还会受介质厚度、散射各向异性、光路类型透射还是反射影响单层相位屏模型更接近“薄散射体一定传播距离”的情形它能复现出“相关随位移衰减”的核心现象但不要指望它精确预测所有实验数据。5. 参数扫描实验散射强度和传播距离如何影响记忆范围5.1 改变相位屏方差散射越强记忆范围越小吗这是我最想让你亲手试的参数实验。把相位屏标准差 sigma 依次取 0.5, 1.5, 3.0分别计算记忆效应曲线画在一起results {} for sigma in [0.5, 1.5, 3.0]: rng np.random.default_rng(2024) phase_local rng.normal(0, sigma, (N, N)) corrs_local [] for s in shifts: E_shift np.roll(E_in, s, axis1) E_screen E_shift * np.exp(1j * phase_local) I_shift np.abs(np.fft.ifft2(np.fft.fft2(E_screen) * prop)) ** 2 corrs_local.append(np.corrcoef(I_ref_local[mask], I_shift[mask])[0, 1]) results[sigma] corrs_local注意一个细节每改一次 sigma都固定随机种子但相位屏本身是重新生成的所以不能用同一个相位屏直接改大小你比较的是“不同散射强度”下的统计行为而不是同一屏的强度缩放。实验结论一般是这样sigma 很小时光场几乎没被扰乱平移一点之后图案依然高度相似曲线很宽sigma 增大后散斑颗粒变小、相位变化更剧烈记忆范围收窄曲线半高全宽显著减小。这符合直觉也符合“散射越强信息被搅拌得越碎能记住原路线的角度范围越小”的物理图像。5.2 改变传播距离散斑放大和失相关的博弈把观察距离 z 从 0.02 m 改到 0.1 m你会发现散斑颗粒明显变大因为 λz/D 变大了。散斑颗粒变大意味着观察平面上相邻区域强度变化变慢理论上记忆效应曲线会变宽一些但同时传播过程会让不同角度的散射光进一步分开造成额外的失相关。实际扫描下来z 增大时曲线峰值还是会下降得更快一点因为传播距离越远原来一个位置的光被分散得越开后续照明位置的散斑图样和参考图的交叠区域变小。这个现象在真实实验里也有对应散射介质离探测器越远透过它成像的记忆效应角度范围往往越窄。5.3 掩模尺寸对曲线的隐藏影响你也许没注意到mask 的半径也会显著改变曲线形状。mask 取太小只取光斑中心一小块相关系数会因为统计样本太少而抖动剧烈mask 取太大把大量暗背景卷进来曲线会被人为抬高看起来好像记忆效应范围很大。一个相对稳的做法mask 半径取照明光斑半径的 1.5 到 2 倍既能覆盖主要光场又不会混入过多背景。这是我在跑模拟时踩过的小坑分享给你免得你以后对着一条“异常平坦”的曲线发呆。6. 从模拟到应用记忆效应能用来做什么6.1 透过散射介质成像的基本逻辑如果散射介质能“记住”入射角度那我们就可以利用这个性质做扫描成像。思路很直接把一束聚焦光打到散射介质上在出射面用相机记录散斑图样然后改变入射角度再记一张散斑图。因为记忆效应的存在入射角在小范围内变化时出射散斑图样只发生整体平移几乎不改变图案结构。这样我们可以通过反卷积从不同角的散斑图里恢复出隐藏在散射介质后面的物体信息。这在生物组织成像中特别有吸引力因为组织是强散射介质传统显微镜根本看不深记忆效应给了我们一个“虽然看不清但散射里藏着可提取的信号”的窗口。6.2 波前整形用反馈优化对抗散射另一个应用方向是波前整形。既然散射过程可以理解为光场被随机相位调制了那我们在入射端加一个空间光调制器给它加载一个精心设计的相位图让经过散射介质后的光在目标位置重新汇聚成一个焦点。这个过程通常用一个迭代反馈算法不断微调相位屏使目标区域的光强达到最大。模拟里复现这个过程的成本并不高你把入射场乘上一个可调节的相位掩模用本文的 propagate 函数算目标点的强度然后用遗传算法或梯度下降法迭代调相位。我第一次把焦点调出来的时候感觉就像在迷宫里找到了一条隐藏路径——散射本来是不可预测的但一旦你把它当成一个可以被“驯服”的线性系统问题就变成了优化问题。6.3 扩展方向从相位屏到多层散射模型单层随机相位屏毕竟是简化版它默认散射全部发生在同一个平面上。真实生物组织、泡沫、乳浊液里的散射是发生在三维体积内的这时记忆效应范围会随介质厚度进一步收窄散斑统计也复杂得多。想走得更远可以把单层相位屏扩展成多层在几个 z 位置分别放置随机相位屏中间用角谱法一段一段传播。这种“多层相位屏模型”既保留了计算效率又能模拟出体积散射的很多关键特征。代码结构也很好改只要把本文里的 propagate 过程包到一个循环里即可for phase_layer in phase_layers: E E * np.exp(1j * phase_layer) E np.fft.ifft2(np.fft.fft2(E) * prop_layer)这种扩展会让记忆效应相关曲线的形状更接近真实实验数据也会让你对散射成像的理解更立体。我最初也只是抱着“把散斑画出来看看”的心态写下这份代码结果越调越觉得有意思。现在回头看这个模拟的价值不止在于复现一个现象它更是一个很好用的物理直觉训练场网格怎么取、掩模怎么选、相关怎么算、曲线为什么长这样每一个选择背后都有物理和数值的约束。你把这份代码从头到尾吃透以后再去翻那些讲散射成像、散斑相关的论文会发现很多公式和图表突然就变得亲切了。最后再分享一个小技巧跑通之后把相位屏的标准差和传播距离都做成可调节参数做一次小规模参数扫描你对记忆效应的理解会比我当初看十篇文章还要牢。
用Python模拟光学记忆效应：基于随机相位屏与角谱法的散斑分析
做光学实验的同学应该都有过这种经历一束激光打到毛玻璃或者生物组织上后面的屏幕上不是干净的光斑而是一片密密麻麻的颗粒状图案这就是散斑。更让人头疼的是当你想透过这块散射介质去看后面的物体图案乱成一团常规成像手段直接失效。但你如果试着把入射光束稍稍平移或偏转一点会发现出射的散斑图样在一定范围内几乎没变只是整体挪了个位置偏转幅度再大一些图样才慢慢失相关。这个“还能记住原光路信息”的有效范围就是光学记忆效应。这篇文章我就用 Python NumPy 从零搭一个二维散射模拟把散射介质简化成随机相位屏用角谱法做衍射传播最后同时给出散斑图样和记忆效应相关曲线并附完整可运行的代码。整个实现用到的只有 numpy 和 matplotlib 两个库非常适合刚接触光学模拟、又想快速看到物理图像的人。你不需要有很强的物理背景但最好对二维数组和傅里叶变换有个基本概念没有也没关系我会把每一步在算什么都讲清楚。1. 模拟之前先搞清楚散射和记忆效应到底是怎么一回事1.1 散射过程在光学里为什么难算光在介质中传播如果介质折射率分布不均匀波前就会被分割成无数个子波源每个子波源向四面八方发射彼此干涉最终在观察面上形成复杂的强度分布。这个过程的严格数学描述需要求解 Maxwell 方程组而实际样品生物组织、雾、毛玻璃的折射率分布又极难精确测量所以直接“硬算”几乎不可能。在实际模拟里我们很少去追每个粒子的散射而是用一个叫“随机相位屏”的模型来做等效把散射介质对光波的影响近似看成在光通过它时给波前叠加了一个随机空间相位。这个相位扰动会让原本平滑的波前变得皱巴巴的之后再经过一段自由空间衍射传播就能形成统计特征上和真实散斑很接近的图案。这个模型当然有一定局限性比如它没有考虑介质内部的多次散射和偏振变化但在描述“记忆效应”、散斑相关这类统计现象时它已经足够准确而且计算效率极高。用电脑做光学实验的优势就在这儿参数随便改改完重跑一次就行不用到光学平台上重新调光路。1.2 记忆效应散射背后的有序性如果散射完全是“混乱”的那散斑图样应该对入射条件极其敏感稍微动一下就面目全非。但实验和理论都发现在相当多散射介质的透射或反射光路里入射光发生小角度偏转或小范围平移时出射的散斑图样会保持很强相关性。这种“记住了原有照射条件”的现象就叫光学记忆效应。通俗点说散射介质像一面打碎的镜子碎片的走向看起来是随机的但如果你把整面镜子轻轻转一个小角度碎片的分布虽然整体变了但碎片之间的相对位置关系还是被“记住”了因此反射图案会有保留地跟着变化。记忆效应之所以重要是因为它给我们指了一条路散射不是完全不可逆的散射中仍然携带着入射光的信息只是被编码得很复杂。后面要讲的透过散射介质成像、波前整形本质上都是在利用这个“有序性”。1.3 模拟模型的取舍我为什么用相位屏而不是真的随机粒子搭建这个模拟时我第一个考虑的问题是要不要真的在二维网格里扔一堆随机折射率小球然后用逐点散射去算答案是不需要也不划算。真实粒子模型的边界条件非常麻烦散射角度、粒子大小、密度分布每个都要调算一次要很久新手很难判断结果到底对没对。随机相位屏虽然看起来“物理味”不够浓但它抓住了散射的核心效应——波前随机化。对一个薄散射样品来说光在样品内部发生的相位积累可以近似写成折射率起伏沿光路的积分这个积分结果在统计特性上就等价于一个随机相位分布。所以在二维模拟里把介质面替换成一张“随机相位图”是完全站得住脚的简化也是散射光学模拟里最常见的做法。2. 环境准备和数值参数设计别让 numpy 装不上拖后腿2.1 安装 numpy 和 matplotlib从命令行到 PyCharm先用最简单的命令行方式。打开终端确保当前 Python 环境是你平时写代码的那个然后执行pip install numpy matplotlib如果你想指定国内镜像加速可以用pip install numpy matplotlib -i https://pypi.tuna.tsinghua.edu.cn/simple安装完以后在 Python 里执行import numpy as np和import matplotlib.pyplot as plt不报错就说明基础环境没问题。如果你用的是 PyCharm更推荐在项目设置里安装。因为 PyCharm 每个项目通常都会创建一个独立虚拟环境你在系统命令行里装的包不一定能被这个项目看到。正确做法是点击File - Settings - Project - Python Interpreter在列表右侧点搜索 numpy 和 matplotlib然后 Install Package。这样装出来的库一定属于当前项目不会出现“明明装了却找不到”的问题。VSCode 用户也容易踩同一个坑右下角选择的 Python 解释器和你在终端里默认用的不是同一个。建议先创建一个.venv虚拟环境然后在 VSCode 命令面板里执行Python: Select Interpreter选到这个虚拟环境再在集成终端里pip install。2.2 为什么 PyCharm 显示有 numpyimport 却报 ModuleNotFoundError这个问题的原因几乎都是解释器错位。你很可能在 PyCharm 的某个虚拟环境里看到 numpy 已经安装了但当前运行脚本用的解释器却指向了另一个环境那个环境里没有 numpy。排查方法很简单在 PyCharm 里打开File - Settings - Project - Python Interpreter看当前环境列表再在编辑区最右下角看当前脚本实际使用的解释器确保两者一致。还有一个非常隐蔽的情况你的项目里有多个虚拟环境比如一个在.venv目录一个在系统 Python 里PyCharm 可能把系统 Python 当作默认解释器导致你在项目环境里装的包看不到。如果上面都排查过了还是报错就在 PyCharm 的 Python Console 里执行一句import sys, numpy print(numpy.__file__) print(sys.executable)如果sys.executable指向的路径和你在 PyCharm 设置里看到的解释器不一致直接把 Interpreter 重新选一遍、重启项目就好。这类环境错位问题占了新手报错里至少三分之一值得一次性搞清楚。2.3 网格、波长与 FFT 的关系物理量不是凭空来的光学模拟里网格参数是物理正确性的基础。我的数组尺寸是 N×N每个网格点的物理尺寸是 pixel那么整个模拟区域边长就是 L N × pixel。这里有个关键约束网格必须足够细能够采样到光场的最小结构。拿波长 λ 532 nm 来说这是绿色激光的典型波长。网格间距取 pixel 2 μm也就是 2 微米远大于波长但又不是大得离谱。光斑半径取 100 μm 左右在 512 × 512 的网格里大约占 100 个像素能保证散斑的细节被充分采样。观察距离 z 我取 0.05 m也就是 5 厘米这比较接近桌面光学实验里的常见光路尺寸。另外还要记住 FFT 的两个特性频域坐标的分辨率是 1/L也就是 1/N×pixel频域能表示的最大频率是 1/(2×pixel)对应 Nyquist 频率。计算角谱传播因子时所有频率分量都要落在这个范围内否则传播结果会出现混叠伪影。本文参数下最大空间频率对应的 λf 约等于 0.133远小于 1所以角谱里的开根号始终有实数解不会出现倏逝波的麻烦。3. 核心代码逐段拆解从随机相位屏到散斑图样3.1 随机相位屏怎么写先固定随机种子保证每次跑出来的散斑图一样方便复现rng np.random.default_rng(2024) phase rng.normal(0, 1.5, (N, N))这里生成的是均值为 0、标准差为 1.5 弧度的高斯随机相位。(N, N)是二维数组每个位置代表散射屏上那个点的相位扰动。标准差 1.5 弧度意味着相位扰动幅度不小足以让波前发生明显畸变但又不会大到让相位出现过于剧烈的跳变这样散斑统计比较接近实验。如果你把标准差调小到 0.1 rad光波基本感觉不到散射出射图案还是接近干净的光斑调大到 10 rad相位在 0 到 2π 之间疯狂跳动也能出散斑但会产生大量高频细节对网格分辨率要求更高。1.5 是我反复试下来视觉和计算都比较舒服的取值。3.2 高斯光束入射场入射光用最简单的高斯光束形式x (np.arange(N) - N // 2) * pixel x2d, y2d np.meshgrid(x, x) w0 60 * pixel E_in np.exp(-(x2d ** 2 y2d ** 2) / w0 ** 2)这里w0是高斯光束的束腰半径取 60 个像素也就是 120 μm。用减半再乘 pixel 的方式是为了让坐标原点落在数组中心。E_in 是二维复场分布现在虚部为 0只有振幅后续做倾斜或平移时会乘上额外的相位因子。为什么要用高斯光束而不是一个圆形光斑因为高斯光束的边缘是平滑衰减的在频域里带宽小不容易在和随机相位屏相乘后产生强烈的边缘衍射。圆形硬边光斑会在边界处产生一圈明显的 Fresnel 环干扰散斑的统计分布。3.3 用角谱法实现衍射传播衍射传播是整个过程的核心。我的做法是把经过相位屏的光场做二维 FFT把它分解成一系列平面波然后让每个平面波各自传播一段距离 z再逆 FFT 合成新的光场。频域坐标要和无 FFT 时的排列顺序对应。直接用np.fft.fftfreq生成的是符合 FFT 输出顺序的频率序列不需要再做 fftshift这样最后乘以传播因子时才不会错位fx np.fft.fftfreq(N, dpixel) fx2d, fy2d np.meshgrid(fx, fx) k0 2 * np.pi / lam prop np.exp(1j * k0 * z * np.sqrt(1 - (lam * fx2d) ** 2 - (lam * fy2d) ** 2))传播因子的物理含义是频率为 (fx, fy) 的平面波沿 z 方向传播距离 z 后相位变化为 k0 * z * sqrt(1 - (λfx)^2 - (λfy)^2)。这里 λfx 和 λfy 代表平面波方向余弦的横向分量开根号后是纵向方向余弦。对近轴光来说这个因子非常接近 1但保留完整形式能让模型在更大角度下也成立。接下来定义单次“照明散射传播”的过程def propagate(E_in, phase, prop): E_screen E_in * np.exp(1j * phase) E_out np.fft.ifft2(np.fft.fft2(E_screen) * prop) return E_out先给入射场乘上随机相位屏等效于光场透过散射介质然后 FFT 到频域乘传播因子再逆 FFT得到观察面上的复场。强度分布就是复场模长的平方def intensity(E): return np.abs(E) ** 23.4 完整代码画出散斑图与记忆效应曲线合并上面的所有片段加上记忆效应计算的循环就得到下面这份可以直接运行的完整代码。我会在下一节详细解释记忆效应那段循环在做什么。import numpy as np import matplotlib.pyplot as plt # 物理参数 lam 532e-9 # 波长 532 nm pixel 2e-6 # 网格采样间距 2 μm N 512 # 网格大小 512x512 z 0.05 # 散射屏到观察面距离 5 cm # 实空间坐标 x (np.arange(N) - N // 2) * pixel x2d, y2d np.meshgrid(x, x) # 频域坐标注意顺序要和 fft2 输出对应 fx np.fft.fftfreq(N, dpixel) fx2d, fy2d np.meshgrid(fx, fx) # 角谱传播因子 k0 2 * np.pi / lam prop np.exp(1j * k0 * z * np.sqrt(1 - (lam * fx2d) ** 2 - (lam * fy2d) ** 2)) # 随机相位屏标准差 1.5 rad rng np.random.default_rng(2024) phase rng.normal(0, 1.5, (N, N)) # 入射高斯光束 w0 60 * pixel E_in np.exp(-(x2d ** 2 y2d ** 2) / w0 ** 2) # 传播函数 def propagate(E_in): E_screen E_in * np.exp(1j * phase) E_out np.fft.ifft2(np.fft.fft2(E_screen) * prop) return E_out # 参考散斑图 E_ref propagate(E_in) I_ref np.abs(E_ref) ** 2 # 计算记忆效应曲线横向平移入射光束观察散斑失相关速度 shifts np.arange(-30, 31, 2) # 平移量单位像素 corrs [] # 限定的相关计算区域光斑附近的圆形掩模避免边缘空区域干扰 mask (x2d ** 2 y2d ** 2) (100 * pixel) ** 2 for s in shifts: E_shift np.roll(E_in, s, axis1) # 横向移动入射光束 I_shift np.abs(propagate(E_shift)) ** 2 c np.corrcoef(I_ref[mask], I_shift[mask])[0, 1] corrs.append(c) # 绘制结果 plt.figure(figsize(12, 4)) plt.subplot(1, 3, 1) plt.imshow(I_ref, cmapinferno) plt.title(Speckle pattern) plt.colorbar() plt.subplot(1, 3, 2) plt.plot(shifts * pixel * 1e6, corrs, o-) plt.xlabel(Lateral shift (um)) plt.ylabel(Correlation) plt.title(Memory effect curve) plt.grid(True) plt.subplot(1, 3, 3) plt.imshow(np.sqrt(I_ref), cmapinferno) plt.title(Speckle amplitude) plt.colorbar() plt.tight_layout() plt.show()这份代码在我机器上跑一遍大约需要十几秒主要时间花在记忆效应循环里对每个平移量做一次完整的 FFT 传播。如果你想跑得更快可以把 N 从 512 降到 256曲线形状基本不变。4. 记忆效应曲线的完整计算如何把“相似程度”量化出来4.1 散斑图解读你看到的并不是噪声跑完代码第一张图是观察面上的散斑强度图。它看起来像随机噪声但其实每个亮斑都有自己的物理尺寸也就是“散斑颗粒”的大小。散斑平均尺寸大约正比于 λz/DD 是照明光斑的直径。在这个参数下λz/D 大约一百多微米换算成像素大概几十个像素所以你看到的亮斑不是单像素闪烁而是一小块一小块的。散斑图有一个重要统计特性强度分布服从负指数分布也就是大部分地方是暗的少数亮斑强度特别高。你可以用np.histogram(I_ref.ravel(), bins100)看一眼会发现暗像素占比很高这和真实激光散斑的统计规律是一致的。如果你做出来的散斑图亮得像块白板多半是相位屏方差太小散射不够强。4.2 相关函数计算为什么不能直接比整张图记忆效应曲线计算的核心是比较参考散斑图 I_ref 和平移后散斑图 I_shift 的相似度。直接在整个 512×512 数组上算相关系数有一个隐患光斑只照亮了画面中央一小块区域其余地方都是接近 0 的黑色背景。黑色背景参与计算会强行拉高相关系数因为无论怎么平移黑色背景都长一个样最后曲线看起来就像一条接近 1 的直线失去区分度。所以我加了圆形掩模 mask只取中心半径 100 μm 范围内的像素做相关。这个范围比照明光斑略大能保证覆盖到主要散斑区域又不会把大片背景卷进来。相关系数用np.corrcoef的计算公式是[ C \frac{\sum (I_{ref} - \bar{I}{ref})(I{shift} - \bar{I}{shift})} {\sqrt{\sum (I{ref} - \bar{I}{ref})^2 \sum (I{shift} - \bar{I}_{shift})^2}} ]这个公式衡量的是两组强度值的线性相关程度。当两张散斑图完全一样时C 等于 1完全无关时C 接近 0。注意散斑强度不服从高斯分布所以这是一种“基于强度的线性相关”不是严格的散斑对比度测量但对于教学演示和定性判断记忆效应范围足够直观。4.3 结果如何与记忆效应理论对应跑出来的记忆效应曲线典型形状是中间一个峰值接近 1两侧快速下降最后在 0 附近波动。峰值半高宽就是“记忆范围”的直观体现入射光束横向平移在这个范围以内散斑图样能保持较强相似性超出这个范围图样才真正“翻篇”。我之前用一组参数实测曲线半高全宽大约在几十微米量级。这个结果和散斑颗粒尺寸在一个量级符合直觉你把照明光斑移走一个散斑颗粒的距离原来那一小块散斑区域已经被完全不同的光场覆盖了相关性自然大幅下降。要注意的是真实散射介质的记忆效应范围还会受介质厚度、散射各向异性、光路类型透射还是反射影响单层相位屏模型更接近“薄散射体一定传播距离”的情形它能复现出“相关随位移衰减”的核心现象但不要指望它精确预测所有实验数据。5. 参数扫描实验散射强度和传播距离如何影响记忆范围5.1 改变相位屏方差散射越强记忆范围越小吗这是我最想让你亲手试的参数实验。把相位屏标准差 sigma 依次取 0.5, 1.5, 3.0分别计算记忆效应曲线画在一起results {} for sigma in [0.5, 1.5, 3.0]: rng np.random.default_rng(2024) phase_local rng.normal(0, sigma, (N, N)) corrs_local [] for s in shifts: E_shift np.roll(E_in, s, axis1) E_screen E_shift * np.exp(1j * phase_local) I_shift np.abs(np.fft.ifft2(np.fft.fft2(E_screen) * prop)) ** 2 corrs_local.append(np.corrcoef(I_ref_local[mask], I_shift[mask])[0, 1]) results[sigma] corrs_local注意一个细节每改一次 sigma都固定随机种子但相位屏本身是重新生成的所以不能用同一个相位屏直接改大小你比较的是“不同散射强度”下的统计行为而不是同一屏的强度缩放。实验结论一般是这样sigma 很小时光场几乎没被扰乱平移一点之后图案依然高度相似曲线很宽sigma 增大后散斑颗粒变小、相位变化更剧烈记忆范围收窄曲线半高全宽显著减小。这符合直觉也符合“散射越强信息被搅拌得越碎能记住原路线的角度范围越小”的物理图像。5.2 改变传播距离散斑放大和失相关的博弈把观察距离 z 从 0.02 m 改到 0.1 m你会发现散斑颗粒明显变大因为 λz/D 变大了。散斑颗粒变大意味着观察平面上相邻区域强度变化变慢理论上记忆效应曲线会变宽一些但同时传播过程会让不同角度的散射光进一步分开造成额外的失相关。实际扫描下来z 增大时曲线峰值还是会下降得更快一点因为传播距离越远原来一个位置的光被分散得越开后续照明位置的散斑图样和参考图的交叠区域变小。这个现象在真实实验里也有对应散射介质离探测器越远透过它成像的记忆效应角度范围往往越窄。5.3 掩模尺寸对曲线的隐藏影响你也许没注意到mask 的半径也会显著改变曲线形状。mask 取太小只取光斑中心一小块相关系数会因为统计样本太少而抖动剧烈mask 取太大把大量暗背景卷进来曲线会被人为抬高看起来好像记忆效应范围很大。一个相对稳的做法mask 半径取照明光斑半径的 1.5 到 2 倍既能覆盖主要光场又不会混入过多背景。这是我在跑模拟时踩过的小坑分享给你免得你以后对着一条“异常平坦”的曲线发呆。6. 从模拟到应用记忆效应能用来做什么6.1 透过散射介质成像的基本逻辑如果散射介质能“记住”入射角度那我们就可以利用这个性质做扫描成像。思路很直接把一束聚焦光打到散射介质上在出射面用相机记录散斑图样然后改变入射角度再记一张散斑图。因为记忆效应的存在入射角在小范围内变化时出射散斑图样只发生整体平移几乎不改变图案结构。这样我们可以通过反卷积从不同角的散斑图里恢复出隐藏在散射介质后面的物体信息。这在生物组织成像中特别有吸引力因为组织是强散射介质传统显微镜根本看不深记忆效应给了我们一个“虽然看不清但散射里藏着可提取的信号”的窗口。6.2 波前整形用反馈优化对抗散射另一个应用方向是波前整形。既然散射过程可以理解为光场被随机相位调制了那我们在入射端加一个空间光调制器给它加载一个精心设计的相位图让经过散射介质后的光在目标位置重新汇聚成一个焦点。这个过程通常用一个迭代反馈算法不断微调相位屏使目标区域的光强达到最大。模拟里复现这个过程的成本并不高你把入射场乘上一个可调节的相位掩模用本文的 propagate 函数算目标点的强度然后用遗传算法或梯度下降法迭代调相位。我第一次把焦点调出来的时候感觉就像在迷宫里找到了一条隐藏路径——散射本来是不可预测的但一旦你把它当成一个可以被“驯服”的线性系统问题就变成了优化问题。6.3 扩展方向从相位屏到多层散射模型单层随机相位屏毕竟是简化版它默认散射全部发生在同一个平面上。真实生物组织、泡沫、乳浊液里的散射是发生在三维体积内的这时记忆效应范围会随介质厚度进一步收窄散斑统计也复杂得多。想走得更远可以把单层相位屏扩展成多层在几个 z 位置分别放置随机相位屏中间用角谱法一段一段传播。这种“多层相位屏模型”既保留了计算效率又能模拟出体积散射的很多关键特征。代码结构也很好改只要把本文里的 propagate 过程包到一个循环里即可for phase_layer in phase_layers: E E * np.exp(1j * phase_layer) E np.fft.ifft2(np.fft.fft2(E) * prop_layer)这种扩展会让记忆效应相关曲线的形状更接近真实实验数据也会让你对散射成像的理解更立体。我最初也只是抱着“把散斑画出来看看”的心态写下这份代码结果越调越觉得有意思。现在回头看这个模拟的价值不止在于复现一个现象它更是一个很好用的物理直觉训练场网格怎么取、掩模怎么选、相关怎么算、曲线为什么长这样每一个选择背后都有物理和数值的约束。你把这份代码从头到尾吃透以后再去翻那些讲散射成像、散斑相关的论文会发现很多公式和图表突然就变得亲切了。最后再分享一个小技巧跑通之后把相位屏的标准差和传播距离都做成可调节参数做一次小规模参数扫描你对记忆效应的理解会比我当初看十篇文章还要牢。