本文介绍如何对非网格、非均匀、甚至缺失值的二维散点数据（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 实现：

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}")