简介:SASfit是一款用于小角度散射(SAXS/SANS)数据分析和结构模型拟合的开源程序,面向材料、生物、胶体等领域从事结构表征的科研人员与工程师。软件可计算回转半径、散射不变性、Porod常数等整体结构参数,并支持多种形状因子与尺寸分布拟合;借助一组通用全局参数,还能同时拟合多条散射曲线,在对比变化实验或极化中子测量中能有效降低参数相关性,提升结果确定性。压缩包内共208个文件、约54.19MB,主要包含66个tbl表格数据、44个so动态库与44个头文件、18个gpf拟合脚本、17个rgb颜色映射表,以及csv、txt、pdf说明和def示例等,可辅助用户快速查看编译产物、依赖接口与使用示例,也便于直接利用示例脚本开展二次开发和定制。目前已有1617人学习或下载,适合希望深入理解模型拟合原理、基于开源代码扩展功能或将其集成到自身数据处理流程中的科研用户。
1. 拟合小角散射数据不是玄学:SASfit把这件事变成了可复现的工作流
做小角散射的人手里都有一堆这样的数据:样品的散射强度I(q)随着散射矢量q增大而快速衰减,想要从中得到粒径、形状和尺寸分布,就必须把合适的结构模型拟合上去。SASfit正是为此而生的开源程序,它支持SAXS和SANS两类数据,内置几十种结构模型,从实心球、核壳球、椭球、圆柱到聚合物链、微乳液结构因子都有,拟合完成直接输出参数、误差、chi²和残差曲线。它最大的价值在于把散射实验数据变成可以写进报告的结构结论——反推出半径、壳层厚度、多分散度、体积分数这些关键参数。这篇文章不是照着界面念说明书,而是把模型选型、参数约束、数据预处理和排错方法串起来,适合刚拿到第一条曲线的研究生,也适合要批量处理几十条数据的工程师。
2. 散射拟合的数学骨架:SASfit到底在算哪个方程
2.1 从I(q)到结构参数:SASfit替你解的那个方程
在稀释、无相互作用的近似下,小角散射强度对q的依赖可以写成简洁的形式:
I(q) = scale × P(q) + background
其中P(q)是形状因子,描述单个粒子的几何对散射的贡献;scale是强度缩放因子,物理上它集中了粒子数密度N、粒子体积V的平方以及粒子与溶剂之间的散射长度密度差Δρ²;background是与q无关的常数背景,用来吸收实验中未能完全扣除的溶剂散射和杂散散射。SASfit的大多数模型都采用这种归一化约定:P(q)在q趋近0时归一化到1,也就是说只描述散射曲线的形状,把强度绝对值全部交给scale。这个约定的实际含义是,即使实验数据没有做绝对强度标定,依然可以用任意单位的强度曲线拟合出尺寸参数,只是不能从scale解读出粒子数密度和体积分数。
如果体系里粒子间的相互作用不可忽略,比如高体积分数的胶体、带电粒子体系,公式变成:
I(q) = scale × P(q) × S(q) + background
S(q)是结构因子,描述粒子间的空间关联。SASfit提供了硬球、带电硬球、黏性硬球等结构因子模型,可以和任意的形状因子组合。以最简单的实心球模型为例,P(q)的解析表达式是 [3(sin(qR) − qR·cos(qR)) / (qR)³]²,其中R是球的半径,qR很小时这一项趋近1。核壳球则在这个基础上叠加壳层与核心的散射长度密度差,表达式里会出现两个不同半径的散射振幅之差,因此核壳结构会在散射曲线上留下特征性的次级振荡峰。
SASfit要做的就是用非线性最小二乘迭代,找到一组参数使计算出的模型曲线与实验数据的加权残差平方和最小。默认使用Levenberg-Marquardt算法,这个算法在参数初值离最优解不远时收敛非常快,但初值离谱或者参数高度相关时很容易陷在局部极小值里。这就是后文里参数边界和初值设置成为整个拟合流程成败关键的根源。
如果内置模型满足不了需求,SASfit还允许用户自定义P(q)的表达式。自定义模型时最关键的是保持P(q)在q→0处归一化到1,否则scale和几何参数之间会产生很强的相关性,拟合结果会变得极不稳定。我在实际使用中遇到过几次自定义模型拟合发散,最终都是因为归一化系数和scale揽在一起扯不清楚。
2.2 模型库选型逻辑:先从简单模型开始,别让参数表吓住
SASfit内置模型用几十种来描述绝不过分,从简单几何体到复杂自组装结构都有。对新手来说,真正的难点不是模型不够用,而是不知道选哪个。我的习惯是遵循一条原则:先根据实验特征缩小候选范围,然后从参数最少的模型开始。
| 实验曲线特征 | 优先考虑的模型 | 选型理由 |
|---|---|---|
| 低q区平台明显,中高q有振荡峰 | 实心球、核壳球 | 球体干涉峰明显,核壳结构会引入次级峰 |
| 低q区有-4的幂律下降 | 界面平滑致密粒子的Porod型模型 | -4幂律是平滑界面的典型特征 |
| 中q区有-2的幂律下降 | 高斯链、排除体积链 | 链状聚合物的典型散射特征 |
| 整条曲线缓慢下降,没有明显平台 | 椭球、圆柱 | 形状的各向异性使平台被抹平 |
| 微乳液、双连续结构 | 微乳液双连续结构模型 | 这类模型专门描述双连续相的峰型 |
选型之前建议先做两步快速分析,而不是直接堆模型。第一步在低q区做回转半径近似,得到整个散射体的回转半径Rg,球体的话半径大约对应sqrt(5/3)·Rg;第二步看高q区的幂律指数,指数接近-4说明粒子界面平滑,明显小于-4说明有表面粗糙或分形层次结构。这两步用任意数据处理软件都能做,做完再回到SASfit选模型,能省掉很多盲目试模型的时间。
2.3 参数分类与约束机制:不约束参数的拟合等于碰运气
SASfit的模型参数看起来五花八门,实际上可以分成四类。第一类是强度相关参数,只有scale和background两个;第二类是几何尺寸参数,球形是半径,圆柱是半径加长度,核壳结构还要加壳厚;第三类是分布参数,即多分散度,SASfit可以在几何参数上叠加Schulz、Gauss或者Lognormal分布,描述样品中尺寸不完全均一的情况;第四类是结构因子参数,包括粒子的体积分数、硬球相互作用半径等。
小角散射数据的信息量是有限的。一条q_max/q_min覆盖两个数量级的曲线,能稳定确定的独立参数一般不要超过5到7个。超过这个数,多出来的参数就会在样品数据里找解释空间,结果看起来拟合得很漂亮,换个数据集就完全不是那个结果。参数相关性是最常见的翻车源头:scale和体积分数高度纠缠,多分散度和背景参数也会互相补偿。在数据上它们各自变化却产生几乎相同的模型曲线,拟合算法自然无法分辨。
所以我在开始拟合前会做三步参数约束。第一步把物理上不可能的值封死,半径必须大于0,多分散度通常在0到0.5之间,体积分数不能超过0.74这个球密堆上限。第二步根据回转半径分析的结果设置几何参数的量级,初值取估算值的约50%作为上下限。第三步先固定所有分布参数和结构因子参数,只让几何参数自由,收敛之后再逐步放开其余参数。这个顺序能显著降低局部极小的风险,尤其适合新手。
分布函数的选择看起来是小事,实际会对参数结果产生明显影响。Gauss分布对平均值对称,适合宽度比较窄的样品;Lognormal分布允许向大尺寸方向延伸长尾,很多工业胶体和生物组装体系的长尾特征正是如此;Schulz分布解析性最好,计算快,对中等多分散度体系是个稳妥默认项。同一个数据集用三种分布拟合出的平均半径可能存在几个纳米的偏差,这并不代表拟合错误,只是分布假设不同。报告时应当写明用了哪种分布,换分布时不要直接对比参数。
提示:SASfit里每个参数都可以单独设置初始值、最小值和最大值,还可以在拟合过程中固定住。把不打算拟合的参数锁定,比删掉它更安全,因为它仍然参与模型计算,只是不参与优化。
3. 用SASfit跑通第一次拟合:数据格式、参数设置和结果检查
3.1 准备一份SASfit能直接读的数据文件
SASfit读取的散射数据通常是纯文本格式,空格或Tab分隔,三列分别是q、I(q)和误差。q的单位常见是Å⁻¹,如果实验站给的是nm⁻¹,换算时注意是10倍关系。误差这一列不能省,拟合的权重完全依赖它,误差缺失或者填0,最小二乘的权重就会变成无穷大或非有限值,结果不可信。
从原始探测器文件到SASfit输入文件,中间通常要经过辐射校正、容器扣除、样品浓度归一化等步骤。这些步骤各自的处理流程由实验站软件完成,SASfit本身不负责背景扣除,background参数只是拟合里的一个常数。我习惯把已经扣除过溶剂散射的数据送入SASfit,残留的微小背景用background参数吸收。下面这段脚本是每次拿到新数据都会先跑一遍的清洗流程:
import numpy as np # 原始数据: 假设前三列是 q, I, dI, 其他列用不到 raw = np.loadtxt('raw_saxs.dat', skiprows=2) q = raw[:, 0] I = raw[:, 1] err = raw[:, 2] # 去掉 q 过小、过大和误差非正的坏点 mask = (q > 2e-3) & (q < 0.8) & (err > 0) q, I, err = q[mask], I[mask], err[mask] # 归一到最大强度,让 scale 初值可以取 1 附近 norm = I.max() I /= norm err /= norm np.savetxt('sasfit_input.dat', np.column_stack((q, I, err)), fmt='%.6e', delimiter=' ', header='q[1/A] I dI', comments='#')脚本把强度归一化到最大值,归一化后scale的初值直接设成1附近即可,这是我的习惯做法。q的下限2e-3和上限0.8是按典型SAXS配置给的,实际要根据数据覆盖范围和信噪比调整。下限附近常常是光束阻挡器边缘污染最严重的区域,宁可多去掉一两个点,也不要让低q区坏点主导拟合。高q端通常没问题,但q太大导致计数很差时误差很大,那些点权重本来就小,留着也无妨。
3.2 配置一次完整拟合:从模型到参数边界
打开SASfit,载入刚生成的sasfit_input.dat之后,操作流程一般如下。首先在模型列表中选择带尺寸分布的球模型Sphere with size distribution,然后设置参数表的初始值:scale填1,background填0或一个很小的正数比如0.001,半径填根据低q区回转半径估算的数值,多分散度填0.1,分布函数选Schulz。接着设置参数边界:scale的上下限取0.1到10,background取-0.01到0.1,半径上下限取估算值的50%到200%,多分散度限在0.01到0.5。最后运行拟合。
这个配置里有两个细节值得说。一是background上限给到0.1,是因为预处理已经把最大强度归一化成1,0.1已经是背景占比很高的上限,给太小会把背景锁死在边界上。二是多分散度下限0.01而不是0,这是为了让参数不落在边界值上,避免某些拟合器对边界值的误差统计处理失真,也保留了一定的优化空间。
拟合结束之后,界面通常同时显示模型曲线、实验数据和残差曲线。残差曲线是重点,它显示每个数据点上实验值与拟合值的差除以误差后的结果。残差在0附近随机分布最理想,如果残差呈现明显的弯曲或周期性的起伏,意味着模型形状与实际散射形态不一致,这时候要回到模型选择,而不是继续调参数。
如果一次拟合从不同初始值出发得到明显不同的参数结果,说明模型存在参数退化性,也就是参数组合不唯一。我一般会重复三次不同初值的拟合,结果不一致就不把这次拟合当作可靠结论,转而简化模型。例如把核壳球换成实心球,或者先不加结构因子,等简单模型跑通再逐层加复杂度。
3.3 拟合质量不能只看chi²:把参数误差和物理合理性一起检验
SASfit输出里的chi²、参数估计和渐近标准误差都很有用,但很多人只盯着chi²。chi²本身是不完善的指标,它依赖误差估计的准确性。实验误差被低估,chi²虚高;误差被高估,chi²虚低。我把chi²当作参考,判断拟合品质主要看三件事:残差分布、参数误差和物理合理性。
下面的脚本可以快速完成前两种检查:
import numpy as np # 读取实验数据和SASfit导出的拟合曲线 exp = np.loadtxt('sasfit_input.dat') fit = np.loadtxt('fit_curve.dat') q = exp[:, 0] I_exp = exp[:, 1] err = exp[:, 2] I_fit = fit[:, 1] # 假设导出文件第二列是模型强度 res = (I_exp - I_fit) / err print('chi2/dof =', np.sum(res**2) / (len(res) - 5)) # 统计超过3倍误差点数 bad = np.sum(np.abs(res) > 3) print('residual outliers > 3sigma:', bad, '/', len(res)) # 检查残差是否在某一段系统性同号 for lo, hi in [(0, 5), (5, 10)]: seg = res[int(len(res)*lo/20):int(len(res)*hi/20)] print(f'residual mean q bin {lo}-{hi}:', seg.mean())这里len(res)-5是自由度,5是带尺寸分布的球模型的参数个数,实际使用时要按模型的真实参数数调整。残差段均值如果显著偏离0,说明这个q区间存在系统偏差,即使整体chi²很好看也要警惕。
物理合理性检查最简单但最容易犯。多分散度拟合出0.49而样品来自已知的窄分布单分散体系,多半是模型选错了;半径拟合出500 Å而动态光散射结果显示水合粒径15 nm,尺寸尺度完全对不上。出现这类矛盾时不要强行解释,先回头检查数据预处理和模型选择。
4. 拟合翻车血泪排查:SASfit常见问题和避坑指南
4.1 低q区拟合曲线压不住,残差在低角度整体上翘
现象:不管选什么模型,拟合曲线在低q段都低于实验数据,残差呈现明显的正峰,而且随着q减小残差越来越大。
原因:低q区上翘通常意味着体系里有大尺度结构,比如几十纳米的粒子聚集成几百纳米的团簇。单个粒子的模型里根本不存在这么大的空间尺度,所以永远压不住。另一个常见原因是样品池边缘的寄生散射没有完全扣除,这种杂散信号集中在低q区。
解决:先看数据本身,把I对q做双对数图,如果最低的几个点在斜率趋势上方突然抬头,优先怀疑是聚集体或寄生散射。实验上对样品做一次离心或过滤,重新测量空池背景并扣除。数据处理上可以从拟合区间剔除最低的几个点,范围不用很大,剔除后如果chi²明显下降,说明问题就集中在低q区。
4.2 拟合发散或参数漂移到边界,模型曲线却看不出明显变化
现象:迭代不收敛,或者半径、多分散度、scale等参数一路漂移到边界并被卡住。此时模型曲线和实验点之间可能看起来已经比较接近了。
原因:参数相关性强。多分散度和scale的组合可以产生非常相似的散射曲线:大尺度因子配宽分布与较小尺度因子配窄分布,在有限的数据区间内几乎无法区分。LM算法在这两个参数的组合谷里来回振荡,最终被顶到边界上。
解决:先固定多分散度,用一个合理的小值比如0.05,只让半径和scale自由拟合,拿到一个稳定的基础解。然后再放开多分散度,边界不要给到0.5以上,多数单峰样品不会超过0.3。如果放开之后多分散度又冲边界,说明样品本身可能是双峰分布,单峰模型根本描述不了它。
4.3 chi²接近1,但残差在中q区出现系统性正负振荡
现象:全局拟合指标很好,残差在某个q区间内连续正负交替,呈波浪状。单独看每一小段残差都不算大,但整个区间明显不是随机分布。
原因:中q区的振荡往往对应粒子间的干涉峰或核壳结构的次级峰。硬球结构因子会在q≈2π/d的位置产生一个宽峰,d是粒子间平均距离。如果模型里没有S(q),这个峰无法被描述,只能被周围参数分摊吸收,于是残差呈现波浪状。
解决:不要把这一现象当作统计噪声。在SASfit模型定义中加入硬球结构因子,拟合粒子的体积分数和作用半径。如果加入S(q)后中q区残差立即平息,说明粒子间关联确实存在。注意加入S(q)会让体积分数和scale产生相关性,需要重新收紧边界,通常先固定体积分数在实验预估浓度附近,拟合收敛后再放开。
4.4 多分散度拟合到接近0,是不是样品真的单分散
现象:多分散度参数在每次拟合中都跑到0.01的下边界,多次更换初值都是同样结果。
原因:低信噪比或者模型本身对多分散度不敏感时,分布宽度无法被数据约束,优化器总是选择更容易降低chi²的方向推参数。多分散度在数据里没有有效梯度信息,被推到边界是统计学上的必然,不能当作样品单分散的证据。
解决:把多分散度固定在几个不同的值分别拟合,比如0、0.05、0.1、0.2,比较各组的chi²。如果chi²随多分散度几乎没有变化,说明数据不支持报告这个参数,报告时应当写“多分散度未被数据约束”,而不是写“样品单分散”。我习惯在写结论时把这种参数的固定值作为假设明确列出,而不是混在拟合结果表里,避免给合作方传递错误的确定性。
5. 从能跑到能信:参数敏感性和多模型比较给SASfit结果把关
5.1 参数敏感性分析:数据到底支持哪个参数
SASfit输出的参数误差是基于协方差矩阵的渐近估计,参数相关性强时这个误差会骗人。更可靠的做法是参数扫描:选定一个参数,把它固定在一系列取值上,每次重新拟合其余参数,记录chi²随固定值的变化。如果曲线是窄而深的单谷,说明数据对这个参数有很强的约束力;如果是大平底,说明它基本是自由的。这个检查比协方差矩阵更直观,也更方便写进报告。
我用一组扫描结果做判断的脚本如下:
import numpy as np # scan.dat 两列: 固定半径值, 对应最优chi2 s = np.loadtxt('scan.dat') i = np.argmin(s[:, 1]) r_best, chi_min = s[i] # 按 chi2 <= chi_min + 1 划定近似的 1-sigma 区间 sel = s[:, 1] <= chi_min + 1 r_lo, r_hi = s[sel, 0].min(), s[sel, 0].max() print(f'最优半径 {r_best:.1f}, 近似不确定区间 {r_lo:.1f}-{r_hi:.1f}') if s[sel].shape[0] > s.shape[0] * 0.8: print('警告: 参数未被有效约束')这里的chi²最小值加1是近似一个标准误的经典做法,实际使用时如果误差列估计不准,区间意义会打折扣,但用于筛选哪些参数可信已经足够。
5.2 多模型比较:用惩罚项挡住过拟合
参数多的模型chi²总是更小,直接比较不公平。我习惯用BIC(贝叶斯信息准则)把参数数量惩罚进去,数值越低越好。SASfit本身不计算BIC,但拿到每次拟合的chi²、自由参数数和数据点数后,几行代码就能算清楚:
import numpy as np def bic(chi2, k, n): return chi2 + k * np.log(n) # 模型A: 实心球 4个自由参数, chi2=320, 数据点200 print('实心球 BIC:', bic(320, 4, 200)) # 模型B: 核壳球 6个自由参数, chi2=298 print('核壳球 BIC:', bic(298, 6, 200))核壳球chi²小了22,但BIC只小了约12,优势并不像chi²差距看起来那么大。实际使用时如果模型间BIC差异小于2,我会视为没有决定性差异,选择参数更少的那个作为最终结论。多模型对比要在同一批数据、同一套误差权重下进行,换了一条误差列,BIC的绝对数值没有可比性。
我在一次批次样品分析中,为了让拟合更完美,叠加了核壳层、硬球结构因子和尺寸分布,结果BIC排序表明最简洁的实心球模型反而是多数样品的最佳选择。从那以后我形成了一个习惯:每个样品至少跑一个简单模型和一个复杂模型,把BIC结果列入汇总表再决定哪些参数进入最后的讨论,而不是挑一条最漂亮的拟合曲线交差。这个习惯帮我挡掉了不少拟合得很漂亮但一推就倒的结论。希望帮到你。
本文还有配套的精品资源,点击获取