9

我正在使用 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]

plot_one

那里有大量数据可以很好地拟合，但curve_fit的结果很差。当我使用合成数据运行时，我很适合：

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

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

plot_two

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

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

编辑 2013 年 9 月 4 日 16:30 PDT

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

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

直方图

在此处输入图像描述

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

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

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

4

1 回答 1

1

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

您可以将其替换a * sin(x - b) + c为a * sin(x) + b * cos(x) + c，因为任何带有偏移的正弦都可以写为正弦和余弦的适当组合（“相量加法”，如傅立叶变换）。

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

于 2012-11-03T16:52:00.493 回答