企业官网建设流程全解析

做传染病模型参数优化的人，大概率都摸过SEIR模型和哈里斯鹰算法这个组合。SEIR的四个舱室关系很清晰，但传播率β、潜伏期倒数σ、恢复率γ真实取值怎么定，永远是绕不开的坑。我最早靠手工调参，后来用Matlab自带的局部优化器，效果都不稳。真正让参数拟合这件事跑顺的，是用哈里斯鹰算法（HHO）把参数搜索自动化。这篇文章完整复盘了HHO-SEIR的参数优化思路、Matlab代码和实操经验，适合需要做数据拟合的学生、研究人员，以及想快速上手的开发者。

1. 问题拆解：SEIR模型的参数为什么难定

1.1 SEIR模型原理与四个关键参数

SEIR模型是在经典SIR模型中间加了暴露舱室E，用来描述“已经接触病毒但还没具备传染性”的人群。四个舱室分别是易感者S、暴露者E、感染者I和康复者R，总人群数N = S + E + I + R。模型用四个常微分方程描述状态变化：

• dS/dt = -βSI/N，易感者因为接触感染者而减少
• dE/dt = βSI/N - σE，进入暴露状态，再以速率σ进入感染状态
• dI/dt = σE - γI，感染者以速率γ康复
• dR/dt = γI，康复者不再参与传播

这里面最关键的参数是β、σ、γ：β是有效接触率，决定传播速度；σ是1/平均潜伏期，决定暴露人群什么时候转成感染者；γ是1/平均传染期，决定感染者多久康复。真实观测里通常只有每日感染人数或者累计确诊数，E和初始潜伏人数E0都不可直接观测。问题就变成了：只知道I曲线的一部分观测值，要反推出一组能让模型尽可能吻合数据的β、σ、γ和E0。这本质上是一个带约束的非线性优化问题，而且目标函数常常是多峰的。

1.2 传统参数估计方法的局限

我最初处理这个问题时，用的是网格枚举加局部最小二乘。网格枚举在参数少的情况下还能跑，但三个参数各取几十个候选值，组合数量轻松上百万，而且等间隔网格很容易漏掉狭窄的最优区域。后来换Matlab的fmincon，局部收敛能力确实强，但目标函数有很多局部极小值，初始点稍微偏一点就会陷进去，跑出来的参数解释不了数据。用梯度类方法还有个麻烦：ODE数值积分本身带有误差，目标函数对参数的梯度近似经常抖动，导致迭代方向反复横跳。对SEIR这种小规模但强非线性的反演问题，全局优化算法往往比局部算法更实用。

1.3 为什么选哈里斯鹰算法

哈里斯鹰优化算法是2019年Heidari等人提出的元启发式算法，灵感来自哈里斯鹰群体合作围猎兔子的行为。它有几个特性很契合SEIR参数优化：不需要目标函数梯度，适用面宽；有动态的探索开发切换机制，能平衡全局搜索和局部精细搜捕；控制参数少，Matlab实现成本比遗传算法和粒子群低不少。我用同一组模拟数据对比了HHO、PSO（粒子群）、GA（遗传算法），三者固定相同的种群规模和迭代次数，HHO找到的均方误差更小，尤其在边界附近的稳定性更直观。简单把几种常见方案对比一下：

这里不是说HHO万能，而是面对“参数少、目标函数多峰、有边界约束”的SEIR反演场景，它是一个非常实用的选择。

2. HHO算法核心机制：从围猎策略到参数搜索

2.1 逃逸能量与探索开发切换

HHO最讨巧的设计是用兔子逃逸能量E来控制算法阶段的切换。E的计算是E = 2E0(1 - t/T)，其中E0是每次迭代随机生成的[-1,1]之间的数，t是当前迭代数，T是最大迭代数。当|E|≥1时，种群处于探索阶段，在参数空间里大范围搜索；当|E|<1时，进入开发阶段，围绕当前最优解精细搜捕。注意E0每次迭代都重新随机，所以同一算法在不同运行里有不同的探索开发节奏，这很像真实围猎中兔子随时改变逃跑方向，敌人也得随机应变。这种机制对多峰目标函数尤其友好，直接把“全局搜索够不够”和“局部挖掘细不细”揉在同一个公式里。

2.2 探索阶段的两种位置更新公式

探索阶段公式有两个分支。第一个分支是随机挑一只鹰的当前位置X_rand，然后生成新位置：X_new = X_rand - r1|X_rand - 2r2X_i|，这里的r1、r2是[0,1]均匀随机数，本质是在随机个体附近做扰动。第二个分支是结合种群均值X_mean和当前最优解X_rabbit：X_new = (X_rabbit - X_mean) - r3(lb + r4(ub-lb))，其中r3、r4是随机数，lb、ub是参数下界和上界向量。这个公式利用全局统计信息把新个体拉到更合理的区域，同时保留了随机跳跃能力，避免整个种群过早聚集到某个局部区域。实际跑SEIR优化时，第二分支对边界附近的搜索贡献很大。

2.3 开发阶段的包围策略与Levy飞行

进入开发阶段后，算法根据随机数r和|E|的大小分成四条路线：软包围、硬包围、软包围加渐进式俯冲、硬包围加渐进式俯冲。前两种相当于在兔子周围收紧包围圈，后两种会引入Levy飞行，即用一个重尾分布生成大尺度随机步长，让鹰在局部搜捕时偶尔跳开很远的距离。Levy飞行是HHO跳出局部坑的关键工具，尤其在SEIR参数拟合里，适应度景观往往存在很多由ODE数值解误差造成的“假谷”，Levy的远距离跳跃能有效摆脱这些假谷。把HHO套用到SEIR参数优化时，每个鹰个体就是一组[β, σ, γ, E0]，适应度函数是模型输出的I序列与观测I序列的均方误差（MSE），算法最终返回的是让MSE最小的X_rabbit。

这里额外说一点：SEIR参数优化不是标准的连续光滑问题。因为目标函数要调用ode45，数值积分误差会随参数变化，适应度景观会出现不规则的起伏。HHO这种无梯度算法不会去计算导数，所以天然能容忍这种不光滑性，这也是我最终选它的重要原因。

3. Matlab实现：HHO-SEIR优化代码拆解

3.1 生成模拟观测数据

下面的代码先用一组已知真参数生成SEIR感染曲线，并加5%高斯噪声作为观测数据。强烈建议先用模拟数据验证算法，再碰真实数据，否则很难判断优化结果可信度。我习惯把这段生成流程放在main脚本开头，每次运行都能重现实验。

% 生成SEIR模拟数据（作为观测值） clear; clc; rng(2024); N = 100000; % 总人口 I0 = 10; % 初始感染者 E0_true = 30; % 真实初始潜伏者 S0 = N - I0 - E0_true; % 初始易感者 R0 = 0; y0_true = [S0; E0_true; I0; R0]; beta_true = 0.45; % 有效接触率 sigma_true = 1/5.2; % 潜伏期约5.2天 gamma_true = 1/10; % 恢复期约10天 tdata = (0:1:60)'; % 观测60天，按天采样 [~, Y_true] = ode45(@(t,y) seir_model(t,y,beta_true,sigma_true,gamma_true), tdata, y0_true); Idata = Y_true(:,3) + randn(length(tdata),1)*max(Y_true(:,3))*0.05; % 加噪声 Idata(Idata < 0) = 0; % 实际数据不可能为负

注意tdata用了列向量，ode45会返回对应时间点的解。噪声强度用峰值感染人数的5%，太小会让优化问题过于简单，无法检验鲁棒性。

3.2 SEIR微分方程函数

Matlab里解SEIR必须有一个ODE函数，这个函数和理论微分方程一一对应。

function dydt = seir_model(t, y, beta, sigma, gamma) % 输入：t时间，y=[S;E;I;R]，输出导数 N = sum(y); S = y(1); E = y(2); I = y(3); R = y(4); dS = -beta * S * I / N; dE = beta * S * I / N - sigma * E; dI = sigma * E - gamma * I; dR = gamma * I; dydt = [dS; dE; dI; dR]; end

这里我用sum(y)计算总人口，而不是外部传入固定N，因为ODE数值解会有极小漂移，动态取和能保证比例关系始终一致。时间单位必须和参数量纲匹配，观测间隔是1天，参数就是“每天的量”。

3.3 目标函数：用ode45算MSE

目标函数是连接HHO和SEIR模型的桥梁。它的作用是传入一组待优化参数，解出整个感染曲线，然后计算模型I序列与观测I序列的均方误差。

function err = seir_fit_objective(params, tdata, Idata, N, I0) % params = [beta, sigma, gamma, E0] beta = params(1); sigma = params(2); gamma = params(3); E0 = params(4); % 边界安全检查 if beta <= 0 || sigma <= 0 || gamma <= 0 || E0 <= 0 || E0 >= N err = 1e10; return; end S0 = N - I0 - E0; if S0 <= 0 err = 1e10; return; end y0 = [S0; E0; I0; 0]; [~, Y] = ode45(@(t,y) seir_model(t,y,beta,sigma,gamma), tdata, y0); Imodel = Y(:,3); err = mean((Imodel - Idata).^2); % MSE end

参数组合里包含E0，是因为初始暴露人群对感染曲线前期的上升斜率有直接影响，不优化E0，模型很难同时匹配“起峰时间”和“峰值强度”。MSE作为目标，对峰值附近的误差更敏感，这正好符合我们对疫情拟合的直观期待。如果某组参数导致ODE求解发散，ode45可能返回NaN，此时把err设成1e10兜底，算法会自然淘汰这些坏解。

3.4 HHO主算法与主脚本

下面这段是HHO的核心实现。标准HHO包含探索、软包围、硬包围、渐进式俯冲四类策略，我这里用嵌套函数实现了Levy飞行，并把越界个体拉回边界内随机位置，避免种群堆积在边界上。

function [best_params, best_err, conv] = hho_seir(lb, ub, dim, Npop, MaxIt, tdata, Idata, N, I0) % HHO优化SEIR参数 % lb/ub：参数下界/上界向量；dim：参数维度；Npop：种群数；MaxIt：迭代数 % 初始化种群 X = repmat(lb', Npop, 1) + rand(Npop, dim) .* repmat((ub - lb)', Npop, 1); fit = zeros(Npop, 1); for i = 1:Npop fit(i) = seir_fit_objective(X(i,:), tdata, Idata, N, I0); end [best_err, idx] = min(fit); rabbit = X(idx, :); conv = zeros(MaxIt, 1); for t = 1:MaxIt % 猎物逃逸能量 E0 = 2 * rand - 1; E = 2 * E0 * (1 - t / MaxIt); J = 2 * (1 - rand); % 跳跃强度 for i = 1:Npop r = rand; if abs(E) >= 1 % 探索阶段 q = rand; if q >= 0.5 k = randi(Npop); X_new = X(k,:) - rand(1,dim) .* abs(X(k,:) - 2*rand(1,dim).*X(i,:)); else X_mean = mean(X, 1); X_new = (rabbit - X_mean) - rand(1,dim) .* (lb + rand(1,dim).*(ub - lb)); end else % 开发阶段 if r >= 0.5 && abs(E) >= 0.5 % 软包围 X_new = rabbit - X(i,:) - E .* abs(J .* rabbit - X(i,:)); elseif r >= 0.5 && abs(E) < 0.5 % 硬包围 X_new = rabbit - E .* abs(rabbit - X(i,:)); elseif r < 0.5 && abs(E) >= 0.5 % 软包围加渐进式俯冲（先试探一次） X_mean = mean(X, 1); X1 = rabbit - E .* abs(J .* rabbit - X_mean); if seir_fit_objective(X1, tdata, Idata, N, I0) < fit(i) X_new = X1; else X2 = rabbit - E .* abs(J .* rabbit - X_mean) + 0.01 * levy(dim); if seir_fit_objective(X2, tdata, Idata, N, I0) < fit(i) X_new = X2; else X_new = X(i,:); end end else % 硬包围加渐进式俯冲 X_mean = mean(X, 1); X1 = rabbit - E .* abs(J .* rabbit - X_mean); if seir_fit_objective(X1, tdata, Idata, N, I0) < fit(i) X_new = X1; else X2 = rabbit - E .* abs(J .* rabbit - X_mean) + 0.01 * levy(dim); if seir_fit_objective(X2, tdata, Idata, N, I0) < fit(i) X_new = X2; else X_new = X(i,:); end end end end % 边界处理：越界后拉回边界内的随机位置 for d = 1:dim if X_new(d) < lb(d) || X_new(d) > ub(d) X_new(d) = lb(d) + rand * (ub(d) - lb(d)); end end % 更新适应度并替换 fnew = seir_fit_objective(X_new, tdata, Idata, N, I0); if fnew < fit(i) X(i,:) = X_new; fit(i) = fnew; end end [best_err, idx] = min(fit); rabbit = X(idx, :); conv(t) = best_err; end best_params = rabbit; function L = levy(dim) beta = 1.5; sigma_num = (gamma(1+beta) * sin(pi*beta/2) / (gamma((1+beta)/2) * beta * 2^((beta-1)/2)))^(1/beta); u = randn(1, dim) * sigma_num; v = randn(1, dim); step = u ./ abs(v).^(1/beta); L = step; end end

两点说明：第一，标准HHO的渐进式俯冲分两步试探，很多简化版本直接省了，但实测模拟数据时，加上这两步能显著提高开发阶段搜捕效率。第二，Levy飞行里维度dim要和参数维度一致，不要写死成标量，否则多维更新时会报尺寸错误。

主脚本调用代码：

% 边界：beta、sigma、gamma、E0的逻辑范围 lb = [0.1, 1/15, 1/20, 1]; ub = [0.9, 1/3, 1/5, 5000]; dim = 4; Npop = 30; MaxIt = 120; [best_params, best_err, conv] = hho_seir(lb, ub, dim, Npop, MaxIt, tdata, Idata, N, I0); fprintf('优化结果：beta=%.4f, sigma=%.4f, gamma=%.4f, E0=%.2f\n', best_params(1), best_params(2), best_params(3), best_params(4)); % 用优化参数重跑模型并画图 [~, Y_best] = ode45(@(t,y) seir_model(t,y,best_params(1),best_params(2),best_params(3)), tdata, [N-best_params(4)-I0; best_params(4); I0; 0]); figure; plot(tdata, Idata, 'o'); hold on; plot(tdata, Y_best(:,3), 'LineWidth', 1.5); grid on; xlabel('天数'); ylabel('感染人数'); legend('观测数据', 'HHO优化拟合', 'Location', 'best'); title('HHO-SEIR参数优化结果');

这里观测数据画成散点，模型曲线画成实线，一眼就能看出拟合优劣。

4. 实验结果与优化参数解读

4.1 一次典型优化的参数还原情况

跑上述代码时，一组典型结果如下：

β和σ的反演精度比较高，γ相对误差略大。原因是恢复速率主要影响感染曲线末尾的衰减斜率，而末端数据通常更平缓，对γ的灵敏度低于峰期参数。如果目标是评估基本再生数R0=β/γ，真值是4.5，优化值是4.33，误差约3.8%，仍然可接受。这种参数间的联合后验相关性是SEIR反演里最常见的现象，优化算法只能给一个点估计，更严格的做法还需要不确定性量化，但那不是本文重点。

4.2 与PSO和GA的收敛速度对比

我用同样的数据、种群规模和迭代次数对比了HHO、PSO、GA。收敛曲线的大致规律是：前30代HHO的误差下降速度中等，但到60代以后还能继续突破几个局部平台，最终误差往往低于其他两者。PSO前期下降快，后期容易早熟；GA依赖交叉变异，参数多了收敛偏慢。需要强调的是，这个结论基于我自己的实现和参数设定，不代表所有场景都如此。做算法对比时，一定要固定相同的函数评估次数，否则结论没有可比性。从实际效果看，HHO在SEIR这种四维参数空间上的稳定性是最好的。

4.3 边界设定对优化效果的影响

SEIR参数边界一定要基于流行病学常识来设，不是越宽越好。σ对应潜伏期，数据按天采样时潜伏期通常3到14天，那么σ大概在1/14到1/3之间，把上界放宽到1/2虽然也能收敛，但会浪费大量函数评估次数。β我设为0.1到0.9，如果数据爆发速度特别快，再适当放宽。E0边界用观测首日感染人数的0到N，并且不要取到0，否则模型初始状态无法触发传播。每次跑完优化，除了看误差和曲线，还要看E0和初始感染人数之间的关系。如果E0被推得很大、S0被压得很低，说明模型正尝试用大量潜伏者解释前期缓慢增长，这往往意味着σ偏小，需要调整σ的边界重新优化。

5. 实操中的常见问题与排查技巧

5.1 ode45求解失败或耗时过长

最常遇到的问题是目标函数传入某些参数组合后，ode45返回NaN或者直接报错。我习惯把问题分成两类。一类是参数直接导致模型发散，比如β过大、σ过小，感染人数短时间内爆炸，数值解法步长自动缩到极小；这类要在目标函数里加有效性判断，要么返回大误差，要么用初始状态合法性过滤。另一类是时间跨度太长或数据点太密，ode45默认容差下反复试探步长，速度很慢。解决办法是调整tdata的采样间隔，观测60天用1天一采样就够，如果数据是小时级再加密。ODE精度不需要调到很高，参数优化本身不是做数值预报，给迭代次数留更多预算反而更划算。

5.2 边界约束处理不当导致陷入边界解

HHO在开发阶段生成的新位置很容易越过边界。如果只是简单截断到边界值，会有大量个体堆积在边界上，算法把边界当成临时“最优解”，很难真正跳出。所以我用“越界后重新映射到边界内随机位置”的方式，这比截断更有效地维持种群多样性。我自己刚开始做时吃过亏：把越界分量直接设为上下界，结果β频繁跑到0.9，拟合曲线明显过冲，但误差曲线却迟迟不下降。

5.3 怎么判断优化结果是否陷入局部最优

连续跑10次HHO，记录每次的MSE和参数组合。如果10次结果参数抖动很大、误差也忽高忽低，大概率是陷入不同局部最优。此时优先增加种群规模到50到80，迭代次数提高到200以上，而不是盲目扩大边界。另一种诊断方法是看收敛曲线，如果误差曲线在前20代就水平贴着直线，后面完全没有下降，说明探索能力不足，可以适当调大Levy飞行步长系数，或者改变逃逸能量的随机幅度。从经验来看，HHO对SEIR这种平滑度不高的低维问题很少完全早熟，但边界设太宽时局部坑依然不少。

5.4 观测数据口径对拟合结果的直接影响

最后说一个不是代码能解决的问题：如果观测数据只有每日新增确诊数，直接用I曲线拟合会有很大偏差。SEIR里的I是当前具有传染性的人数，而不少统计口径的“确诊”是累计数，需要先变换成现有感染人数，或者直接引入累计确诊的扩展模型。我在模拟数据里特意把I定义成模型状态I，并按时间点采样，就是为了避免口径错位。拿到真实数据的第一件事不是写优化代码，而是把数据口径理清楚，否则再好的算法也白搭。

6. 扩展思路与个人体会

6.1 可扩展的模型变体与优化目标

HHO-SEIR这套框架稍作改动就能扩展到其他模型。把SEIR改成SEIRD，加入死亡舱室D和死亡率δ，目标函数变成同时拟合感染序列和死亡序列；改成SEIQR则加入隔离舱室Q，用多目标加权MSE来拟合。换模型时只需要修改seir_model函数里的微分方程和目标函数里的观测映射，HHO部分基本不动。目标函数也不一定用MSE，可以把峰值时间误差、峰值高度误差加进去，这样优化出的参数对峰值预测更准确。如果还要评估干预措施，可以把β改成随时间变化的分段函数，比如前30天β=β1、后30天β=β2，HHO优化的维度增加几个，流程完全兼容。

6.2 个人体会：算法很重要，但参数解释更重要

HHO不是那种“每篇论文都能涨点”的新算法，但用来做ODE参数反演非常顺手。元启发式算法的价值不在学术炫技，而在实际问题里快速拿到一个足够好的参数解。要承认它的局限：不能给出参数后验分布，不确定性只能靠多次运行近似；性能依赖种群随机性，复现时需要固定随机种子。如果想做更严谨的统计推断，可以在HHO得到初值后再接MCMC采样，这样既减少冷启动时间，又能获得不确定性量化。

我自己的经验是，参数优化的目标不是让曲线“完美贴合”每一个点，而是在保持流行病学解释力的前提下做到拟合。每次跑完HHO，我都会把优化参数带回模型做一次传播推演，看R0、峰值时间、峰值强度是否符合实际情况。如果算法给出来的参数拟合误差很小但R0高到离谱，我会宁可接受误差稍大一点但参数更合理的解。希望这套流程能帮你少走几天弯路。