我使用 python/numpy/scipy 来实现此算法，用于根据地形坡向和坡度对齐两个数字高程模型 (DEM)：
“用于量化冰川厚度变化的卫星高程数据集的联合配准和偏差校正”，C. Nuth 和 A. Kääb，doi：10.5194/tc-5-271-2011 http://dx.doi.org/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]
curve_fit 的结果 [-19.66719631 -49.6673076 -3.12198723]
不确定这是否与使用屏蔽数组、 curve_fit 的限制有关，或者我是否只是忽略了一些简单的事情。感谢您的任何建议。
==========================
编辑 2013 年 9 月 4 日 16:30（太平洋夏令时间）
正如@Evert 和其他人所建议的，问题肯定与异常值有关。去除异常值后，我能够获得更好的拟合效果。看看我的旧代码，似乎我计算了每个方面的中值绝对偏差，然后在拟合之前删除了 2*mad 之外的任何内容。
我在 2012 年 11 月生成了一些额外的图：
但再次查看这些，我几乎肯定它们是针对不同的输入数据生成的。这是我现在能找到的所有内容，因此我将它们作为有偏差抽样的案例示例放在这里。对于这样的情况，这种 DEM 对齐方法注定会失败 - 并且它与 scipy 的曲线拟合能力无关。
我最终开发了一种不同的对齐方法，涉及两个屏蔽 2D numpy 数组的归一化互相关、子像素细化和垂直偏移去除。它速度更快，并且始终提供更好的结果。尽管这种方法已经被 Oleg Alexandrov 开发的迭代最近点 (ICP) 工具 (pc_align) 所取代，作为NASA 艾姆斯立体声管道 http://ti.arc.nasa.gov/tech/asr/intelligent-robotics/ngt/stereo/.
感谢您的所有回复，对于放弃这个问题我深表歉意。
