本文介绍如何对非网格、非均匀、甚至缺失值的二维散点数据（x, y → z）进行双线性模型 z = a·x + b·y + c·x·y + d 的最小二乘拟合，直接求解闭式系数 [a, b, c, d]，无需插值或黑盒模型。
双线性拟合本质上是关于参数 $ a, b, c, d $ 的线性回归问题——尽管模型在变量 $ x $ 和 $ y $ 上呈现乘积项（$ xy $），但它对未知系数仍是线性的。因此，无需复杂优化或迭代算法，可直接通过正规方程（Normal Equation） 求得最小二乘意义下的唯一最优解。
给定 $ N $ 个观测样本 $ (x_i, y_i, z_i) $，我们最小化残差平方和：
$$ 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 $ 分别求偏导并令其为零，得到一组线性方程组：
$$ \begin{bmatrix} \sum x_i^2 & \sum x_i y_i & \sum x_i^2 y_i & \sum x_i \ \sum x_i y_i & \sum y_i^2 & \sum x_i y_i^2 & \sum y_i \ \sum x_i^2 y_i & \sum x_i y_i^2 & \sum x_i^2 y_i^2 & \sum x_i y_i \ \sum x_i & \sum y_i & \sum x_i y_i & N \end{bmatrix} \begin{bmatrix} a \ b \ c \ d \end{bmatrix}
\begin{bmatrix} \sum x_i z_i \ \sum y_i z_i \ \sum x_i y_i z_i \ \sum z_i \end{bmatrix} $$
该方程组系数矩阵对称正定（只要数据足够分散），可用 numpy.linalg.solve 稳健求解。
立即学习“Python免费学习笔记（深入）”；
以下为完整、可直接运行的 Python 实现：
import numpy as np
def bilinear_fit(data):
"""
对散点数据 (x, y, z) 进行双线性最小二乘拟合：
z = a*x + b*y + c*x*y + d
返回系数 [a, b, c, d]
Parameters:
-----------
data : list of [x, y, z] or ndarray of shape (N, 3)
Returns:
--------
np.ndarray: [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 和右侧向量 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])
# 求解线性系统：A @ [a,b,c,d] = RHS
coeffs = np.linalg.solve(A, RHS)
return coeffs
# 示例数据（来自用户提供的 30 个样本）
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"a = {a:.12f}")
print(f"b = {b:.12f}")
print(f"c = {c:.12f}")
print(f"d = {d:.12f}")
# 验证拟合效果
predictions = a * np.array(D)[:, 0] + b * np.array(D)[:, 1] + c * np.prod(np.array(D)[:, :2], axis=1) + d
errors = predictions - np.array(D)[:, 2]
print("\n{'x':>9} {'y':>9} {'z':>9} {'fit':>9} {'error':>9}")
for i in range(len(D)):
print(f"{D[i][0]:9.3f} {D[i][1]:9.3f} {D[i][2]:9.3f} {predictions[i]:9.3f} {errors[i]:9.3f}")✅ 优势说明：
- ✅ 完全适配任意散点：不要求网格结构、等间距或完整采样；
- ✅ 解析解、无迭代：数值稳定，计算高效（$ O(N) $ 预处理 + $ O(1) $ 矩阵求解）；
- ✅ 可解释性强：直接输出物理/工程意义明确的系数 $ a,b,c,d $；
- ✅ 兼容 SciPy/Scikit-learn 生态：可作为 sklearn.base.BaseEstimator 封装复用。
⚠️ 注意事项：
- 若数据高度共线性（如所有 $ x $ 值几乎相同），矩阵可能接近奇异，建议先做中心化或使用 np.linalg.lstsq（自动处理秩亏）；
- 如需不确定性估计（标准误、置信区间），可进一步计算协方差矩阵 $ \sigma^2 (A^\top A)^{-1} $；
- 对于更高阶交互（如二次项 $ x^2, y^2 $），只需扩展设计矩阵维度，方法论完全一致。
该方案规避了 LinearNDInterpolator 无法导出系数、sklearn.LinearRegression 需手动构造特征矩阵的繁琐步骤，是科学计算与工程建模中双线性关系建模的简洁、可靠首选。
