稀疏突发计数数据预测：SARIMAX与负二项回归在漏洞活动预测中的实战对比
1. 项目概述：当漏洞数据变得“吝啬”时
在安全运营中心（SOC）或者漏洞管理团队里，我们常常面临一个尴尬的局面：警报要么不来，一来就是一大堆。但更多时候，我们面对的是另一种更令人头疼的数据——稀疏且突发。想象一下，你负责监控一个大型企业网络，每天从各种扫描器、蜜罐、入侵检测系统（IDS）和威胁情报源涌来的原始日志可能有TB级，但经过过滤、去重和验证，真正能确认的、针对特定资产或漏洞的“有效攻击尝试”事件，可能一天就只有寥寥几条，甚至连续几天风平浪静。然后，在某个毫无征兆的下午，突然因为一个刚曝出的高危漏洞（比如某个流行中间件的远程代码执行漏洞），相关攻击事件在几小时内激增到上百条，随后又迅速归于平静。
这种数据模式，就是典型的“稀疏与突发”。它不像服务器CPU利用率或者网站日活用户数那样平滑、连续，可以用经典的时间序列模型（如ARIMA）轻易拿捏。传统安全事件预测模型，如果直接套用常规方法，往往会“水土不服”：在平静期（零值或接近零值）预测出毫无意义的低风险，而在爆发期又严重滞后，无法在攻击浪潮真正来临前给出有效预警。
我最近投入了大量时间，专门研究如何基于这种“吝啬”的数据来预测漏洞活动。核心思路是，不把漏洞攻击事件数当作一个普通的连续数值来预测，而是正视其“计数”的本质——它永远是非负整数（0， 1， 2， …），并且存在大量零值。这自然引出了两个技术路线的对比：一方是经过改造以适应计数数据的时间序列模型，另一方则是天生为计数而生的计数模型。本文将深入拆解这两种思路，结合具体的模型（如SARIMAX和泊松回归）进行实战对比，分享从数据预处理、模型选型、调参到效果评估的全过程，以及我踩过的那些坑。
无论你是安全分析师、威胁情报研究员，还是对时序预测感兴趣的数据科学家，这篇文章都将为你提供一个处理稀疏突发计数数据的完整框架和实操指南。
2. 核心思路与模型选型背后的逻辑
面对漏洞活动预测，我们首先要明确预测目标。通常，我们可能想预测：“未来24小时内，针对某特定CVE编号漏洞的攻击尝试次数”或“未来一周内，企业内部被扫描探测的端口次数”。这个目标变量是一个时间序列，但每个时间点（如每小时、每天）的值是一个计数。
2.1 为什么通用时间序列模型会“失灵”？
经典时间序列模型，如自回归积分滑动平均模型（ARIMA）及其季节性版本（SARIMA），假设数据背后是一个连续的、近似正态分布的随机过程。它们通过历史值（AR）、历史误差（MA）和差分（I）来捕捉趋势和周期性。但当数据是稀疏计数时，问题接踵而至：
- 负值预测：ARIMA模型可能会预测出负的攻击次数，这显然没有物理意义。
- 对零值不敏感：模型将零视为一个普通的低数值，无法有效建模“长时间无事件”这种状态背后的机制（可能是防御有效，也可能是攻击者暂时未关注）。
- 突发处理能力差：当出现一个巨大的离群值（爆发）时，模型需要很长时间才能“学习”到这个突变，导致预测严重滞后，且可能过度拉高后续平静期的预测值。
为了解决这些问题，我们引入了SARIMAX模型。这里的“X”代表外生变量。核心改进思路是：我们仍然使用时序模型的结构来捕捉自相关性和季节性，但通过连接函数（Link Function）和特定的误差分布，使其适应计数数据。一种常见做法是假设数据服从泊松分布或负二项分布，并使用对数连接函数，确保预测值为正。这本质上是将通用时间序列模型“包装”成一个广义线性模型（GLM）的框架。它的优势在于能直接利用成熟的时间序列建模能力处理季节性和趋势。
2.2 为什么计数模型是“天选之子”？
另一条路是直接使用专为计数数据设计的回归模型，其代表是泊松回归和负二项回归。它们从根上就承认数据是非负整数的特性。
- 泊松回归：假设在给定的时间区间内，事件发生的次数服从泊松分布。其核心特征是期望（均值）等于方差。它通过一个对数连接函数，将预测的均值与一系列特征（包括时间特征）线性关联起来。
- 负二项回归：可以看作是泊松回归的扩展。它放松了“均值=方差”这个强假设，引入了一个离散参数来处理“过离散”现象——即数据的方差远大于均值，这在突发性攻击数据中极其常见。
对于时间序列预测，我们需要将这些计数模型“时间序列化”。方法很简单：将时间特征作为核心的自变量。例如，对于日度数据，我们可以构造以下特征：
day_of_week(周一至周日，独热编码)is_weekend(是否周末)month(月份，考虑季节性)time_index(一个简单的线性趋势项)lag_1,lag_2,lag_7(过去1天、2天、7天的攻击次数，用来捕捉自相关性)
这样，计数模型就具备了处理时间依赖的能力。它的优势在于模型假设更贴近数据本质，对零值和过离散有天然的良好处理能力，模型解释性也更强。
2.3 模型选型决策树
在实际项目中，我的选型逻辑通常遵循以下路径：
graph TD A[开始：稀疏突发计数数据] --> B{数据探索：检查过离散程度}; B -- 方差 ≈ 均值 --> C[选择泊松回归]; B -- 方差 >> 均值（常见） --> D[选择负二项回归]; C & D --> E{是否存在强季节性或复杂趋势？}; E -- 是， 且季节模式稳定 --> F[优先尝试SARIMAX（计数版本）]; E -- 否， 或季节模式复杂/多变 --> G[优先尝试计数模型+时间特征]; F --> H[模型训练与评估]; G --> H; H --> I{验证集效果对比}; I -- SARIMAX更优 --> J[采用SARIMAX， 关注外生变量]; I -- 计数模型更优 --> K[采用计数模型， 可尝试加入滞后项];注意：这个决策树不是绝对的。最可靠的方法是进行基准测试。我会同时构建SARIMAX（如基于泊松分布的）和负二项回归模型，在一个预留的验证集上比较它们的预测性能。性能指标不能只看均方误差（MSE），因为对计数数据不敏感，应重点关注平均绝对误差（MAE）、泊松偏差或零膨胀部分的预测准确率。
3. 数据准备与特征工程的魔鬼细节
模型再高级，垃圾数据进去，出来的也是垃圾。对于漏洞活动数据，预处理和特征工程至关重要，且与传统连续值时序数据有所不同。
3.1 数据收集与聚合
数据通常来自多个源头：
- 漏洞扫描日志：Nessus, OpenVAS, Qualys等工具的每日扫描结果，按CVE或资产聚合。
- 网络流量日志：IDS/IPS（如Suricata, Snort）警报，防火墙拒绝日志。
- 蜜罐交互数据：记录攻击者探测、利用尝试的次数。
- 威胁情报订阅：关于特定漏洞被活跃利用的IOC（失陷指标）数量。
第一步是时间对齐与聚合。确定预测的最小时间单位（如1小时、1天）。将原始事件日志按此单位进行聚合，计算每个时间桶内的“攻击事件计数”。这里就会产生大量的零值桶。
3.2 关键特征构建
除了上一节提到的基础时间特征，对于安全领域，以下特征往往有奇效：
漏洞生命周期特征：
days_since_disclosure: 漏洞公开至今的天数。通常刚公开时活动激增，随后衰减。has_exploit_public: 布尔值，表示是否有公开的利用代码（PoC/Exp）。这是最强的活动驱动因子之一。severity_score: CVSS基础评分。高分漏洞更受关注。asset_value: 受影响资产的重要性权重（需内部定义）。
外部威胁情报特征：
ti_mentions_24h: 过去24小时内在特定威胁情报论坛、Twitter、暗网提及该漏洞的次数（可通过API获取）。malware_family_activity: 已知利用该漏洞的恶意软件家族的近期活动水平。
滞后特征与滚动统计：
lag_1,lag_2,lag_3,lag_7: 直接滞后项，捕捉短期自相关。rolling_mean_7: 过去7天的移动平均，反映近期基线。rolling_max_7: 过去7天的最大值，捕捉近期是否已有过爆发。zero_run_length: 当前零值持续的长度。这个特征对于区分“常态平静”和“爆发后的平静”很有用。
3.3 处理“零膨胀”问题
我们的数据可能包含两种零：一种是“真零”（攻击者确实没来），另一种是“漏报零”（由于检测能力限制没发现）。当零值比例极高时（如>80%），就需要考虑零膨胀模型，如零膨胀泊松回归（ZIP）或零膨胀负二项回归（ZINB）。这些模型内部包含两个过程：一个逻辑回归过程判断当前时间点是否“必然为零”（即处于非活动状态），另一个计数过程（泊松/负二项）在非“必然为零”时生成计数。
实操心得：不要一上来就用零膨胀模型。先拟合一个普通的负二项回归，检查残差。如果模型在零值区域系统性预测不佳，再考虑升级到ZINB。零膨胀模型参数更多，更易过拟合，尤其是在数据量不大的情况下。
4. 模型实战：SARIMAX vs. 负二项回归
我们以一个模拟数据集为例，假设我们预测某企业面向互联网的SSH服务每日遭受的暴力破解尝试次数。数据具有每周周期性（周末较低），且偶尔因新的暴力破解工具发布而突发激增。
4.1 基于Statsmodels的计数SARIMAX实现
Python的statsmodels库的SARIMAX模块支持通过exog参数引入外生变量，并通过指定mle_regression=False和使用log连接函数来近似实现计数特性。但更严谨的做法是使用GeneralizedLinearModel（GLM）框架下的Poisson或NegativeBinomial家族，并手动加入ARIMA结构（这比较复杂）。这里展示一个使用statsmodels的GLM与滞后项结合的方法，模拟计数时序模型。
import pandas as pd import numpy as np import statsmodels.api as sm from statsmodels.genmod.generalized_linear_model import GLM from statsmodels.genmod.families import NegativeBinomial # 假设 df 是一个包含‘date’和‘attack_count’的DataFrame df['time_index'] = range(len(df)) df['day_of_week'] = df['date'].dt.dayofweek # 创建滞后特征 for lag in [1, 2, 7]: df[f'lag_{lag}'] = df['attack_count'].shift(lag) # 创建周末标识 df['is_weekend'] = df['day_of_week'].isin([5, 6]).astype(int) # 删除因滞后产生的NaN行 df_model = df.dropna().copy() # 定义特征（X）和目标（y） # 这里将滞后项、时间趋势、周几都作为外生变量 exog_features = ['time_index', 'is_weekend'] + [f'lag_{i}' for i in [1, 2, 7]] # 为周几创建虚拟变量（独热编码） exog_with_dummies = pd.get_dummies(df_model[['day_of_week'] + exog_features], columns=['day_of_week'], drop_first=True) X = sm.add_constant(exog_with_dummies) # 添加常数项 y = df_model['attack_count'] # 使用负二项回归（GLM框架） # 注意：statsmodels的NegativeBinomial家族默认使用log连接，且需要指定alpha参数（离散参数） # 这里使用‘nb2’方差形式，并让模型通过MLE估计alpha try: nb_model = sm.GLM(y, X, family=sm.families.NegativeBinomial(alpha=1.0)) # 初始化alpha nb_result = nb_model.fit(method='lbfgs', maxiter=1000) # 使用优化器 print(nb_result.summary()) except Exception as e: print(f"GLM拟合失败: {e}") # 回退到泊松回归，它更稳定但可能不适合过离散数据 poisson_model = sm.GLM(y, X, family=sm.families.Poisson()) poisson_result = poisson_model.fit() print(poisson_result.summary())这段代码的关键在于将时间序列结构（滞后项、趋势、季节性）通过特征工程的方式，融入到广义线性模型中。这本质上就是一个“计数模型”的思路，但它清晰地利用了时间序列的依赖关系。
4.2 基于Scikit-learn与标准化的负二项回归
对于更复杂的特征工程和管道化操作，scikit-learn的兼容性更好。我们可以使用statsmodels的离散模型，或scikit-learn的PoissonRegressor和NegativeBinomialRegressor（需sklearn版本支持或使用glmnet等库）。这里以类scikit-learn的方式组织：
from sklearn.linear_model import PoissonRegressor from sklearn.preprocessing import StandardScaler, OneHotEncoder from sklearn.compose import ColumnTransformer from sklearn.pipeline import Pipeline from sklearn.model_selection import TimeSeriesSplit # 准备特征和数据 features = ['time_index', 'is_weekend', 'lag_1', 'lag_2', 'lag_7', 'day_of_week'] X_raw = df_model[features] y = df_model['attack_count'] # 定义预处理：数值特征标准化，分类特征独热编码 numeric_features = ['time_index', 'lag_1', 'lag_2', 'lag_7'] categorical_features = ['day_of_week', 'is_weekend'] preprocessor = ColumnTransformer( transformers=[ ('num', StandardScaler(), numeric_features), ('cat', OneHotEncoder(drop='first'), categorical_features) ]) # 创建泊松回归管道 poisson_pipeline = Pipeline(steps=[ ('preprocessor', preprocessor), ('regressor', PoissonRegressor(alpha=0.01, max_iter=1000)) # alpha是L2正则化强度 ]) # 使用时间序列交叉验证 tscv = TimeSeriesSplit(n_splits=5) mae_scores = [] for train_idx, test_idx in tscv.split(X_raw): X_train, X_test = X_raw.iloc[train_idx], X_raw.iloc[test_idx] y_train, y_test = y.iloc[train_idx], y.iloc[test_idx] poisson_pipeline.fit(X_train, y_train) y_pred = poisson_pipeline.predict(X_test) mae = np.mean(np.abs(y_test - y_pred)) mae_scores.append(mae) print(f"泊松回归平均MAE: {np.mean(mae_scores):.2f}")注意事项：对于计数模型，特别是泊松回归，特征缩放很重要。虽然模型本身对线性变换的尺度不敏感（因为连接函数是非线性的），但正则化项（如L2）会受到特征尺度影响。使用
StandardScaler是个好习惯。另外，PoissonRegressor预测的是计数的期望值（均值），通常是浮点数，需要根据业务决定是四舍五入还是保留小数。
4.3 模型评估：超越MAE的指标
对于稀疏突发数据，MAE可能被大量的零值所主导，掩盖了模型对突发峰值预测的好坏。我们需要多维度评估：
- 按计数区间评估：将测试集的真实值分成几个区间，如
[0],[1-5],[6-20],[20+]，分别计算每个区间内的MAE或预测偏差。这能看出模型在平静期和活跃期的表现。 - 峰值捕获能力：定义一个“峰值”阈值（如历史计数的95分位数）。计算模型预测值超过该阈值的时刻中，有多少比例与真实峰值时刻重合（精确率），以及真实峰值有多少被预测出来（召回率）。
- 概率校准评估：对于泊松/负二项模型，它们本质上是概率模型。我们可以检查预测分布的可靠性。例如，对于所有预测均值为λ的点，真实值落在预测分布（如泊松(λ)）的某个置信区间（如90%）内的比例是否接近90%。
- 零值预测准确率：单独计算模型预测为零且真实值为零的准确率，这对运营成本（减少误报警）很有意义。
5. 常见陷阱、调试心得与进阶思考
在实际操作中，我遇到了无数坑，这里分享几个最典型的。
5.1 过离散与模型选择失误
问题：一开始使用了泊松回归，模型训练似乎正常，但预测值方差远小于真实值方差，对突发峰值的预测严重不足。诊断：绘制真实值的均值-方差关系。计算np.var(y_true) / np.mean(y_true)，如果远大于1（比如>2），则存在明显的过离散。解决：立即切换到负二项回归。负二项回归引入了一个离散参数（在statsmodels中常称为alpha），专门捕捉超出泊松假设的变异。如果使用scikit-learn的PoissonRegressor，可以考虑换成支持负二项分布的库，或者尝试使用准泊松回归（它允许方差与均值成比例，但不指定完整分布）。
5.2 外生变量的“未来泄露”
问题：在构建“过去24小时威胁情报提及数”这个特征时，不小心使用了包含预测时间点在内的数据，导致模型在训练和验证时看到了“未来”信息，得到虚假的高精度。解决：对于任何需要聚合计算的特征，必须严格进行时间窗滞后。例如，预测t时刻的攻击数，使用的特征数据只能来自t-1时刻及以前。在代码中，这通常意味着对特征列使用.shift(1)操作。在管道中，可以使用sklearn的FunctionTransformer自定义特征计算，确保时序安全。
5.3 突发导致的历史数据“污染”
问题：一次巨大的攻击爆发（一个离群值）会显著拉高移动平均（rolling_mean）等特征，导致爆发过后很长一段时间，模型仍会预测较高的活动水平，形成“长尾假警报”。解决：
- 使用稳健的滚动统计：用中位数（
rolling_median）代替均值，或用修剪均值。 - 对特征进行缩尾处理：将滚动统计值中超过99分位数的值截断到99分位数。
- 引入“爆发后衰减”特征：创建一个标志位，在爆发后的N个时间段内设为1，并作为一个特征输入模型，让模型学习爆发后的衰减模式。
5.4 模型持续学习与更新
漏洞威胁 landscape 变化很快。一个静态模型很快就会过时。必须建立模型重训机制。
- 滑动窗口训练：始终只用最近N个月的数据训练模型。
- 性能监控与触发重训：持续监控模型在最新数据上的预测误差（如MAE）。当误差连续多日超过阈值，或遇到概念漂移（如新型攻击模式出现）时，自动触发模型重训。
- A/B测试框架：将新训练的模型与当前生产模型进行A/B测试，在影子模式下运行一段时间，确认其性能提升后再切换。
5.5 从预测到决策：设定动态阈值
模型输出的是未来一段时间内攻击次数的期望值。安全运营需要的是“是否报警”的二元决策。直接将预测值四舍五入后与固定阈值（如>5次）比较，效果往往不好。更好的做法：根据预测值，计算攻击次数超过某个业务可接受阈值（如k）的概率。例如，对于泊松分布，P(X > k) = 1 - sum_{i=0}^{k} (exp(-λ) * λ^i / i!)。然后，根据运营成本（误报成本 vs. 漏报成本），设定一个概率阈值（如P(X > 5) > 0.3则报警）。这个概率阈值可以通过历史数据上的成本效益分析来优化。
最终，选择SARIMAX还是计数模型，没有银弹。在我的多数实践中，特征工程良好的负二项回归（或零膨胀负二项回归）因其灵活性、解释性和对过离散数据的天然适配性，往往成为首选。而SARIMAX在具有非常强且稳定的季节性模式（如每日、每周规律极其明显）时，可能更有优势。最重要的永远是：理解你的数据，构建合理的评估框架，然后用实验结果说话。这个从稀疏突发的安全日志中挖掘预测价值的过程，本身就是一场对抗数据不确定性的精彩攻防。