Matlab实现法诺共振拟合与Q因子提取全流程
2026/9/8 6:13:21 网站建设 项目流程

上周我在整理介质超表面样品的反射谱时,又撞见那个老问题:共振峰明显不对称,左边陡得几乎是直线掉下去,右边却拖出一条平缓的长尾巴。用最常见的洛伦兹线型去拟合,残差永远呈“S”形弯在零轴两侧,怎么调都压不平。干我们这行的一看就懂——这是典型的法诺共振(Fano resonance),离散亮态和连续暗态干涉产生的非对称线型。想把这组谱线里的共振位置和线宽准确提取出来,进而计算Q因子,光凭眼睛读峰位肯定不行,数据一多也容易出错。我最终用Matlab搭了一套法诺共振拟合流程,从谱线导入到Q因子输出,整套脚本固化下来。这篇博文就是把方法原原本本写清楚,公式怎么理解、初值怎么定、用哪个函数最稳、算Q因子有哪些想不到的坑,一步不落。

1. 法诺共振到底在拟合什么:公式拆解与物理含义

1.1 Fano线型的数学表达:一个公式里的非对称来源

先看最通用的法诺共振拟合公式,我习惯写成:

I(x) = y0 + A · [(q + ε)² / (1 + ε²)]

其中 ε = 2(x - Er) / Γ

这五个参数各有各的物理身份:

  • x 是横轴,可以是光子能量(eV)、频率(THz)或者波长(nm),取决于你的实验设备和数据单位;
  • Er 是共振中心位置,也就是离散态的本征能量或频率,拟合出来的第一个核心结果;
  • Γ 是共振线宽,严格说就是半高全宽(FWHM),它反映了共振态寿命和损耗大小,第二个核心结果;
  • q 是法诺参数,描述线型不对称程度,q 越大线型越接近对称的洛伦兹峰,q 越接近 0,谷越深、非对称越明显;
  • A 是振幅系数,控制整体强度大小,可正可负,正负号决定了谱线是峰还是谷;
  • y0 是背景平移项,用来吸收测量中的直流偏移和连续背景。

那这个公式为什么能产生非对称形状?关键在分子的 (q + ε)²。当 ε 从负到正变化时,(q + ε)² 会在 ε = -q 处取到 0,也就是说谱线会有一个强度为零的点。与此同时分母 1 + ε² 又在 ε = 0 附近限制着整体强度,两者叠加,就会在共振位置一侧形成尖锐的谷、另一侧形成缓和的肩,或者反过来。这个“谷”并不是仪器的噪声坑,而是连续态和离散态相消干涉的真实物理结果。

这里我想特别提示一个容易被新手忽略的点:法诺公式在不同论文里有好几种等价写法,有些版本会在括号后面减 1,写成 A[(q + ε)²/(1 + ε²) - 1],有些版本会乘以额外常数。不同写法会让拟合出来的 A 和 y0 数值完全不同,但 Er、Γ、q 这几个关键参数是稳定的。所以我建议固定使用不加“-1”的版本,这样和大多数文献、开源代码交流起来最省事。

1.2 Q因子不是照抄谱线宽度:定义、换算与物理意义

Q因子,全称品质因子(Quality factor),在共振现象里是衡量共振“尖锐程度”的核心指标,通用定义是:

Q = Er / Γ

分子是共振中心位置,分母是线宽。这里的核心前提是:Er 和 Γ 的单位要一致。也就是说,如果拟合时横轴用的是电子伏特(eV),那么 Er 和 Γ 都用 eV,除出来的 Q 是一个无量纲数;如果横轴是纳米(nm),就用 nm 单位下的中心波长和线宽相除。很多同学第一次算 Q 时容易犯的错,是把中心用 eV、线宽用纳米直接混着除,结果错误离谱自己还没发现。

Q 因子的物理意义其实很直接:它表示共振系统在振荡一个周期内存储的能量与损耗能量的比值。Q 越高,说明系统损耗越低,能量在腔体或结构里“存活”的时间越长。换算成光子寿命的话:

τ = 2Q / ω0

其中 ω0 = 2πf0 是共振角频率。这个关系式在做时间分辨测量或者估算慢光效应时经常用到。高 Q 的微纳结构,光子寿命可达皮秒甚至纳秒量级,对应极窄的线宽,这也是为什么超表面连续域束缚态(BIC)结构总能刷新高 Q 记录的原因。

需要多说一句:如果拟合出来的法诺线型不对称性很强,你直接用肉眼在原始数据上量“半高宽”,得到的数值会和拟合出的 Γ 有明显出入,因为非对称线型的视觉宽度和理论宽度定义并不完全一致。这也是为什么要相信拟合参数而不是目测取值。

1.3 什么时候必须用法诺拟合,什么时候洛伦兹就够了

在实际工作中,我一般先做一次洛伦兹拟合作为预判。如果残差在共振位置两侧呈现明显的“S”形波动,那基本可以判定需要切换到法诺模型。反过来,如果你的谱线虽然不对称,但洛伦兹拟合的残差已经接近噪声水平,那说明不对称性弱到可以忽略,用洛伦兹并不影响后续 Q 值提取。

还有一个经验性的判断方法:看谱线的谷底是否接近零强度。法诺共振的相消干涉会让谷底压得特别低,甚至趋近于零(q 接近 0 时),而单纯的洛伦兹谷通常不会低到那个程度。另外,法诺线型的一个显著特征是不对称峰的两侧极值位置并不对称地分布在中心两侧,一翼平缓一翼陡峭,这个肉眼就能分辨。

我的建议是:只要谱线肉眼可见地不对称,直接上法诺模型。因为法诺模型包含洛伦兹极限(q 很大时),你拟合出来的 q 如果非常大,自然就退化成洛伦兹了。多两个参数换来的是更普适的模型,代价只是初值要稍微花点心思。

2. 拟合方案选型与Matlab工具箱准备

2.1 为什么是Matlab:批处理、自定义模型和工程衔接

每次有学生问我“能不能用 Origin 做”,我都会说能,但很痛苦。Origin 的全局拟合功能虽然也能自定义函数,但你要处理几十条谱线时,每一条都要手动设置初值、手动导出结果,这种重复劳动太消耗精力了。而 Matlab 的优势在于脚本化:定义好模型函数、写一个循环,几十条谱线丢进去,几分钟后一张参数表格就出来了,还能自动生成全部拟合对比图。

另一个很现实的原因是团队协作。我们实验室的大部分光学模型、时域有限差分仿真数据、甚至设备控制程序都用 Matlab 写,谱线拟合和后续数据处理在同一套环境下完成,省去了跨语言搬运数据的麻烦。Matlab 的脚本可读性好,后来接手的人也好维护。如果你是全 Python 栈,用 scipy 的 curve_fit 也能做类似的事,但在这篇博文里我以 Matlab 为主线讲全套思路。

2.2 优化工具箱与统计工具箱的版本差异

法诺拟合主要依赖两个工具箱:Optimization Toolbox 提供核心的 lsqcurvefit 非线性最小二乘求解器,Statistics and Machine Learning Toolbox 提供 fitnlm 这个更顺手的非线性回归接口,以及 nlparci 参数区间估计函数。

上机之前先确认工具箱是否安装到位,两条命令:

ver('optim') ver('stats')

如果输出里显示版本信息就说明没问题。也可以用 license 检测:

license('test', 'Optimization_Toolbox') license('test', 'Statistics_Toolbox')

返回 1 表示有许可证可用。

这里有几个版本相关的坑值得提一句:早期版本(R2015a 之前)用的是 optimset,新版用 optimoptions,老语法在新版里会报错或警告。另外,lsqcurvefit 和 fitnlm 在高版本 Matlab 里的输出结构基本稳定,但如果你的代码要给别人在旧版本上跑,最好加一段版本判断或者统一使用更稳定的 fitnlm,它的接口变化相对小一些。至于评论里很多人提到的工具箱下载、安装包之类的问题,我建议直接用学校或公司提供的正版授权,省去一堆环境兼容的麻烦。

2.3 核心算法选择:lsqcurvefit和fitnlm的取舍

法诺拟合本质上是一个五参数非线性最小二乘问题,数据点几千个、参数只有五个,问题规模不大,Levenberg-Marquardt(LM)算法就够用。

两个函数怎么选,我给出自己的判断:

  • 如果你只需要快速得到参数,不太关心置信区间,那就用 lsqcurvefit。它对自定义模型函数的定义方式最直接,输出自由度大,还能设置参数边界、添加额外约束。
  • 如果你想拿到参数标准误、置信区间、残差统计量,写论文时直接引用误差棒,那 fitnlm 更省事,它内部帮你做了误差传播和线代运算。
  • 如果初值给不准、担心陷入局部极小值,先用 particleswarm 或 GlobalSearch 做一轮全局预搜索,再用结果作为初值交给 lsqcurvefit 精修。

我自己的主力方案是:用 lsqcurvefit 做核心拟合,配合 nlparci 算 95% 置信区间。原因很简单,lsqcurvefit 的边界约束写起来特别灵活,可以针对物理上不可能的参数组合直接封死,比如把 Γ 限制在大于 0 的区间,防止拟合优化器跑到负线宽这种荒谬结果。fitnlm 虽然也有参数界,但处理起来不如 lsqcurvefit 那么顺手。

3. 从数据到Q因子:法诺拟合完整流程演示

3.1 构造带噪声的模拟法诺谱,先跑通脚本

这里我用一段模拟数据来演示,最大好处是“真值已知”:我可以验证拟合算法能不能恢复出预先设定的参数,确认整个流程没有 bug 之后再套用到实验数据上。

先定义法诺模型函数,在 Matlab 里新建一个 fano_model.m 文件:

function y = fano_model(x, p) % 法诺共振模型 % p(1): Er 共振中心 % p(2): Gamma 线宽(FWHM) % p(3): q 法诺参数 % p(4): A 振幅 % p(5): y0 背景平移 Er = p(1); Gamma = p(2); q = p(3); A = p(4); y0 = p(5); f = (x - Er) ./ (Gamma / 2); y = y0 + A .* (q + f).^2 ./ (1 + f.^2); end

然后在主脚本里构造带噪声的数据:

% 构造模拟数据 rng(1); % 固定随机种子,保证结果可复现 x = linspace(0.6, 1.8, 1000)'; % 光子能量,单位 eV p_true = [1.25, 0.08, -2.5, -1.5, 0.9]; % 真值 y_clean = fano_model(x, p_true); y_noisy = y_clean + 0.03 * randn(size(x)); figure('Color', 'w'); plot(x, y_noisy, '.', 'MarkerSize', 5); hold on; plot(x, y_clean, 'k-', 'LineWidth', 1.2); xlabel('光子能量 (eV)'); ylabel('反射率 (a.u.)'); legend('含噪声数据', '真实法诺线型', 'Location', 'southeast');

运行之后你能在图上看到一个典型的非对称谷:左侧相对光滑地下降,右侧收尾更快或者更慢。为什么加噪声?因为实验数据永远伴随噪声,如果不先在模拟数据里加上噪声检验拟合算法的抗噪能力,直接上实验数据会很容易被各种异常结果打得措手不及。

3.2 初始参数估计:五维参数逐个定位,不靠猜

非线性拟合最怕的就是初值乱给。五个参数你不可能全凭运气。我的经验是按顺序逐个估计:

第一步,Er 的初值。对于法诺共振,谱线极值点不完全等于 Er,尤其是在 q 绝对值接近 2 这种中等不对称情况下,半峰点和极值点都会有偏移。所以我更推荐用谱线重心法来估:

y_shifted = y_noisy - min(y_noisy); Er0 = sum(x .* y_shifted) / sum(y_shifted);

这相当于把谱线当作质量分布求质心,作为 Er 的初值即使不完美,也差不到哪去。

第二步,Γ 的初值。在谱图上找到谷底和肩部极值之间的横轴距离,估算一个值。比如你看到谷在 1.25 eV 附近,肩峰在 1.32 eV 附近,距离是 0.07 eV,那么 Γ 的初值可以取 0.06 到 0.1 之间。

第三步,q 的初值。看谷的深度和方向:如果右侧平缓左侧陡峭,q 往往是正值;反过来是负值。谷越深、q 绝对值越小,通常 q 在 -1 到 -5 区间;如果线型几乎对称,q 绝对值可能在 10 以上。

第四步,A 和 y0。y0 取远离共振位置的背景平均值,A 取谱线最低点与 y0 的差。如果你把数据归一化过,这两个参数会小一些。

综合起来,对这条模拟数据,我大致给的初值是:

p0 = [1.3, 0.06, -1.0, -1.0, 0.95];

虽然不是特别准,但离真值已经比较近,足以让 LM 算法收敛。

3.3 约束非线性拟合、参数输出与残差诊断

接下来调用 lsqcurvefit 执行拟合。这里我强烈建议加上边界约束,理由很简单:把 Γ 限制在正区间可以避免优化器为了压低残差跑出负线宽;把 q 限制在一个合理范围内可以避免它跑到几百上千的退化情况,那种时候你根本没法解释拟合结果。

% 定义边界,单位都和横轴一致 lb = [1.0, 0.005, -20, -10, -0.5]; ub = [1.5, 0.300, 20, 10, 2.0]; opts = optimoptions('lsqcurvefit', ... 'Display', 'final', ... 'MaxFunctionEvaluations', 1e4, ... 'MaxIterations', 2000, ... 'FunctionTolerance', 1e-10, ... 'StepTolerance', 1e-10); [p_fit, resnorm, residual, exitflag] = lsqcurvefit(@fano_model, p0, x, y_noisy, lb, ub, opts);

拟合完成后,立刻画对比图和残差图:

x_fit = linspace(x(1), x(end), 2000)'; y_fit = fano_model(x_fit, p_fit); figure('Color', 'w'); plot(x, y_noisy, '.', 'MarkerSize', 5); hold on; plot(x_fit, y_fit, 'r-', 'LineWidth', 1.6); xlabel('光子能量 (eV)'); ylabel('反射率 (a.u.)'); legend('实验数据', '法诺拟合', 'Location', 'southeast'); title(sprintf('Er=%.4f, \\Gamma=%.4f, q=%.2f', p_fit(1), p_fit(2), p_fit(3))); figure('Color', 'w'); plot(x, residual, 'k.'); xlabel('光子能量 (eV)'); ylabel('拟合残差');

判断拟合质量有三个标准:第一看残差是否围绕零轴均匀分布,如果在某个区域系统性偏高、某个区域系统性偏低,说明模型形式或拟合区间有问题;第二看参数是否全部落在你预设的物理区间内;第三看拟合线在数据密集区是否穿过最多数据点。我通常还会顺手算一个决定系数 R²:

SS_res = sum(residual.^2); SS_tot = sum((y_noisy - mean(y_noisy)).^2); R2 = 1 - SS_res / SS_tot;

R² 在 0.99 以上基本是很漂亮的拟合了,但要注意一点:法诺模型本身自由度多,R² 高不代表参数一定准确,必须结合残差形态一起判断。

3.4 计算Q因子及不确定度:一次拟合给完整结论

拟合收敛后,Q 因子就很简单了:

Er_fit = p_fit(1); Gamma_fit = p_fit(2); Q = Er_fit / Gamma_fit;

如果横轴是 eV,那么 Q 就是无量纲数。为了在论文里写误差棒,我用 nlparci 提取 95% 置信区间。lsqcurvefit 的第七个输出是 Jacobian 矩阵,可以直接喂给 nlparci:

[~, ~, ~, ~, ~, ~, J] = lsqcurvefit(@fano_model, p_fit, x, y_noisy); ci = nlparci(p_fit, residual, 'jacobian', J); % 提取 Er 和 Gamma 的标准误 sigma_Er = (ci(1,2) - ci(1,1)) / (2 * 1.96); sigma_Gamma = (ci(2,2) - ci(2,1)) / (2 * 1.96); % 误差传播计算 Q 的标准误 sigma_Q = Q * sqrt((sigma_Er / Er_fit)^2 + (sigma_Gamma / Gamma_fit)^2);

这个误差传播公式本质上是把 Q 看作 Er 和 Γ 的比值,对两个变量分别求偏导后合起来的近似估计。当两个参数的标准误都远小于参数本身时,这个公式非常可靠。如果你的 Γ 拟合误差超过 10%,就得小心了,说明拟合可能不稳定,需要回头检查初值和数据质量。

如果嫌 lsqcurvefit 输出置信区间麻烦,可以直接用 fitnlm:

fano_fun = @(b, x) b(5) + b(4) .* (b(3) + (x - b(1)) ./ (b(2)/2)).^2 ./ (1 + ((x - b(1)) ./ (b(2)/2)).^2); mdl = fitnlm(x, y_noisy, fano_fun, p0); Er_est = mdl.Coefficients.Estimate(1); Gamma_est = mdl.Coefficients.Estimate(2); Q_est = Er_est / Gamma_est; SE_Er = mdl.Coefficients.SE(1); SE_Gamma = mdl.Coefficients.SE(2); SE_Q = Q_est * sqrt((SE_Er / Er_est)^2 + (SE_Gamma / Gamma_est)^2);

两种方法的结果基本一致,fitnlm 还能直接给出残差分析的多个统计量,适合用于正式的统计分析。我平时为了统一流程,主要用 lsqcurvefit,但需要快速出标准误时会切到 fitnlm。

4. 拟合结果的物理解读与场景扩展

4.1 不同微纳结构里的典型Q值范围

把拟合流程跑通之后,对拟合出的 Q 值有一个合理的物理区间预判,能帮你判断结果是否正常。我根据自己做过的结构和看到的文献,整理了下面这个典型范围表:

结构类型典型Q值范围线型特征常见场景
金属等离激元纳米颗粒(偶极共振)5~30接近对称洛伦兹局域表面等离激元传感
等离激元法诺纳米结构(如七聚体/圆盘环)10~80明显非对称法诺折射率传感、表面增强光谱
介电超表面准连续域束缚态(quasi-BIC)100~10^4线宽窄、可对称可非对称滤波、激光、非线性光学
光子晶体微腔10^3~10^6多为对称洛伦兹腔量子电动力学、慢光
分子/激子与等离激元强耦合100~500可表现法诺特征室温量子光学、极化激元

这张表可以当作拟合结果的“体检表”。比如你在介电超表面上测到 Q = 3,那基本说明哪里出了问题,要么数据噪声太大,要么拟合区间选取不当,要么结构本身损耗确实很大需要重新审视。

4.2 法诺参数与耦合、损耗、对称性的关系

拟合得到的 q 值并不仅仅是一个波形参数,它背后有明确的物理含义。q 反映的是离散态与连续态之间的相位关系和耦合强度:|q| 接近 0 时相消干涉最强,谷最深,非对称性最明显;|q| 很大时连续态贡献可以忽略,线型退化为对称的洛伦兹峰,这通常意味着离散态很弱地耦合到连续态背景,或者连续态背景本身很弱。

我在做超表面结构参数扫描时经常发现一个规律:当结构的对称性被微小破坏时(比如圆盘变椭圆),q 值会从几百急剧降到几个,同时线宽明显变窄、Q 值大幅上升。这是准连续域束缚态的典型行为——对称性破缺程度直接控制着辐射损耗,进而控制法诺线型和 Q 值。所以拟合出的 q 值可以反过来指导你判断结构的对称性状态,这在做样品质量控制时特别好用。

另一个常用场景是环境折射率变化。当待测介质折射率变化时,Er 会线性漂移,而 Γ 一般来说变化不大,于是 Q 值基本保持不变。这就是为什么共振传感领域更关注的是波长移动量和谱线宽度的比值,而不是单单看 Q 值。

4.3 从拟合参数到传感灵敏度与光子寿命

拟合参数除了算 Q,还能延伸出几个很实用的物理量。最直接的是光子寿命 τ = 2Q/ω0。比如你在 1.25 eV 处拟合出 Q = 15,那么 ω0 ≈ 1.90 × 10^15 rad/s,τ ≈ 15.8 fs。这个时间尺度在金属等离激元结构里很常见,对于超表面准连续域束缚态,Q 可以到 10^3 以上,光子寿命就进入皮秒量级。

对于传感应用,还可以计算品质因数 FOM:

FOM = S · Q / λ0

其中 S 是折射率灵敏度(nm/RIU),λ0 是中心波长。举个例子,某个结构 S = 400 nm/RIU,Q = 50,中心波长 800 nm,那么 FOM = 25。在很多传感文献里,FOM 比 Q 值更能体现一个结构的实际探测能力,因为它同时考虑了共振位置的漂移量和谱线本身的尖锐程度。

这些延伸计算都不需要额外测量,直接从拟合参数里就能推算出来。这也是为什么我强烈建议把拟合脚本封装成函数的原因:参数一出来,Q、τ、FOM 全都自动算好,省去大量重复劳动。

5. 经验复盘:拟合中常见问题和排查技巧

5.1 参数发散或掉进局部极小值

这是最常遇到的问题,表现为:拟合结果显示 q 跑到几百甚至几千,或者 Γ 接近边界值,或者 Er 明显偏离谱线中心。这类问题的根源基本都在初值或边界上。

我的排查顺序是:先把数据缩小到共振附近区域重新拟合,排除远端背景点对拟合的拉扯;然后检查初值是否离真值太远,尤其是 Er,建议用重心法重新估一遍;最后给参数加上合理的边界约束。如果边界约束加上还是发散,那多半是模型不适合这段数据,比如谱线本身有多重共振,单法诺模型描述不了。

如果手上有足够时间,我会用全局优化做一次预搜索:

rng(0); lb_g = [0.8, 0.001, -50, -20, -1]; ub_g = [1.8, 0.500, 50, 20, 3]; p_g = particleswarm(@(p) sum((y_noisy - fano_model(x, p)).^2), 5, lb_g, ub_g); p0 = p_g;

particleswarm 虽然慢一点,但能帮你跳出初值的困境,之后再用 lsqcurvefit 精修一波,成功率非常高。

5.2 单位换算、基线漂移和背景处理

Q 因子的计算严重依赖单位一致。用 eV 拟合就用 eV 算 Q,用 nm 拟合就用 nm 算 Q。如果横轴本来是波长(nm),但你读到另一篇文献的 Γ 用 meV 标定,想对比宽度时就需要换算;如果你做的是 1550 nm 通信波段器件,1 meV 大约相当于 0.8 nm 的波长宽度,这类换算是常有的事。

基线漂移是另一个高频问题。很多时候光谱仪测出来的反射谱并不在一个平坦背景上,而是叠了一个缓变的倾角。这时我建议在模型里加一个线性背景项:

I(x) = y0 + k·x + A · [(q + ε)² / (1 + ε²)]

也就是从五参数变成六参数。这个线性项能有效吸收样品表面倾斜、光源光谱分布不均、探测器响应不平坦等系统误差。但注意别把 k 加得过大,否则它会和 A、y0 之间产生很强的参数相关性,拟合出的参数误差会变大。我通常经过多次试算,只把 k 限制在一个很小的范围内。

5.3 多峰重叠与高噪声数据的处理

如果谱图里两个共振靠得很近,线型会相互叠加,单峰模型拟合出来的 Er 会偏到一个“折中”的位置,Γ 也会偏大,Q 自然失真。处理思路有两种:一种是强行用双法诺模型去拟合,也就是把两组法诺项相加,共享一个 y0;另一种是把两个峰分别用不同数据区间拟合,各自提取参数。我的经验是,如果两个峰中心距离小于 Γ 的三倍,先用双峰模型,否则优先做区间隔离,因为双峰模型的参数太多,对初值依赖非常大。

高噪声数据是拟合误差的主要来源之一。我的做法,第一步是测量时多做几次平均,这是最好的降噪手段;第二步才是在处理时做适度平滑,比如 3~5 点滑动平均,千万别用窗口大的平滑,共振峰会被抹平,Γ 会虚高;第三步是拟合时考虑使用鲁棒拟合。fitnlm 支持鲁棒拟合参数:

mdl = fitnlm(x, y_noisy, fano_fun, p0, 'Robust', 'bisquare');

鲁棒拟合能自动降低离群点的影响,对于偶尔出现的跳点特别有效。

最后再分享一个我一直在用的习惯:把整套拟合流程封装成一个函数,输入是横轴数据、纵轴数据和初始参数估计,输出是五参数、Q 值和误差。批量处理样品时,写一个循环把文件夹里几十条谱线全部算完,自动导出 Excel 表。这个习惯帮我节省了大量重复劳动,也减少了手动操作引入的错误。尝试一次,你会回来感谢这个决定的。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询