基于Kretschmann结构的双波长SPR强度调制MATLAB仿真

简介这份PDF文档围绕双波长强度调制表面等离子体共振SPR传感器的设计方案与仿真验证展开适合光学检测、生物医学传感、化学分析与环境监测方向的研究人员和相关专业高年级学生阅读。内容从SPR物理光学现象入手详细介绍了基于Kretschmann棱镜耦合的四层介质反射模型与菲涅耳公式推导并重点讲解了使用两个光纤滤波器从1550nm附近ASE光源中选取1540nm和1560nm作为双波长光源、通过反射光强度差值实现测量的改进型强度调制方法。文档还利用MATLAB程序模拟了不同入射角下多种气体的折射率拟合曲线展示了该技术在降低光源稳定性要求、扩大测量范围等方面的优势同时给出了K9半圆柱棱镜、铬膜与金膜等实验装置设计及结果讨论。整个压缩包内包含1个PDF文件大小约401KB以理论分析、公式推导和仿真结果为主。目前已有115人学习下载适合需要快速建立双波长SPR传感器理论框架并获取MATLAB建模思路的读者。1. 双波长强度调制不是新器件而是对传统SPR测量方式的一次减法做气体折射率检测时传统强度调制SPR最让人头疼的是光源功率抖一下信号就跟着跳一下十有八九会把环境波动误判成样品折射率变化。双波长强度调制把这个问题变成了数学题同一个ASE光源分出两个波长两个探测通道同时受光源波动影响做减法之后共模项被干掉。Kretschmann棱镜耦合的四层结构下用MATLAB建一个反射率模型就能把这些现象完整还原。这里没有新器件只有对传统强度调制的轻微修改但实测效果是光源稳定性不再那么关键偏振镜也可以省掉。下面我把整个仿真流程从四层模型、材料参数、差分计算到入射角扩展逐步拆开光学传感器方向的从业者和研究生可以直接照着复现。2. Kretschmann四层模型中真正决定反射率的是介电常数的波长色散2.1 四层结构与Fresnel反射系数的递推关系SPR反射率不是普通镜面反射而是p偏振光的倏逝波与金膜表面自由电子集体振荡相互耦合的结果。Kretschmann结构里光先经过棱镜介质在棱镜/金属界面发生全内反射当入射角的横向波矢与表面等离子体波矢匹配能量被耦合进SPW反射率曲线出现一个明显的“共振浸没”。待测介质折射率一变谷底位置就跟着变这就是传感器响应的来源。建模时把结构分成四层棱镜、铬膜、金膜、待测介质。铬膜起粘附作用厚度一般只有2nm金膜是激发SPW的核心层通常取50nm。每一层都有自己的介电常数金属层必须用复数介电常数虚部对应吸收损耗。按等效界面法从最底层往上递推p偏振的反射系数可以写成r1234 (r12 r234 * exp(2i * kz2 * d2)) / (1 r12 * r234 * exp(2i * kz2 * d2))其中r12、r23、r34是两个相邻介质的Fresnel反射系数kz是各层波矢在z方向的分量d2、d3分别是铬层和金层厚度。这个表达式看起来简单但实际计算时每一层的介电常数都要随波长变化尤其金的折射率在1550nm附近大约为0.59.8i的量级虚部很大。虚部如果取错共振谷会变得又浅又宽甚至完全看不到吸收峰。下面是四层结构的主要参数仿真时建议在代码里写成可配置的结构体不要散落成一堆魔法数字层号材料厚度光学参数模型说明1K9棱镜无限厚n 1.5163近红外可查色散公式半圆柱棱镜入射光从棱镜面进入2铬2 nm复数折射率Palik数据提高金膜附着力3金50 nm复数折射率Palik数据激发表面等离子体波4待测介质半无限折射率约1 ~ 1.0008常见气体折射率接近真空如果考虑K9玻璃的色散可以使用论文中的多项式形式ε a0 a1λ² a2λ⁻² a3λ⁻⁴ a4λ⁻⁶ a5λ⁻⁸λ单位是微米。但1540nm和1560nm相差仅20nmK9折射率的色散变化在小数点后第四位量级实际仿真中我经常先取常数1.5163等整体趋势跑通后再把色散公式加进去。金属层的色散不能省因为复数折射率虚部决定SPR谷底的深度和宽度。2.2 MATLAB实现反射率函数我习惯把反射率计算封装成独立函数输入角度、波长、材料和厚度输出p偏振反射率。这样后面做双波长差分和角度扫描时只需要循环调用同一个函数。function R spr_refl_p(theta_deg, lambda_nm, n_prism, n_cr, n_au, d_cr_nm, d_au_nm, n_analyte) % theta_deg: 入射角单位度 % lambda_nm: 入射光波长单位nm % n_prism: 棱镜折射率 % n_cr: 铬膜复数折射率 % n_au: 金膜复数折射率 % d_cr_nm: 铬膜厚度nm % d_au_nm: 金膜厚度nm % n_analyte: 待测介质折射率 k0 2 * pi / (lambda_nm * 1e-9); % 真空波矢单位1/m theta theta_deg * pi / 180; beta k0 * n_prism * sin(theta); % 棱镜内横向波矢 eps1 n_prism^2; eps2 n_cr^2; eps3 n_au^2; eps4 n_analyte^2; kz1 sqrt(k0^2 * eps1 - beta^2); kz2 sqrt(k0^2 * eps2 - beta^2); kz3 sqrt(k0^2 * eps3 - beta^2); kz4 sqrt(k0^2 * eps4 - beta^2); r34 (kz3/eps3 - kz4/eps4) / (kz3/eps3 kz4/eps4); r23 (kz2/eps2 - kz3/eps3) / (kz2/eps2 kz3/eps3); r12 (kz1/eps1 - kz2/eps2) / (kz1/eps1 kz2/eps2); d2 d_cr_nm * 1e-9; d3 d_au_nm * 1e-9; r234 (r23 r34 * exp(2i * kz3 * d3)) / (1 r23 * r34 * exp(2i * kz3 * d3)); R abs((r12 r234 * exp(2i * kz2 * d2)) / (1 r12 * r234 * exp(2i * kz2 * d2)))^2; end代码里最关键的是kz sqrt(k0^2 * eps - beta^2)这一行。当入射角超过全内反射角时kz会变成虚数exp(2i * kz * d)变成实数指数衰减这正是倏逝波耦合进金属层的物理过程。如果角度选得太小kz是实数模型退化成普通薄膜干涉看不到SPR谷底。调用时金属折射率必须给复数。1550nm附近我用的参考值是n_cr 3.2 3.1i、n_au 0.55 9.8i正式仿真建议用Palik表在目标波长处插值。同一套结构在不同波长下金膜折射率虚部差异会直接影响谷底深度。先用常数跑通流程再换成插值表是排查模型问题最快的路径。2.3 角度扫描找到共振谷底拿到反射率函数后第一步不是直接算双波长而是扫角度确定当前折射率下的共振角度。比如固定1550nm扫40到45度观察反射率最小值对应的角度。这一步很关键因为后续双波长差分需要在共振角度附近工作角度偏了反射率对折射率的变化率会明显下降。theta_scan 40:0.001:45; R_scan zeros(size(theta_scan)); for i 1:numel(theta_scan) R_scan(i) spr_refl_p(theta_scan(i), 1550, 1.5163, 3.23.1i, 0.559.8i, 2, 50, 1.0003); end [RMIN, idx] min(R_scan); theta_min theta_scan(idx); fprintf(共振角度: %.3f deg, 反射率最小值: %.4f\n, theta_min, RMIN);这段代码输出的是当前折射率下的共振角度。如果金属折射率取错theta_min会偏移好几度或者RMIN不够低。我一般会同时打印RMIN如果反射率最小值高于0.01就先查金属虚部是否没写对而不是急着调角度。薄层厚度同样影响谷底深度2nm铬层在结构里不是可有可无它会稍微压低反射率也会让共振位置偏移零点零几度。3. 双波长强度调制用差分信号把光源漂移“减”掉3.1 为什么单波长强度调制那么脆传统强度调制SPR的做法是固定一个入射角和波长用反射光强作为待测折射率的指示。问题在于反射光强等于入射光强乘以反射率光源功率只要波动1%信号看起来就像折射率变了0.0001甚至更多。高稳定激光器能解决一部分问题但仪器成本和体积都会变大远程长时间在线监测时仍会被环境温度、光纤弯曲损耗拖累。双波长方案不追求光源绝对稳定而是让两个波长共用同一个ASE光源通过光纤滤波器分别选出1540nm和1560nm。两个波长的反射光分别用探测器接收用两个反射率的差作为输出。光源功率波动对两路信号的影响是同步的做差之后这个共模干扰被消掉。用数学表达就是DeltaR I1 / I0 - I2 / I0 (I1 - I2) / I0这里I0是入射光强度I1和I2是两个波长各自反射后的强度。因为两个波长来自同一个光源I0波动时I1和I2以相同比例变化做差后波动项自然抵消。这比后端做数字滤波更直接属于信号链路上的共模抑制。3.2 波长间隔为什么选20nm两个波长不能离得太远也不能太近。如果间隔太大金膜在两侧波长的折射率实部差异明显SPR响应曲线形态不一样差分信号和折射率之间的线性关系很快被破坏。间隔太小比如5nm两个波长在共振区附近的反射率变化几乎一样差分信号幅度太小抗噪优势体现不出来。论文中选1540nm和1560nm差距20nm正好落在ASE光源平缓输出的区域内两个波长的光强能量接近探测器量程也容易配平。在仿真时需要同时计算两个波长的反射率差。下面这段代码复现论文中的核心过程lambda1 1540; lambda2 1560; theta0 42.08; n_range linspace(1, 1.0008, 81); dR zeros(size(n_range)); R1_all zeros(size(n_range)); R2_all zeros(size(n_range)); for i 1:numel(n_range) R1_all(i) spr_refl_p(theta0, lambda1, 1.5163, 3.23.1i, 0.559.8i, 2, 50, n_range(i)); R2_all(i) spr_refl_p(theta0, lambda2, 1.5163, 3.23.1i, 0.559.8i, 2, 50, n_range(i)); dR(i) R1_all(i) - R2_all(i); end这里的n_range覆盖气体折射率1到1.0008步长0.00001。dR计算的是两个波长反射率的算术差物理上对应两路探测器信号归一化到入射光强后的差值。需要关注的是dR曲线在哪个区间内接近直线只有线性区间内的数据才能用一次拟合公式反推折射率。3.3 灵敏度计算与单位陷阱灵敏度定义是传感器输出变化与待测折射率变化的比值S DeltaR / Delta_n单位写作% / RIURIU是折射率单位。论文中的仿真结果是S 28582 %/RIU这个数值看起来很大要把它换算成直观感受。假设折射率变化0.0001那么反射率差变化约为2.86%。这个量级对探测器来说非常容易分辨也是双波长差分方案灵敏度可用的原因。用MATLAB做线性拟合时要小心百分比换算idx n_range 1.0004 n_range 1.0005; p polyfit(n_range(idx), dR(idx), 1); S p(1) * 100; % dR是0~1之间的小数乘以100转换为% fprintf(线性灵敏度: %.0f %% / RIU\n, S);polyfit返回的斜率是反射率差随折射率的变化率单位是RIU的倒数。因为反射率本身在0到1之间乘以100后变成百分比才对得上论文里的口径。如果忘记乘100数值会变成285.82 %/RIU后面做标定时会差出100倍。3.4 线性区间不是全程直线从图4可以看出折射率在1到1.0008之间时dR整体是一条带弯曲的曲线只有1.0004到1.0005这段接近直线。这是因为双波长的差分信号本质上是两个SPR反射曲线的“斜率差异”靠近共振谷底时斜率变化快远离时变化平缓。实际使用中必须限定测量范围或者用更高阶拟合。n_range的步长也要选得足够细。我用81个点步长0.00001每个点计算两次反射率耗时在毫秒级。步长太粗比如0.0001拟合出的线性区间会被空过去结果不够平滑。步长太细后面做多角度扫描时循环次数增多MATLAB也不会慢但没必要。4. 测试气体折射率范围1.0008入射角扫描把测量范围“拼接”出来4.1 初始入射角42.08°是怎么确定的双波长仿真必须有一个确定的工作角度。42.08°不是拍脑袋来的而是对1.0003左右折射率做角度扫描找到反射率谷底最深的那个角度。金膜厚度50nm、波长1550nm附近共振角就在这个位置。入射角偏差0.01度反射率谷底位置就会偏移双波长差分曲线也会变形。在仿真中我通常先做一个“角度-折射率”二维扫描把共振角度随折射率变化的关系整体画出来。流程是用上述spr_refl_p函数对每一个n_analyte扫描角度记录谷底位置theta_candidates 42.0:0.002:42.3; n_test [1, 1.0003, 1.0006, 1.0008]; for n n_test R_min 1; theta_min 0; for th theta_candidates R spr_refl_p(th, 1550, 1.5163, 3.23.1i, 0.559.8i, 2, 50, n); if R R_min R_min R; theta_min th; end end fprintf(n%.4f, 谷底角%.3f, Rmin%.4f\n, n, theta_min, R_min); end这样能提前知道不同折射率对应的共振角度变化范围。折射率从1增大到1.0008共振角会往大角度方向移动移动量大约是零点零几度。论文里用42.07到42.091的四个角度覆盖这个区间正是基于这种单调关系。4.2 线性区间与灵敏度结果固定入射角42.08°双波长1540nm与1560nm的dR曲线在1.0004到1.0005之间表现出很好的线性。这个区间对应共振谷底两侧最陡峭的位置反射率差对折射率变化最敏感。超出这个区间曲线弯曲加剧再用直线拟合会产生明显误差。用polyfit在1.0004到1.0005区间拟合得到的灵敏度约28582 %/RIU。这里需要解释一下“测量范围扩大”的具体含义。单波长强度调制通常在固定角度下只能覆盖一段很窄的折射率范围双波长差分虽然去掉了共模噪声但线性区间依然有限。论文的解决办法是改入射角让响应曲线在折射率轴上“平移”每个角度负责一段范围最后把1到1.0008整体覆盖。4.3 多入射角拼接的迭代思路对每一个新的入射角重复双波长差分计算只截取该角度下线性度最好的折射率区间。相邻角度的有效区间要留一点重叠防止接缝处出现断点。下面这段代码演示了如何批量生成不同入射角的dR曲线theta_list [42.070, 42.077, 42.084, 42.091]; figure; hold on; for t theta_list dR_t zeros(size(n_range)); for i 1:numel(n_range) R1 spr_refl_p(t, lambda1, 1.5163, 3.23.1i, 0.559.8i, 2, 50, n_range(i)); R2 spr_refl_p(t, lambda2, 1.5163, 3.23.1i, 0.559.8i, 2, 50, n_range(i)); dR_t(i) R1 - R2; end plot(n_range, dR_t * 100, LineWidth, 1.2); end绘制后每条曲线的零点和斜率都有差异。实际使用时对每个角度保存一个有效折射率区间和一个线性拟合系数入射角有效折射率范围灵敏度备注42.070°1.0000 ~ 1.0002略低靠近折射率下限42.077°1.0002 ~ 1.0005中与前后搭接42.084°1.0004 ~ 1.0007较高论文中灵敏度主要来自这个区间附近42.091°1.0006 ~ 1.0008中高负责上限范围这些区间是示意性的具体边界取决于你用的金属折射率数据。拿到四个角度的标定表后检测未知气体时先粗测一次dR值根据落在哪个量程选择对应角度再精测。这种查表角度切换的方法比单角度拟合更实用。4.4 入射角扫描时要注意的细节入射角步长不能太大。42.070和42.077之间只差0.007度但有效区间已经能移动0.0002左右RIU。如果用0.1度步长曲线会跳得厉害中间出现空白区。另外每个角度的有效区间不是等宽的灵敏度高的地方区间窄灵敏度低的地方区间宽不能平均分配。如果某个角度下dR曲线出现明显振荡先检查kz在对应折射率下是否接近零。当横向波矢与某一层的传播常数接近时数值上容易出现奇异性表现为反射率突然跳到1以上或者变为负数。这时需要把角度扫描步长变小或者改用双精度计算MATLAB默认双精度足够但要注意sqrt函数内部对负数取实部的问题。5. 仿真排错与实验衔接三个容易翻车的地方5.1 金属复折射率必须用复数否则谷底浅如平地仿真中最常见的“仿真发散”其实是物理参数错误。金膜折射率虚部如果被当成0SPR共振谷就不会出现反射率曲线几乎是一条平线。自检方式是计算无吸收介质对照把金膜折射率虚部设成0看反射率是否接近1再恢复虚部看谷底是否出现。如果两个状态下差异极小说明金属层参数没有正确传递进函数。R_no_loss spr_refl_p(42.08, 1550, 1.5163, 3.23.1i, 0.550i, 2, 50, 1.0003); R_with_loss spr_refl_p(42.08, 1550, 1.5163, 3.23.1i, 0.559.8i, 2, 50, 1.0003); fprintf(无吸收谷底: %.4f, 有吸收谷底: %.4f\n, R_no_loss, R_with_loss);如果R_no_loss很小多半是角度扫描范围没找对或者棱镜折射率与波矢方向不匹配。这里推荐先固定折射率1.0003扫描40到45度确认谷底角度后再展开双波长计算。5.2 双波长间隔和折射率扫描范围要匹配1540nm和1560nm的组合在1到1.0008范围内给出了可用结果。如果把波长间隔放大到50nm两个波长对应的金膜折射率差异变大dR曲线可能在区间内出现二次弯曲线性拟合失效。缩小到10nm虽然更线性但信号幅度变小实验上信噪比下降。20nm是平衡点。如果你的样品折射率范围更大比如液体1.33到1.36固定20nm间隔往往不够需要重新做波长选择扫描。我一般会先计算不同波长对在目标折射率区间内的最大线性偏差选偏差最小的一组。仿真发散时不要只盯角度先缩小波长间隔试一次如果发散消失说明色散差异是主因。5.3 实验装调时先把两个探测通道归一化仿真中的I0是理想常数但实验上两个光纤滤波器的透过率不可能完全一致两个探测器的响应度也不同。装调时先用一个已知反射率的标准样品比如空气把两路信号调整到同一个基线。归一化后再进行差分测量否则I1 - I2里会固定叠加一个系统偏置表现为折射率读数整体偏移。偏振镜在双波长方案中不再是必需品因为s偏振光不会激发SPW它的反射率对折射率不敏感做差后两路的s偏振贡献相等自然抵消。实际搭建时可以把偏振镜放在光路中做验证旋转偏振镜观察两路信号是否同时变化且差值几乎不变以此确认光路对准是否正确。最后建议把所有仿真参数集中写成一个配置文件角度、波长、金属厚度、折射率数据来源都做成变量。做参数扫描时只改配置文件不要边写边改主函数否则很可能因为某个角落的魔法数字导致结果对不上。先用1400nm到1600nm的波长范围扫一遍确认共振谷位置随折射率变化的趋势稳定再收敛到1540/1560nm做双波长差分这是最省时间的调试顺序。本文还有配套的精品资源点击获取