基于双TIFF栅格影像的四阶多项式系数优化线性回归咨询

一、解决栅格尺寸不匹配问题 #

两幅影像宽度一致（608列），仅高度差1行（232 vs 231），可通过以下两种方式对齐：

• 直接裁剪对齐：将行数多的影像裁剪为231行，与另一幅尺寸匹配。使用GDAL命令行：

  gdal_translate -outsize 608 231 input_232.tif output_cropped.tif
  

  或GDAL Python API实现：

  from osgeo import gdal
  
  src_ds = gdal.Open("input_232.tif")
  dst_ds = gdal.GetDriverByName("GTiff").Create(
      "output_cropped.tif", 608, 231, 
      src_ds.RasterCount, src_ds.GetRasterBand(1).DataType
  )
  dst_ds.SetGeoTransform(src_ds.GetGeoTransform())
  dst_ds.SetProjection(src_ds.GetProjection())
  
  for band_idx in range(src_ds.RasterCount):
      src_band = src_ds.GetRasterBand(band_idx + 1)
      dst_band = dst_ds.GetRasterBand(band_idx + 1)
      data = src_band.ReadAsArray(0, 0, 608, 231)
      dst_band.WriteArray(data)
  
  src_ds = None
  dst_ds = None
  

• 重采样对齐：若尺寸差异源于地理范围不一致，先确认两幅影像坐标系/投影一致，再用gdalwarp将其中一幅重采样到目标尺寸：

  gdalwarp -ts 608 231 -r bilinear input.tif output_warped.tif
  

二、处理NA值并优化系数 #

1. 提取有效像素对 #

读取对齐后的影像，筛选出x、y均不为NA的像素：

import numpy as np
from osgeo import gdal

def read_raster(path):
    ds = gdal.Open(path)
    band = ds.GetRasterBand(1)
    no_data = band.GetNoDataValue()
    data = band.ReadAsArray()
    ds = None
    return data, no_data

x_data, x_no_data = read_raster("aligned_x.tif")
y_data, y_no_data = read_raster("aligned_y.tif")

# 生成有效像素掩码
valid_mask = (x_data != x_no_data) & (y_data != y_no_data)
x_valid = x_data[valid_mask].flatten()
y_valid = y_data[valid_mask].flatten()

2. 四次多项式系数拟合 #

使用最小二乘法求解系数，可直接调用numpy的polyfit：

# 拟合四次多项式，返回系数顺序为[a5, a4, a3, a2, a1]
coeffs = np.polyfit(x_valid, y_valid, 4)
# 转换为方程所需顺序：a1, a2, a3, a4, a5
a1, a2, a3, a4, a5 = coeffs[::-1]

或手动构建设计矩阵求解（适合自定义约束）：

# 构建设计矩阵：每行对应[1, x, x², x³, x⁴]
X = np.vstack([
    np.ones_like(x_valid), x_valid, 
    x_valid**2, x_valid**3, x_valid**4
]).T
# 最小二乘法求解
coeffs, _, _, _ = np.linalg.lstsq(X, y_valid, rcond=None)
a1, a2, a3, a4, a5 = coeffs

3. 拟合效果验证 #

可通过计算R²评估拟合质量：

y_pred = a1 + a2*x_valid + a3*x_valid**2 + a4*x_valid**3 + a5*x_valid**4
r_squared = 1 - np.sum((y_valid - y_pred)**2) / np.sum((y_valid - np.mean(y_valid))**2)

内容的提问来源于stack exchange，提问作者Aurélien Lengrand