本文介绍如何对非网格、不规则分布的二维输入（x, y）与输出（z）数据，使用最小二乘法拟合双线性模型 z = a·x + b·y + c·x·y + d，并准确求解系数 a、b、c、d。方法纯 numpy 实现，无需网格化或插值，适用于缺失值、不等采样密度的真实实验数据。
双线性拟合本质上是带交叉项的线性回归：虽然模型在 x 和 y 上是非线性的（含乘积项 x·y），但关于待估参数 a、b、c、d 是线性的。因此，可直接套用多元线性最小二乘法——即最小化残差平方和
$$ S = sum_{i=1}^N left( a x_i + b y_i + c x_i y_i + d – z_i right)^2 $$
对 a、b、c、d 分别求偏导并令其为 0，可导出标准的正规方程组（Normal Equations）： $$ mathbf{A} cdot begin{bmatrix}a b c dend{bmatrix} = mathbf{b} $$ 其中设计矩阵 $mathbf{A}$ 和右侧向量 $mathbf{b}$ 的元素由数据的各阶统计量构成（如 $sum x_i$, $sum x_i y_i$, $sum x_i y_i z_i$ 等）。
以下是一个高效、可复用的 bilinear_fit 函数实现，完全基于 NumPy，无外部依赖：
import numpy as np def bilinear_fit(data): """ 对散点数据 (x, y, z) 进行双线性最小二乘拟合：z = a*x + b*y + c*x*y + d Parameters: ----------- data : list of [x, y, z] or ndarray of shape (N, 3) 输入数据点集合，支持任意数量、非均匀、非网格分布 Returns: -------- tuple: (a, b, c, d) 拟合系数 """ data = np.asarray(data) x, y, z = data[:, 0], data[:, 1], data[:, 2] N = len(x) # 计算必要统计量（避免循环，更高效） Sx = np.sum(x) Sxx = np.sum(x * x) Sy = np.sum(y) Syy = np.sum(y * y) Sxy = np.sum(x * y) Sxxy = np.sum(x * x * y) Sxyy = np.sum(x * y * y) Sxxyy = np.sum(x * x * y * y) Sz = np.sum(z) Sxz = np.sum(x * z) Syz = np.sum(y * z) Sxyz = np.sum(x * y * z) # 构建正规方程组 A @ [a,b,c,d] = RHS A = np.array([ [Sxx, Sxy, Sxxy, Sx ], [Sxy, Syy, Sxyy, Sy ], [Sxxy, Sxyy, Sxxyy, Sxy], [Sx, Sy, Sxy, N ] ]) RHS = np.array([Sxz, Syz, Sxyz, Sz]) try: coeffs = np.linalg.solve(A, RHS) return coeffs[0], coeffs[1], coeffs[2], coeffs[3] except np.linalg.LinAlgError as e: raise ValueError("正规方程组矩阵奇异，请检查数据是否共线（如 x 或 y 全为常数）或样本量不足（N < 4）") from e # 示例：使用用户提供的真实数据 D = [ [1056, 8, 50.89124679], [1056, 16, 61.62827273], [1056, 32, 78.83079982], [1056, 48, 92.90073197], [1056, 64, 105.103744 ], [1056, 80, 116.0303753], [1056, 96, 126.0130906], [1056, 112, 135.2610439], [1056, 128, 143.9159512], [1056, 144, 152.0790946], [2048, 8, 63.71675604], [2048, 16, 77.15971099], [2048, 32, 98.69757849], [2048, 48, 116.313387 ], [2048, 64, 131.5917779], [2048, 80, 145.2721136], [2048, 96, 157.7706532], [2048, 112, 169.3492575], [2048, 128, 180.1853546], [2048, 144, 190.4057615], [4096, 8, 86.7357654 ], [4096, 16, 105.0352703], [4096, 32, 134.3541477], [4096, 48, 158.334033 ], [4096, 64, 179.1320602], [4096, 80, 197.7547066], [4096, 96, 214.7686034], [4096, 112, 230.5302193], [4096, 128, 245.2810877], [4096, 144, 259.193829 ] ] a, b, c, d = bilinear_fit(D) print(f"拟合结果：") print(f"a (x 系数) = {a:.10f}") print(f"b (y 系数) = {b:.10f}") print(f"c (xy 交叉项) = {c:.10f}") print(f"d (截距) = {d:.10f}") # 验证拟合效果 z_pred = a * np.array(D)[:,0] + b * np.array(D)[:,1] + c * np.array(D)[:,0] * np.array(D)[:,1] + d residuals = z_pred - np.array(D)[:,2] print(f"nR² 决定系数 = {1 - np.sum(residuals**2) / np.sum((np.array(D)[:,2] - np.mean(D, axis=0)[2])**2):.4f}")✅ 关键优势说明：
- ✅ 完全适配散点数据：不要求 x-y 构成矩形网格，允许任意分布、缺失点、重复 x/y 值；
- ✅ 解析解、高精度：通过正规方程直接求解，比迭代法更稳定、更快；
- ✅ 可解释性强：返回明确物理意义的系数 a, b, c, d，便于后续分析或嵌入公式；
- ✅ 轻量可靠：仅依赖 NumPy，无 sklearn/scipy 插值器的黑盒限制（如无法导出系数、边界外推异常等）。
⚠️ 注意事项：
- 数据量应 ≥ 4（参数个数），否则方程组欠定；
- 若 x 或 y 列全为常数（如所有 y=8），会导致设计矩阵秩亏，需预处理或改用正则化（如 np.linalg.lstsq 加 rcond）；
- 对强非线性或存在显著离群点的数据，建议先做标准化（如 (x - x_mean)/x_std）以提升数值稳定性；
- 如需不确定性评估（如系数标准误），可进一步计算协方差矩阵 $(mathbf{X}^top mathbf{X})^{-1} sigma^2$，其中 $mathbf{X}$ 为设计矩阵 $[x_i,, y_i,, x_i y_i,, 1]$。
该方法已在多个工程场景（如传感器响应建模、GPU性能预测、图像几何校正）中验证有效，是处理双线性关系最直接、透明且可控的解决方案。
立即学习“Python免费学习笔记（深入）”；