1.5.11 科学计算的总结练习

总结练习主要使用Numpy、Scipy 和 Matplotlib。他们提供了一些使用Python进行科学计算的真实例子。现在，已经介绍了Numpy和Scipy的基本使用，邀请感兴趣的用户去做这些练习。

练习：

1.5.11.13 Sprogø气象站的最大风速预测

1.5.11.14 非线性最小二乘曲线拟合：地形机载激光雷达数据中的点抽取

1.5.11.15 图像处理应用：计数气泡和未融化的颗粒

提议的解决方案:

1.5.11.16 图像处理练习：玻璃中的未融化颗粒的答案例子

1.5.11.13 Sprogø气象站的最大风速预测

这个练习的目的是预测每50年的最大风速，即使在一个时间段内有记录。可用的数据只是位于丹麦的Sprogø气象站的21年的测量数据。首先，将给出统计步骤，接着将用scipy.interpolae模块中的函数来解释。在最后，将邀请感兴趣的读者用不同的方法从原始数据计算结果。

1.5.11.13.1 统计方法

假设年度最大值符合正态概率密度函数。但是，这个函数不能用来预测，因为它从速度最大值中给出了概率。找到每50年的最大风速需要相反的方法，需要从确定的概率中找到结果。这是百分位数函数的作用而这个练习的目的是找到它。在当前的模型中，假设每50年出现的最大风速定义为高于2%百分位数。

根据定义，百分位数函数是累积分布函数的反函数。后者描述了年度最大值的概率分布。在这个练习中，给定年份$i$的累积概率$p_i$被定义为$p_i = i/(N+1)$，其中$N = 21$，测量的年数。因此，计算每个测量过的风速最大值的累积概率是可以行的。从这些实验点，scipy.interpolate模块将对拟合百分位数函数非常有用。最后，50年的最大值将从累积概率的2%百分位数中预估出来。

1.5.11.13.2 计算累积概率

计算好的numpy格式的年度风速最大值存储在examples/max-speeds.npy文件中, 因此，可以用numpy加载：

In [4]:

  1. import numpy as np
  2. max_speeds = np.load('data/max-speeds.npy')
  3. years_nb = max_speeds.shape[0]

下面是前面板块的累积概率定义$p_i$，对应值将为：

In [5]:

  1. cprob = (np.arange(years_nb, dtype=np.float32) + 1)/(years_nb + 1)

并且假设他们可以拟合给定的风速：

In [6]:

  1. sorted_max_speeds = np.sort(max_speeds)

1.5.11.13.3 用UnivariateSpline预测

在这个部分，百分位数函数将用UnivariateSpline类来估计，这个类用点代表样条。 默认行为是构建一个3度的样条，不同的点根据他们的可靠性可能有不同的权重。相关的变体还有InterpolatedUnivariateSpline和LSQUnivariateSpline，差别在于检查误差的方式不同。如果需要2D样条，可以使用BivariateSpline家族类。所有这些1D和2D样条使用FITPACK Fortran 程序，这就是为什么通过splrep和splev函数来表征和评估样条的库更少。同时，不使用FITPACK参数的插值函数也提供更简便的用法（见interp1d, interp2d, barycentric_interpolate等等）。 对于Sprogø最大风速的例子，将使用UnivariateSpline，因为3度的样条似乎可以正确拟合数据：

In [7]:

  1. from scipy.interpolate import UnivariateSpline
  2. quantile_func = UnivariateSpline(cprob, sorted_max_speeds)

百分位数函数将用评估来所有范围的概率：

In [8]:

  1. nprob = np.linspace(0, 1, 1e2)
  2. fitted_max_speeds = quantile_func(nprob)

在当前的模型中，每50年出现的最大风速被定义为大于2%百分位数。作为结果，累积概率值将是：

In [9]:

  1. fifty_prob = 1. - 0.02

因此，可以猜测50年一遇的暴风雪风速为：

In [10]:

  1. fifty_wind = quantile_func(fifty_prob)
  2. fifty_wind

Out[10]:

  1. array(32.97989825386221)

现在，结果被收集在Matplotlib图片中：

1.5.11 科学计算的总结练习 - 图1

答案：Python源文件

1.5.11.13.4 Gumbell分布练习

现在邀请感兴趣的读者用21年测量的风速做一个练习。测量区间为90分钟（原始的区间约为10分钟，但是，为了让练习的设置简单一些，缩小了文件的大小）。数据以numpy格式存储在文件examples/sprog-windspeeds.npy中。 在完成练习后，不要看绘图的源代码。

  • 第一步将是通过使用numpy来找到年度最大值，然后将它们绘制为matplotlibe条形图。 1.5.11 科学计算的总结练习 - 图2

答案：Python源文件

  • 第二步将是在累积概率$p_i$使用Gumbell分布，$p_i$的定义是$-log( -log(p_i) )$用来拟合线性百分位数函数（记住你可以定义UnivariateSpline的度数）。 绘制年度最大值和Gumbell分布将生产如下图片。

1.5.11 科学计算的总结练习 - 图3

答案：Python源文件

  • 最后一步将是找到在每50年出现的最大风速34.23 m/s。

1.5.11.14 非线性最小二乘曲线拟合：地理雷达数据中的点抽取应用

这个练习的目的是用模型去拟合一些数据。这篇教程中的数据是雷达数据，下面的介绍段落将详细介绍。如果你没有耐心，想要马上进行联系，那么请跳过这部分，并直接进入加载和可视化。

1.5.11.14.1 介绍

雷达系统是光学测距仪，通过分析离散光的属性来测量距离。绝大多数光学测距仪向目标发射一段短光学脉冲，然后记录反射信号。然后处理这个信号来抽取雷达系统与目标间的距离。

地形雷达系统是嵌入在飞行平台的雷达系统。它们测量平台与地球的距离，以便计算出地球的地形信息（更多细节见[1]）。

[1] Mallet, C. and Bretar, F. Full-Waveform Topographic Lidar: State-of-the-Art. ISPRS Journal of Photogrammetry and Remote Sensing 64(1), pp.1-16, January 2009 http://dx.doi.org/10.1016/j.isprsjprs.2008.09.007

这篇教程的目的是分析雷达系统记录到的波形数据[2]。这种信号包含波峰，波峰的中心和振幅可以用来计算命中目标的位置和一些特性。当激光柱的脚步距离地球表面1m左右，光柱可以在二次传播时击中多个目标（例如，地面和树木或建筑的顶部）。激光柱的击中每个目标的贡献之和会产生一个有多个波峰的复杂波，每一个包含一个目标的信息。

一种从这些数据中抽取信息的先进方法是在一个高斯函数和中分解这些信息，每个函数代表激光柱击中的一个目标的贡献。

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

1.5.11.14.2 加载和可视化

加载第一个波形：

In [1]:

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

接着可视化：

In [2]:

  1. import matplotlib.pyplot as plt
  2. t = np.arange(len(waveform_1))
  3. plt.plot(t, waveform_1)
  4. plt.show()

1.5.11 科学计算的总结练习 - 图4

你可以注意到，这个波形是单峰80个区间的信息。

1.5.11.14.3 用简单的高斯模型拟合波形

这个信号非常简单，可以被建模为一个高斯函数，抵消相应的背景噪音。要用函数拟合这个信号，我们必须：

  • 定义一个模型
  • 给出初始解
  • 调用scipy.optimize.leastsq
1.5.11.14.3.1 模型

高斯函数定义如下：

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

在Python中定义如下：

In [3]:

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

其中

  • coeffs[0] is $B$ (noise)
  • coeffs[1] is $A$ (amplitude)
  • coeffs[2] is $\mu$ (center)
  • coeffs[3] is $\sigma$ (width)
1.5.11.14.3.2 初始解

通过观察图形，我们可以找到大概的初始解，例如：

In [5]:

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

scipy.optimize.leastsq最小化作为参数给到的函数的平方和。本质上来说，函数最小化的是残差（数据与模型的差异）：

In [6]:

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

因此，让我们通过下列参数调用scipy.optimize.leastsq来求解：

  • 最小化的函数
  • 初始解
  • 传递给函数的额外参数

In [7]:

  1. from scipy.optimize import leastsq
  2. x, flag = leastsq(residuals, x0, args=(waveform_1, t))
  3. print x
  1. [  2.70363341  27.82020741  15.47924562   3.05636228]

答案可视化：

In [8]:

  1. plt.plot(t, waveform_1, t, model(t, x))
  2. plt.legend(['waveform', 'model'])
  3. plt.show()

1.5.11 科学计算的总结练习 - 图5

备注：从scipy v0.8及以上，你应该使用scipy.optimize.curve_fit，它使用模型和数据作为参数，因此，你不再需要定义残差。

1.5.11.14.4 更进一步

  • 试一下包含三个波峰的更复杂波形（例如data/waveform_2.npy）。你必须调整模型，现在它是高斯函数之和，而不是只有一个高斯波峰。

1.5.11 科学计算的总结练习 - 图6

  • 在一些情况下，写一个函数来计算Jacobian，要比让leastsq从数值上估计它来的快。创建一个函数来计算残差的Jacobian，并且用它作为leastsq的一个输入。
  • 当我们想要识别信号中非常小的峰值，或者初始的猜测离好的解决方案太远时，算法给出的结果往往不能令人满意。为模型参数添加限制可以确保克服这些局限性。我们可以添加的先前经验是变量的符号（都是正的）。

用下列初始解：

In [9]:

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

添加了边界限制之后比较一下scipy.optimize.leastsq与scipy.optimize.fmin_slsqp的结果。

[2] 本教程的数据部分来自于FullAnalyze software的演示数据，由 GIS DRAIX 友情提供。

1.5.11.15 图像处理应用：计数气泡和未融化的颗粒

1.5.11 科学计算的总结练习 - 图7

1.5.11.15.1 问题描述

  • 打开图像文件MV_HFV_012.jpg并且浏览一下。看一下imshow文档字符串中的参数，用“右”对齐来显示图片（原点在左下角，而不是像标准数组在右上角）。

这个扫描元素显微图显示了一个带有一些气泡（黑色）和未溶解沙（深灰）的玻璃样本（轻灰矩阵）。我们想要判断样本由三个状态覆盖的百分比，并且预测沙粒和气泡的典型大小和他们的大小等。

  • 修建图片，删除带有测量信息中底部面板。

  • 用中位数过滤稍稍过滤一下图像以便改进它的直方图。看一下直方图的变化。

  • 使用过滤后图像的直方图，决定允许定义沙粒像素，玻璃像素和气泡像素掩蔽的阈限。其他的选项（家庭作业）：写一个函数从直方图的最小值自动判断阈限。

  • 将三种不同的相用不同的颜色上色并显示图片。

  • 用数学形态学清理不同的相。

  • 为所有气泡和沙粒做标签，从沙粒中删除小于10像素的掩蔽。要这样做，用ndimage.sum或np.bincount来计算沙粒大小。

  • 计算气泡的平均大小。

1.5.11.16 图像处理练习：玻璃中的未融化颗粒的答案例子

In [1]:

  1. import numpy as np
  2. import pylab as pl
  3. from scipy import ndimage

1.5.11 科学计算的总结练习 - 图8

  • 打开图像文件MV_HFV_012.jpg并且浏览一下。看一下imshow文档字符串中的参数，用“右”对齐来显示图片（原点在左下角，而不是像标准数组在右上角）。

In [3]:

  1. dat = pl.imread('data/MV_HFV_012.jpg')
  • 修建图片，删除带有测量信息中底部面板。

In [4]:

  1. dat = dat[60:]
  • 用中位数过滤稍稍过滤一下图像以便改进它的直方图。看一下直方图的变化。

In [5]:

  1. filtdat = ndimage.median_filter(dat, size=(7,7))
  2. hi_dat = np.histogram(dat, bins=np.arange(256))
  3. hi_filtdat = np.histogram(filtdat, bins=np.arange(256))

1.5.11 科学计算的总结练习 - 图9

  • 使用过滤后图像的直方图，决定允许定义沙粒像素，玻璃像素和气泡像素掩蔽的阈限。其他的选项（家庭作业）：写一个函数从直方图的最小值自动判断阈限。

In [6]:

  1. void = filtdat <= 50
  2. sand = np.logical_and(filtdat > 50, filtdat <= 114)
  3. glass = filtdat > 114
  • 将三种不同的相用不同的颜色上色并显示图片。

In [7]:

  1. phases = void.astype(np.int) + 2*glass.astype(np.int) + 3*sand.astype(np.int)

1.5.11 科学计算的总结练习 - 图10

  • 用数学形态学清理不同的相。

In [8]:

  1. sand_op = ndimage.binary_opening(sand, iterations=2)
  • 为所有气泡和沙粒做标签，从沙粒中删除小于10像素的掩蔽。要这样做，用ndimage.sum或np.bincount来计算沙粒大小。

In [9]:

  1. sand_labels, sand_nb = ndimage.label(sand_op)
  2. sand_areas = np.array(ndimage.sum(sand_op, sand_labels, np.arange(sand_labels.max()+1)))
  3. mask = sand_areas > 100
  4. remove_small_sand = mask[sand_labels.ravel()].reshape(sand_labels.shape)

1.5.11 科学计算的总结练习 - 图11

  • 计算气泡的平均大小。

In [10]:

  1. bubbles_labels, bubbles_nb = ndimage.label(void)
  2. bubbles_areas = np.bincount(bubbles_labels.ravel())[1:]
  3. mean_bubble_size = bubbles_areas.mean()
  4. median_bubble_size = np.median(bubbles_areas)
  5. mean_bubble_size, median_bubble_size

Out[10]:

  1. (2416.863157894737, 60.0)