矩估计方法用于高斯混合模型的理论分析与遥感实证

Theoretical Analysis and Remote Sensing Empirical Study of the Method of Moments for Gaussian Mixture Models

摘要: 本文系统介绍并实现了基于矩估计(Method of Moments, MoM)用于一维三成分高斯混合模型(Gaussian Mixture Model, GMM)参数估计的方法，推导单分量原点矩的闭式表达，利用1阶到8阶的原点矩和中心距关系，构造了8个关于模型参数的非线性方程，用矩阵最小二乘法得到混合模型的权重的初值，采用Newton-Raphson迭代法求出非线性方程组的解。最后使用中国西部高原纳木错部分区域冰、水、雪、陆地多种物质的光谱进行了验证，针对估计出的每种物质的混合模型概率密度函数，结合距离函数进行了初步分类，分类结果显示了该方法的有效性。

Abstract: This paper systematically introduces and implements the method of moments (MoM) for parameter estimation of a one-dimensional three-component Gaussian mixture model (GMM). We derive the closed-form expressions of raw moments for single Gaussian components, and by leveraging the relationship between the first to eighth raw moments and central moments, we construct eight nonlinear equations with respect to the model parameters. The initial values of the mixture weights are obtained using a matrix least squares method, and the solutions to the nonlinear system are then found via the Newton-Raphson iterative method. Finally, we validate the approach using spectral data of ice, water, snow, and land from the Nam Co region on the Qinghai-Xizang Plateau in western China. For each estimated material’s mixture probability density function, preliminary classification is performed in combination with a distance function. The classification results demonstrate the effectiveness of the proposed method.

1. 引言

对于某一纯净地物类别，在单波段下像素值的直方图分布是单峰的正态分布。因此，如果针对某一类别物质以单一光谱取波段值，则该类别的概率密度函数为单峰的正态分布——只有一个成分[1]。但若针对某一类别物质在至少三个不同光谱上取波段值，则该类别的概率密度函数应为由三个正态分布构成的混合分布，也就是三个正态分布之和，并带有若干比例系数。一个三正态分布混合的概率密度函数可以是单峰、双峰或三峰的。混合中的三个成分分别对应三个光谱波段。混合中每个比例(权重)表示对应成分对混合模型的贡献。在许多遥感地物应用中，单波段的正态分布无法很好地表示一种地物。相比之下，三个光谱波段混合的正态分布一定程度改善了这一缺陷[2]。

目前学术界对高斯混合模型广泛采用的是无监督算法中的EM (Expectation-Maximum)算法用来估计GMM的参数，EM算法通过先猜想缺失数据的分布(E步)，再基于此猜想更新模型参数(M步)，如此反复迭代，逐步逼近最优参数估计。但是因为EM算法严重依赖于初始参数的选择，不同的初始值会导致算法收敛到不同的局部极大值点。在本文中，我们给出一种根据观测数据构造该三成分高斯混合模型的办法，通过多元矩估计的方法建立方程组，最后通过求解方程组来确定混合模型参数的方法。该混合模型在多光谱分类等遥感应用中，可能比现有的多光谱分类方法或其它用途更有用。特别是当缺失数据较多或者“缺失”部分所占的信息量很大时，算法收敛速度可能非常缓慢。

矩估计法是统计学中的一种参数估计技术。Pearson等人[3]引入该法以估计两个一维正态总体混合的均值、方差与比例。基于该技术，Rao等[4]在假设两个一维正态混合分量具有相等方差的前提下给出了估计。Day则将该思路扩展到具有共同协方差矩阵的两个多元正态混合的情形[5]。Cohen用数值方法估计了两个一维正态混合的矩估计[6]。Tan和Chang通过重新参数化得到了矩估计[7]。我们提出将三个正态分布的混合物作为遥感中由三个不同波段(像素)的光谱反射率定义的区域一类的模型。参数估计问题是在所有参数都不受约束的情况下，通过迭代分组数据矩法对三种正态分布的混合物进行参数估计。

2. 高斯混合模型

三种正态分布混合的概率密度函数的形式为：

$f\left(x;{\mu }_1,{\mu }_2,{\mu }_3,{\sigma }_1,{\sigma }_2,{\sigma }_3,p_1,p_2,p_3\right)={\displaystyle {\sum }_{i=1}^3{p}_i{f}_i\left(x;{\mu }_i,{\sigma }_i\right)}$ (1)

$x\in ℝ$ , ${\mu }_i\in ℝ$ , ${\sigma }_i\in \left(0,+\infty \right)$ , $p_i\in \left(0,1\right)$

其中的$p_i$ 代表混合成分的比例，满足${\displaystyle {\sum }_{i=1}^3{p}_i=1}$ ；混合模型的成分函数$f_i\left(x;{\mu }_i,{\sigma }_i\right)$ 的函数形式为：

$f_i\left(x;{\mu }_i,{\sigma }_i\right)=\frac{1}{\sqrt{2\pi {\sigma }_i}}{\text{e}}^{-\frac{{\left(x-{\mu }_i\right)}^2}{2{\sigma }_i^2}}$ (2)

在上式中，x是随机变量，代表某类物质像素维度的三个波谱的数值反射率结果，$f\left(x\right)$ 是特定物质的不同波谱的概率密度函数，$\mu$ 代表不同波谱的均值，$\sigma$ 代表不同波谱的标准差。

3. 参数估计方法

参数估计的过程分为三个步骤。第一步，构建方程组以估计中的参数。第二步，计算合适的初始值，用于求解第一步中的方程组。第三步，采用牛顿迭代法，利用第二步求得的初始值，找到第一步构建的方程组的解。

3.1. 线性方程组的构建

矩母函数也被成矩生成函数，当X是离散变量时，他的形式为：

$M_x\left(t\right)={\displaystyle \sum {\text{e}}^{tx}f\left(x\right)}$ (3)

当x是连续变量时，他的形式为：

$M_x\left(t\right)={\displaystyle \int {\text{e}}^{tx}f\left(x\right)\text{d}x}$ (4)

当$f\left(x\right)$ 为上式(2)的正态分布的形式时，可推导出

$\begin{array}{c}{M}_x\left(t\right)={\displaystyle \int \frac{1}{\sqrt{2\pi }\sigma }{\text{e}}^{-\frac{{\left(x-\mu \right)}^2}{2{\sigma }^2}}}\cdot {\text{e}}^{tx}\text{d}x\\ =\frac{1}{\sqrt{2\pi }\sigma }{\displaystyle \int {\text{e}}^{-\frac{{\left(x-\mu \right)}^2}{2{\sigma }^2}+tx}}\text{d}x\end{array}$ (5)

令$\omega =\frac{x-\mu }{\sigma }$ ，于是$\text{d}x=\sigma \text{d}\omega$

$\begin{array}{c}{M}_x\left(t\right)=\frac{1}{\sqrt{2\pi }}{\displaystyle \int {\text{e}}^{-\frac{{\omega }^2}{2}+t\sigma \omega +t\mu }}\text{d}\omega \\ =\frac{1}{\sqrt{2\pi }}{\text{e}}^{t\mu }{\displaystyle \int {\text{e}}^{-\frac{{\omega }^2}{2}+t\sigma \omega +t\mu }}\text{d}\omega \\ =\frac{1}{\sqrt{2\pi }}{\text{e}}^{t\mu }{\displaystyle \int {\text{e}}^{-\frac{{\omega }^2-2t\sigma \omega +t^2{\sigma }^2}{2}}\cdot {\text{e}}^{\frac{t^2{\sigma }^2}{2}}}\text{d}\omega \\ =\frac{1}{\sqrt{2\pi }}{\text{e}}^{t\mu +\frac{t^2{\sigma }^2}{2}}{\displaystyle \int {\text{e}}^{-\frac{{\left(\omega -t\sigma \right)}^2}{2}}}\text{d}\omega \\ ={\text{e}}^{t\mu +\frac{t^2{\sigma }^2}{2}}\end{array}$ (6)

当$t=0$ 时，对应的是随机变量x的原点距，于是样本r阶原点矩为$M_r={M_X^{\left(r\right)}\left(t\right)|}_{t=0}$ ，便得到了8个对应的原点距公式$M_r=E\left(X^r\right)$ ，($r=1,2,\cdots ,8$ )。

设样本总体均值$\overline{\mu }=M_1$ ，则样本r阶中心距$A_r=E\left({\left(X-\overline{\mu }\right)}^r\right)$ ，再把模型的原点矩估计和样本的中心矩估计对应起来[8]，得到8个方程$g_r\left(\Theta \right)=M_r\left(\Theta \right)-A_r=0$ ，$r=1,\cdots ,8$ 。详细的形式为

$\begin{array}{l}{g}_1={\displaystyle \sum _{i=1}^3{p}_i{\mu }_i-A_1=0;}\\ g_2={\displaystyle \sum _{i=1}^3{p}_i\left({\mu }_i^2+{\sigma }_i^2\right)-A_2=0;}\\ g_3={\displaystyle \sum _{i=1}^3{p}_i\left({\mu }_i^3+3{\sigma }_i^2{\mu }_i\right)-A_3=0;}\\ g_4={\displaystyle \sum _{i=1}^3{p}_i\left({\mu }_i^4+6{\sigma }_i^2{\mu }_i^2+3{\sigma }_i^4\right)-A_4=0;}\\ g_5={\displaystyle \sum _{i=1}^3{p}_i\left({\mu }_i^5+10{\sigma }_i^2{\mu }_i^3+15{\sigma }_i^4{\mu }_i\right)-A_5=0;}\\ g_6={\displaystyle \sum _{i=1}^3{p}_i\left({\mu }_i^6+15{\sigma }_i^2{\mu }_i^4+45{\sigma }_i^4{\mu }_i^2+15{\sigma }_i^6\right)-A_6=0;}\\ g_7={\displaystyle \sum _{i=1}^3{p}_i\left({\mu }_i^7+21{\sigma }_i^2{\mu }_i^5+105{\sigma }_i^4{\mu }_i^3+105{\sigma }_i^6{\mu }_i\right)-A_7=0;}\\ g_8={\displaystyle \sum _{i=1}^3{p}_i\left({\mu }_i^8+28{\sigma }_i^2{\mu }_i^6+210{\sigma }_i^4{\mu }_i^4+420{\sigma }_i^6{\mu }_i^2+105{\sigma }_i^8\right)-A_8=0;}\end{array}$ (7)

其中

$\begin{array}{l}{A}_1=m_1;\\ A_2=m_1^2+m_2;\\ A_3=m_1^3+3m_1{m}_2+m_3;\\ A_4=m_1^4+6m_1^2{m}_2+4m_1{m}_3+m_4;\\ A_5=m_1^5+10m_1^2{m}_3+10m_1^3{m}_2+5m_1{m}_4+m_5;\\ A_6=m_1^6+15m_1^4{m}_2+20m_1^3{m}_3+15m_1^2{m}_4+6m_1{m}_5+m_6;\\ A_7=m_1^7+21m_1^5{m}_2+35m_1^4{m}_4+35m_1^3{m}_4+21m_1^2{m}_5+7m_1{m}_6+m_7;\\ A_6=m_1^8+28m_1^6{m}_2+56m_1^5{m}_3+70m_1^4{m}_4+56m_1^3{m}_5+28m_1^2{m}_6+8m_1{m}_7+m_8;\end{array}$ (8)

在上式(8)中，$m_r=\frac{1}{n}{\displaystyle \sum _{i=1}^n{\left(x_i-\overline{x}\right)}^r}$ ， $r=2,3,\cdots ,8$ . $m_1=\overline{x}=\frac{1}{n}{\displaystyle \sum _{i=1}^n{x}_i}$

在(7)式中构建非线性方程组中，其中$p_3=1-p_1-p_2$ ，方程组中的未知数个数包括${\mu }_1$ ，${\mu }_2$ ，${\mu }_3$ ，${\sigma }_1$ ，${\sigma }_2$ ，${\sigma }_3$ ，$p_1$ ，$p_2$ 共8个，结合8个方程式，在合理设置初始值的情况下，利用牛顿迭代法便可以求出唯一解，下面分析如何确定初始值。

3.2. 确定牛顿迭代法的合理初始值

在选择迭代起始值时，混合模型的均值和标准差参数应分别使用样本波谱维度的均值和标准差。对于混合比例，则使用通过最小二乘法得出的样本分组的均值和标准差。

3.2.1. 为${\mu }_1$ ，${\mu }_2$ ，${\mu }_3$ 和${\sigma }_1$ ，${\sigma }_2$ ，${\sigma }_3$ 寻找合适的初始值

样本的均值和标准差可作为(7)式中参数${\mu }_1$ ，${\mu }_2$ ，${\mu }_3$ 和${\sigma }_1$ ，${\sigma }_2$ ，${\sigma }_3$ 的合适初始值(估计值)，根据样本的三个不同波段的光谱数据$x_i,\left(i=1,2,3\right)$ ，计算每个光谱的均值和方差结果${\mu }_{i0},{\sigma }_{i0}\left(i=1,2,3\right)$ ，令${\mu }_{10}$ ，${\mu }_{20}$ ，${\mu }_{30}$ 和${\sigma }_{10}$ ，${\sigma }_{20}$ ，${\sigma }_{30}$ 分别表示${\mu }_1$ ，${\mu }_2$ ，${\mu }_3$ 和${\sigma }_1$ ，${\sigma }_2$ ，${\sigma }_3$ 的初始值[9]。

3.2.2. 为$p_1$ ，$p_2$ 寻找合适初始值

令$p_1$ 和$p_2$ 合适的初始值分别为$p_{10}$ 和$p_{20}$ ，但我们并不知道$p_{10}$ 和$p_{20}$ 的值。我们希望找到这些值。为此，首先将值${\mu }_{10}$ ，${\mu }_{20}$ ，${\mu }_{30}$ 和${\sigma }_{10}$ ，${\sigma }_{20}$ ，${\sigma }_{30}$ 分别代入${\mu }_1$ ，${\mu }_2$ ，${\mu }_3$ 和${\sigma }_1$ ，${\sigma }_2$ ，${\sigma }_3$ 。接着，也将$p_{10}$ ， $p_{20}$ 和$p_{30}$ $\left(p_{30}=1-p_{10}-p_{20}\right)$ 分别代入(7)式中的$p_1$ ，$p_2$ 和$p_3$ 。如此，我们得到：

(9)

将3.2.1中已知的${\mu }_{10}$ ，${\mu }_{20}$ ，${\mu }_{30}$ 和${\sigma }_{10}$ ，${\sigma }_{20}$ ，${\sigma }_{30}$ 带入(9)式后得到的方程组转化为矩阵形式，

$\left(\begin{array}{cc}{a}_{11}& a_{12}\\ ⋮& ⋮\\ a_{81}& a_{82}\end{array}\right)\left(\begin{array}{c}{p}_{10}\\ p_{20}\end{array}\right)=\left(\begin{array}{c}{A}_1^{\ast }\\ ⋮\\ A_8^{\ast }\end{array}\right)$ (10)

结合最小二乘法，可以得到$p_{10}$ ，$p_{20}$ 的值为$P_0={\left(X^{\text{T}}X\right)}^{-1}{X}^{\text{T}}A$ 。

3.3. 牛顿迭代法求解非线性方程组

牛顿–拉弗森迭代法被用于求解(9)中构建的方程组，其中使用的合适初始值是通过3.2节中的方法找到的。

我们使用python程序用于参数估计和混合模型构建过程。该程序通过牛顿-拉弗森迭代过程求解包含八个未知数的八个联立非线性方程。详细的算法流程如下，

Step1：初值初始化。设当前解向量$x_k=x_0$ 。

Step2：迭代循环。

while (k < MAXITER) do:

计算函数值和雅可比矩阵

计算当前点$x_k$ 处的函数值向量$F_k=F\left(x_k\right)$ 和雅可比矩阵$J_k=J\left(x_k\right)$ 。

求解线性方程组，构建线性系统以求解释向量$\Delta x_k$ ，

$J_k\Delta x_k=-F_k$

更新未知数向量：$x_{k+}{}_1=x_k+\Delta x_k$

Step3：检查收敛性

如果$\text{norm}\left(\Delta x_k\right)<\epsilon$ ，则退出循环。这表示解的变化已经非常微小，趋于稳定。

Step4：递增计数器

$k=k+1$

Step5：输出结果

4. 方法对比

为了验证本文所提出的方法相比于EM方法在估计混合模型参数上的优势，在相同的高斯混合分布离散数据下，我们分别使用本文提出的矩估计方法和经典EM方法进行了参数对比，并量化对比了不同混合分布的参数偏离程度，

具体的量化函数为平均绝对误差(MAE)和均方根误差(RMSE)，详细的计算公式为

$\text{MAE}\left(X,h\right)=\frac{1}{m}{\displaystyle \sum _{i=1}^m\left|h\left(x_i\right)-y_i\right|}$ . (11)

$\text{RMSE}\left(X,h\right)=\sqrt{\frac{1}{m}{\displaystyle \sum _{i=1}^m{\left(h\left(x_i\right)-y_i\right)}^2}}$ . (12)

其中$h\left(x_i\right)$ 和$y_i$ 分别代表预测值和真实值，我们用上述两个评价指标来衡量预测的混合模型参数与真实之间的偏移差距。

具体的对比结果如下：

参数	EM估计结果	矩估计结果	真实参数值
$p_1$	0.33082458	0.3416441	0.333333333
$p_2$	0.33590012	0.32502256	0.333333333
$p_3$	0.3332753	0.33333333	0.333333333
$m_1$	32.16006551	32.21282772	32
$m_2$	39.81166837	39.94827327	40
$m_3$	62.10483519	62.07487003	62
$s_1$	1.65688031	1.68577679	1.6
$s_2$	2.09770345	1.89546207	2
$s_3$	0.78789696	0.80384005	0.8
MAE (EM)	0.625052703
MAE (Ours)	0.550200793
RMSE (EM)	0.08502089
RMSE (Ours)	0.072015535

通过上述结果，可以发现，在MAE和RMSE的指标度量下，本方法在预测值和真实结果的差距都更为接近，证明基于本文矩估计方法的优势。

5. 实验验证

为了验证使用该方法估计出的物质的高斯混合模型的结果，我们选取了中国西部青藏高原纳木错区域的卫星影像，选取的冰、水、雪、陆地四种物质的波谱，分别使用上述方法对每种物质的多光谱影像的第2、3、4索引位置的波段(Blue、Green、Red)进行高斯混合模型的参数估计，最后使用估计的参数进行影像分类[10]-[12]，分类的规则是：对于待分类物质的${\mu }_1$ ，${\mu }_2$ ，${\mu }_3$ 和${\sigma }_1$ ，${\sigma }_2$ ，${\sigma }_3$ 参数，与参考物质的混合模型参数${\mu }_1{}^i$ ，${\mu }_2{}^i$ ，${\mu }_3{}^i$ 和${\sigma }_1{}^i$ ，${\sigma }_2{}^i$ ，${\sigma }_3{}^i$ 进行距离计算，距离计算的方法使用的是曼哈顿距离计算，根据待分类像素的波谱混合模型与参考物质的波谱混合模型的距离相似度进行分类[13] [14]，把每个像素位置分类到与参考波谱的混合模型参数最接近的类别中。距离计算方法为：

$D_j={\displaystyle \sum _{i=1}^3\left|{\mu }_i^j-{\mu }_i\right|}+\left|{\sigma }_i^j-{\sigma }_i\right|$ . ($j=1,2,3,4$ ) (13)

经过实验计算，四种物质的波谱混合模型参数分别为：

冰：$\mu =\left[1763.12,2043.05,2172.60\right]$ ；$\sigma =\left[322.17,302.81,321.88\right]$ ；$p=\left[0.72,0.03,0.25\right]$

雪：$\mu =\left[6893.41,7099.73,7266.87\right]$ ；$\sigma =\left[404.18,424.25,415.19\right]$ ；$p=\left[0.51,0.01,0.49\right]$

水：$\mu =\left[361.96,678.21,741.25\right]$ ；$\sigma =\left[7.18,13.86,15.15\right]$ ；$p=\left[0.33,0.34,0.33\right]$

陆：$\mu =\left[1020.85,1110.67,1742.20\right]$ ；$\sigma =\left[60.55,46.38,95.13\right]$ ；$p=\left[0.28,0.39,0.33\right]$

对应的四种物质的高斯混合模型分布图如下图1所示。

Figure 1. Mixture model distributions of four reference land-cover types

图1. 四种参考地物的混合模型分布图

基于(11)式的距离相似度分类结果如下图2所示。

Figure 2. Multispectral remote sensing land-cover separation results (Nam Co)

图2. 多光谱遥感图像地物分离结果(纳木错)

6. 结论

本文重现并实现了基于矩估计(method of moments)的一维三成分高斯混合模型参数估计方法，并将其应用于纳木错区域的遥感多光谱数据中以对冰、水、雪、陆地进行建模与初步分类。方法要点包括：利用单个高斯成分的原点矩闭式表达构建混合分布的1至8阶原点矩，再通过原点矩与总体均值间的转换得到关于样本中心矩的8个非线性方程；采用矩阵最小二乘法获得混合权重的初值，并以牛顿–拉弗森迭代求解方程组。实证结果表明，该方法能较好拟合各类物质的边缘亮度分布并提供具有解释性的模型参数(均值、方差与成分权重)，在参数估计充足且成分分离较好的条件下，其拟合结果与EM方法具有可比性。基于估计的混合模型进行的像素分类(采用模型对像素对数似然比较或模型间距离度量)在实验中实现了对冰、水、雪与陆地的区分，验证了方法在遥感分类任务中的实用性。