我正在使用 python/numpy/scipy 来实现此算法，以根据地形坡向和坡度对齐两个数字高程模型 (DEM):

“用于量化冰川厚度变化的卫星高程数据集的配准和偏差校正”，C. Nuth 和 A. Kääb，doi:10.5194/tc-5-271-2011

我已经设置了一个框架，但是 scipy.optimize.curve_fit 提供的拟合质量很差。

def f(x, a, b, c):
    y = a * numpy.cos(numpy.deg2rad(b-x)) + c
    return y

def compute_offset(dh, slope, aspect):
    import scipy.optimize as optimization

    idx = random.sample(range(dh.compressed().size), 10000)
    xdata = numpy.array(aspect.compressed()[idx], float)
    ydata = numpy.array((dh/numpy.tan(numpy.deg2rad(slope))).compressed()[idx], float)

    #Generate synthetic data to test curve_fit
    #xdata = numpy.arange(0,360,0.01)
    #ydata = f(xdata, 20.0, 130.0, -3.0) + 20*numpy.random.normal(size=len(xdata))

    print xdata
    print ydata

    x0 = numpy.array([0.0, 0.0, 0.0])

    fit = optimization.curve_fit(f, xdata, ydata, x0)[0]
    #optimization.leastsq(f, x0[:], args=(xdata, ydata))
    genplot(xdata, ydata, fit)
    return fit

def genplot(x, y, fit):
    a = (numpy.arange(0,360))
    f_a = f(a, fit[0], fit[1], fit[2])
    idx = random.sample(range(x.size), 10000)
    plt.figure()
    plt.xlabel('Aspect (deg)')
    plt.ylabel('dh/tan(slope) (m)')
    plt.plot(x[idx], y[idx], 'r.')
    plt.axhline(color='k')
    plt.plot(a, f_a, 'b')
    plt.ylim(-80,80)
    plt.show()

#Input DEMs
dem1_fn = sys.argv[1]
dem2_fn = sys.argv[2]

dem1_ds = gdal.Open(dem1_fn, gdal.GA_ReadOnly)
dem2_ds = gdal.Open(dem2_fn, gdal.GA_ReadOnly)

#Extract band 1 from each dataset as masked array using internal nodata value
dem1 = getperc_new.gdal_getma(dem1_ds, 1)
dem2 = getperc_new.gdal_getma(dem2_ds, 1)

#Produce slope and aspect maps using gdaldem and load into masked arrays
dem1_slope = gdaldem_slope(dem1_fn)
dem1_aspect = gdaldem_aspect(dem1_fn)

#Compute common mask and apply to all products
common_mask = dem1.mask + dem2.mask + dem1_slope.mask + dem1_aspect.mask

diff_euler = numpy.ma.array(dem2-dem1, mask=common_mask)
dem1_slope.__setmask__(common_mask)
dem1_aspect.__setmask__(common_mask)

#Compute relationship between elevation difference, slope and aspect
fit = compute_offset(diff_euler, dem1_slope, dem1_aspect)
print fit

这是适合我的数据，最初由大约 200 万个点组成，但我随机抽样用于测试/绘图目的:

[-14.9639559 216.01093596 -41.96806735]

那里有很多数据可以很好地拟合，但 curve_fit 的结果很差。当我使用合成数据运行时，我得到了一个很好的拟合:

原始输入参数[20.0, 130.0, -3.0]

曲线拟合结果 [-19.66719631 -49.6673076 -3.12198723]

不确定这是否与使用掩码数组、curve_fit 的限制有关，或者我是否只是忽略了一些简单的事情。感谢您的任何建议。

==========================

编辑 9/4/13 16:30 PDT

正如@Evert 和其他人所建议的，问题肯定与异常值有关。去除异常值后，我能够获得更好的拟合度。查看我的旧代码，似乎我计算了每个方面的中值绝对偏差，然后在拟合之前删除了 2*mad 之外的任何内容。

我在 2012 年 11 月生成了一些额外的图:

但再次查看这些，我几乎可以肯定它们是针对不同的输入数据生成的。这是我现在能找到的所有内容，所以我将它们包括在这里作为有偏差抽样的案例示例。对于此类情况，这种 DEM 对齐方法必然会失败 - 它与 scipy 的曲线拟合能力无关。

我最终开发了一种不同的对齐方法，涉及两个蒙版 2D numpy 阵列的归一化互相关、子像素细化和垂直偏移移除。它更快并且始终如一地提供更好的结果。尽管这种方法已被 Oleg Alexandrov 作为 NASA Ames Stereo Pipeline 的一部分开发的迭代最近点 (ICP) 工具 (pc_align) 所取代。 .

感谢您的所有回复，对于放弃这个问题，我深表歉意。

最佳答案

如果您只是想获得具有相位偏移的正弦波，则不需要非线性拟合。

您可以将 a * sin(x - b) + c 替换为 a * sin(x) + b * cos(x) + c，因为任何带有偏移量的正弦可以写成正弦和余弦的适当组合(“相量加法”，就像在傅里叶变换中一样)。

如果这给出了相同的结果，那么问题就不是“非线性”拟合。

关于python - 我知道 scipy curve_fit 可以做得更好，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/12397412/