4. 响应面模型(多项式回归)

响应面模型的核心思想,是用一个显式的多项式函数去逼近输入变量与响应输出之间的真实关系。它不依赖复杂的物理方程,而是通过有限次试验或仿真采样,把“黑箱”关系拟合成可解析、可求导、可优化的代理模型。本节围绕一阶与二阶多项式响应面的数学形式、最小二乘拟合、系数显著性检验、模型诊断与预测区间估计展开,并用一个工程案例完成完整拟合过程。

一、一阶与二阶多项式响应面的数学形式

设输入变量为 x₁, x₂, …, x_k,响应为 y。一阶响应面(线性模型)假设响应随各变量线性变化,其形式为:

y = β₀ + β₁x₁ + β₂x₂ + … + β_k x_k + ε

其中 β₀ 为截距,β_i 为第 i 个变量的线性效应系数,ε 为随机误差。一阶模型适用于变量影响近似线性、交互作用较弱的场景,常用于筛选试验后的初步建模。

二阶响应面(二次模型)在一阶基础上加入平方项和两两交互项,能够刻画曲率和变量间的耦合效应:

y = β₀ + Σβ_i x_i + Σβ_ii x_i² + ΣΣβ_ij x_i x_j + ε   (i < j)

以两个变量为例,完整展开为:

y = β₀ + β₁x₁ + β₂x₂ + β₁₁x₁² + β₂₂x₂² + β₁₂x₁x₂ + ε

二阶模型是响应面法中最常用的形式,因为它能描述极值点、鞍点和脊线,适合后续的优化与寻优。当变量数 k 较大时,二阶模型待估系数个数为 (k+1)(k+2)/2,例如 k=3 时需估计 10 个系数,因此试验点数量必须足够。

二、最小二乘拟合

把 n 次试验的输入与响应写成矩阵形式:设计矩阵 X(n×(p+1),p 为变量项数)、响应向量 y、系数向量 β。最小二乘的目标是最小化残差平方和:

S(β) = (y − Xβ)ᵀ(y − Xβ)

对 β 求导并令其为零,得到正规方程 XᵀXβ = Xᵀy,其解为:

β̂ = (XᵀX)⁻¹Xᵀy

实际计算中通常不直接求逆,而是用 QR 分解或奇异值分解求解,以提高数值稳定性。拟合完成后,预测值为 ŷ = Xβ̂,残差为 e = y − ŷ。最小二乘不要求误差服从正态分布即可给出系数估计,但要做区间估计和显著性检验,则需假设 ε 独立同分布且服从 N(0, σ²)。

三、系数显著性检验

拟合出系数后,需要判断哪些项真正对响应有影响。对每个系数做 t 检验,原假设 H₀: β_j = 0。检验统计量为:

t_j = β̂_j / se(β̂_j)

其中 se(β̂_j) = σ̂·√[(XᵀX)⁻¹]_jj,σ̂² = SSE/(n − p − 1) 为误差方差估计。若 |t_j| 超过 t 分布临界值,或 p 值小于显著性水平(常取 0.05),则拒绝原假设,认为该项显著。对模型整体,用 F 检验判断回归是否显著:F = (SSR/p)/(SSE/(n−p−1))。工程实践中,常按 p 值从大到小逐步剔除不显著项,但需注意层级原则——若交互项显著,其主效应通常应保留。

四、模型诊断与预测区间估计

模型诊断主要检查残差是否满足假设。常用手段包括:残差正态概率图(点应近似落在直线上)、残差与预测值散点图(应无喇叭形或弯曲趋势)、残差与各变量散点图(应无结构)。若发现异方差或非线性,可考虑变换响应或增加项。决定系数 R² 和调整 R² 用于衡量拟合优度,R² 接近 1 说明模型解释了大部分变异,但过高且伴随大量不显著项时需警惕过拟合。

预测区间估计给出新输入点 x₀ 处单次观测的区间:

ŷ₀ ± t_(α/2, n−p−1)·σ̂·√(1 + x₀ᵀ(XᵀX)⁻¹x₀)

若估计的是响应均值区间,则去掉根号内的 1。预测区间比均值区间宽,因为它同时包含模型不确定性和单次观测的随机误差。区间越窄,说明模型在该点的预测越可靠。

二阶响应面拟合与诊断流程 试验数据 x₁, x₂, y(n 组) 构造设计矩阵 X 最小二乘拟合 β̂ = (XᵀX)⁻¹Xᵀy QR / SVD 求解 系数显著性检验 t 检验 / F 检验 剔除不显著项 模型诊断 残差正态图 R² / 调整 R² 诊断不通过:变换响应 / 增加项 / 补点 预测区间估计 ŷ₀ ± t·σ̂·√(1 + x₀ᵀ(XᵀX)⁻¹x₀) 均值区间去掉根号内 1,预测区间更宽 工程案例:某换热器翅片间距 x₁(mm) 与风速 x₂(m/s) 对换热量 y(W) 的影响(示例数据) 试验号 x₁ x₂ y 拟合项 系数估计 p 值 1~9 2~6 1~5 观测值 β₀ / β₁ / β₂ 截距、线性项 <0.05 显著 中心点 4 3 重复 3 次 β₁₁ / β₂₂ / β₁₂ 平方、交互项 判断曲率与耦合

五、工程案例:翅片换热量的二阶响应面拟合

某换热器设计中,关注翅片间距 x₁(2~6 mm)和迎面风速 x₂(1~5 m/s)对换热量 y(W)的影响。采用中心复合设计安排 9 组试验,中心点(x₁=4,x₂=3)重复 3 次以估计纯误差。以下数据为教学示例,用于演示拟合流程,并非实测结果。

试验号x₁ (mm)x₂ (m/s)y (W)
121182
261205
325268
465301
523221
663248
741196
845289
9~1143243, 246, 240

将变量编码为中心化形式:u₁=(x₁−4)/2,u₂=(x₂−3)/2,使各因子取值在 [−1, 1] 内。构造二阶模型:

y = β₀ + β₁u₁ + β₂u₂ + β₁₁u₁² + β₂₂u₂² + β₁₂u₁u₂ + ε

用最小二乘求解,得到一组示例系数估计:β̂₀≈243.0,β̂₁≈12.5,β̂₂≈38.2,β̂₁₁≈−4.8,β̂₂₂≈−6.1,β̂₁₂≈2.5。对系数做 t 检验,若取 α=0.05,β₁、β₂ 的 p 值远小于 0.05,说明翅片间距和风速的线性效应显著;β₂₂ 的 p 值约 0.03,说明风速存在明显二次曲率;β₁₂ 的 p 值约 0.21,交互项不显著,可考虑剔除后重新拟合。剔除交互项后,模型简化为含两个平方项的二次模型,调整 R² 从 0.972 变为 0.975,说明简化未损失拟合能力。

诊断阶段,绘制残差正态概率图,点近似沿直线分布;残差与预测值散点图无明显喇叭形,说明等方差假设基本成立。在 x₁=4 mm、x₂=3 m/s 处,预测换热量约 243 W,单次预测区间约为 243±18 W(示例值),均值区间约为 243±7 W。可见预测区间明显宽于均值区间,工程上若用于单台设备换热量的预估,应使用预测区间而非均值区间。

建模提示:二阶模型至少需要 (k+1)(k+2)/2 个不同试验点,且中心点重复 3~5 次,才能同时估计纯误差和检验失拟。若设计点不足,XᵀX 接近奇异,系数估计会极不稳定。
注意:R² 高不等于模型可用。若残差图出现弯曲趋势,说明二阶项仍不足以描述真实关系,此时应检查变量范围是否过大,或考虑对数变换、Box-Cox 变换后再拟合。

综上,一阶模型用于快速判断主效应,二阶模型用于刻画曲率与交互并支撑优化;最小二乘给出系数估计,t 检验与 F 检验判断显著性,残差诊断验证假设,预测区间量化单点预测的不确定性。这五步构成响应面建模的完整闭环,也是后续寻优与稳健设计的基础。