1.5.11.2. Non linear least squares curve fitting: application to point extraction in topographical lidar data

这个练习的目的是为一些数据拟合一个模型。本教程中使用的数据是激光雷达数据,并在下面的介绍性段落中详细描述。如果您不耐烦,现在想练习,请跳过它,直接进入Loading and visualization

1.5.11.2.1. Introduction

激光雷达系统是分析散射光的性质以测量距离的光学测距仪。它们中的大多数向目标发射短的光脉冲并记录反射的信号。然后处理该信号以提取激光雷达系统和目标之间的距离。

地形激光雷达系统是嵌入在机载平台中的这样的系统。它们测量平台和地球之间的距离,以便传递关于地球地形的信息(更多细节参见[1])。

[1]Mallet,C.和Bretar,F.全波形地形激光雷达:最先进的。ISPRS摄影测量和遥感杂志 64(1),第1-16页,2009年1月http://dx.doi.org/10.1016/j.isprsjprs.2008.09.007

在本教程中,目标是分析激光雷达系统记录的波形[2]这样的信号包含其中心和振幅允许计算命中目标的位置和一些特性的峰值。当激光束的覆盖区在地球表面上大约1m时,该波束可以在双向传播期间(例如地面和树或建筑物的顶部)撞击多个目标。激光束击中的每个目标的贡献的总和然后产生具有多个峰的复合信号,每个峰包含关于一个目标的信息。

从这些数据提取信息的现有技术方法的一种状态是以高斯函数的和来分解它们,其中每个函数表示由激光束击中的目标的贡献。

因此,我们使用scipy.optimize模块将波形拟合为一个或高斯函数的和。

1.5.11.2.2. Loading and visualization

使用以下方法加载第一个波形:

>>> import numpy as np
>>> waveform_1 = np.load('data/waveform_1.npy')

并可视化:

>>> import matplotlib.pyplot as plt
>>> t = np.arange(len(waveform_1))
>>> plt.plot(t, waveform_1)
[<matplotlib.lines.Line2D object at ...>]
>>> plt.show()
../../_images/waveform_1.png

正如你可以注意到的,这个波形是一个单峰的80 bin长度的信号。

1.5.11.2.3. Fitting a waveform with a simple Gaussian model

信号非常简单,并且可以被建模为单个高斯函数和对应于背景噪声的偏移。要使用函数拟合信号,我们必须:

  • 定义模型
  • 提出一个初步解决方案
  • 调用scipy.optimize.leastsq

1.5.11.2.3.1. Model

定义的高斯函数

B + A \exp\left\{-\left(\frac{t-\mu}{\sigma}\right)^2\right\}

可以在python中定义:

>>> def model(t, coeffs):
... return coeffs[0] + coeffs[1] * np.exp( - ((t-coeffs[2])/coeffs[3])**2 )

其中

  • coeffs[0] isB(噪声)
  • coeffs[1] isA(振幅)
  • coeffs[2] is\mu(中央)
  • coeffs[3] is\sigma(宽度)

1.5.11.2.3.2. Initial solution

从查看图形可以找到的近似初始解决方案是例如:

>>> x0 = np.array([3, 30, 15, 1], dtype=float)

1.5.11.2.3.3. Fit

scipy.optimize.leastsq最小化作为参数给出的函数的平方和。基本上,最小化的函数是残差(数据和模型之间的差):

>>> def residuals(coeffs, y, t):
... return y - model(t, coeffs)

所以让我们通过调用scipy.optimize.leastsq()来获得我们的解决方案:

  • 函数最小化
  • 初始解决方案
  • 传递给函数的附加参数
>>> from scipy.optimize import leastsq
>>> x, flag = leastsq(residuals, x0, args=(waveform_1, t))
>>> print(x)
[ 2.70363341 27.82020742 15.47924562 3.05636228]

并可视化解决方案:

>>> plt.plot(t, waveform_1, t, model(t, x)) 
[<matplotlib.lines.Line2D object at ...>, <matplotlib.lines.Line2D object at ...>]
>>> plt.legend(['waveform', 'model'])
<matplotlib.legend.Legend object at ...>
>>> plt.show()

注意:从scipy v0.8及以上版本,你应该使用scipy.optimize.curve_fit(),它将模型和数据作为参数,因此不需要定义残差了。

1.5.11.2.4. Going further

  • 尝试使用包含三个显着峰的更复杂的波形(例如data/waveform_2.npy)。您必须调整现在是高斯函数之和而不是仅一个高斯峰的模型。

    ../../_images/waveform_2.png
  • 在某些情况下,编写明确的函数来计算雅可比,比使用leastsq数值估计更快。创建一个函数来计算残差的雅可比,并将其用作leastsq的输入。

  • 当我们想要在信号中检测非常小的峰值时,或者当初始猜测距离好的解决方案太远时,算法给出的结果通常是不令人满意的。向模型的参数添加约束使得能够克服这种限制。我们可以添加的先验示例是我们的变量的符号(都是正数)。

    使用以下初始解决方案:

    >>> x0 = np.array([3, 50, 20, 1], dtype=float)
    

    比较scipy.optimize.leastsq()的结果以及在添加边界约束时可以使用scipy.optimize.fmin_slsqp()获得的结果。

[2]本教程中使用的数据是FullAnalyze软件可用的演示数据的一部分,由GIS DRAIX友情提供。