CN103049653A - 基于em算法的g0分布参数最大似然估计方法 - Google Patents

Abstract

本发明公开了一种基于EM算法的G0分布参数最大似然估计方法，首先使用矩估计法对G0分布参数进行估计，然后是将矩估计法得到的参数估计值作为初始值，采用EM算法，以迭代求解的方式估计G0分布参数。本发明所设计的基于EM算法的G0分布参数最大似然估计方法具有较高的参数估计精度。

Description

基于EM算法的G0分布参数最大似然估计方法

技术领域

本发明属于合成孔径雷达图像解译领域，涉及一种基于EM算法的G0分布参数最大似然估计方法。

背景技术

由于合成孔径雷达（synthetic aperture radar，SAR）图像数据具有噪声严重、背景杂波复杂的特点，因此基于统计模型来展开SAR图像解译工作得到了广泛的关注。统计模型对实测SAR图像数据统计特性描述的准确性会在很大程度上影响SAR图像解译的性能，为此各国学者发展了很多用于描述SAR图像数据的统计模型。其中，G0分布，具有适用范围广、建模能力强优点，近年来被广泛地应用于SAR图像的解译中。

G0分布在SAR图像解译中的应用离不开对其进行参数估计。所谓参数估计，就是利用观测得到的SAR图像数据来估计G0分布的参数。设观测得到的SAR图像幅度数据为：y={y_t,t=1,…,T}，

且独立同分布，则y_t的G0分布表达式为：

$p\left(y_t\right|\mathrm{\λ})=\frac{2^{1+\mathrm{\β}}{y}_t^{N-1}{\mathrm{\σ}}^{\mathrm{\β}}}{\mathrm{\Γ}(N/2)\mathrm{\Γ}\left(\mathrm{\β}\right)}{(y_t^2+2\mathrm{\σ})}^{-(N/2+\mathrm{\β})}\mathrm{\Γ}(N/2+\mathrm{\β}),\mathrm{\β},\mathrm{\σ},N,y_t>0---\left(1\right)$

其中，参数N为SAR图像数据视数的2倍，可以先验获得不需要估计。因此，参数估计即为利用观测得到的SAR图像幅度数据y={y_t,t=1,…,T}来估计G0分布中的参数β和σ。

目前，G0分布最常用的参数估计方法是矩估计法（method of moments,MoM）。矩估计法，计算相对比较简单，但是精度不高，并且由于其方法自身固有的限制，不能实现G0分布全范围的参数估计。G0分布另外一种参数估计方法是基于Mellin变化的，该方法可以实现G0分布全范围的参数估计，但是参数估计的精度同样不是很高。

最大似然（Maximum Likelihood，ML）估计是统计意义上最优的参数估计方法，但是由于G0分布表达式的复杂性，最大似然估计方法给出的方程组具有高度的非线性，肯定没有解析解。而如果运用数值方法来求解最大似然估计方程组，则面临着计算量巨大和可能无法收敛到正确解的困境。这些都限制了最大似然估计方法在G0分布上的运用。但是，作为统计意义上最优的方法，最大似然估计法可能会消耗更多的计算量，理论上却是一定会给出更准确的参数估计。本专利给出了一种基于EM算法的G0分布参数最大似然估计方法。EM算法是一种迭代的方法来寻找统计模型的最大似然估计，常常用于最大似然估计的方程无法直接求解的情况。

发明内容

本发明所要解决的技术问题是提供一种具有较高的参数估计精度的基于EM算法的G0分布最大似然参数估计方法。

本发明为解决上述技术问题采用如下技术方案：本发明设计了一种基于EM算法的G0分布参数最大似然估计方法，包括如下具体步骤：

步骤（1）：采用矩估计法，为观测得到的SAR图像生成参数估计初始值；

步骤（2）：根据步骤（1）中的参数估计初始值，采用迭代算法估计G0分布参数。

作为本发明的一种优化方法：所述步骤（1）还包括如下具体步骤：

步骤（11）：设观测得到的SAR图像幅度数据为：y={y_t,t=1,…,T}，

且独立同分布，则y_t的G0分布表达式为：

$p\left(y_t\right|\mathrm{\λ})=\frac{2^{1+\mathrm{\β}}{y}_t^{N-1}{\mathrm{\σ}}^{\mathrm{\β}}}{\mathrm{\Γ}(N/2)\mathrm{\Γ}\left(\mathrm{\β}\right)}{(y_t^2+2\mathrm{\σ})}^{-(N/2+\mathrm{\β})}\mathrm{\Γ}(N/2+\mathrm{\β}),\mathrm{\β},\mathrm{\σ},N,y_t>0$

其中；N定义为等效视数，σ定义为形状参数，β定义为尺度参数,t定义为SAR图像幅度数据的序号，T定义为SAR图像幅度数据的个数；

步骤（12）：利用矩估计法来估计参数初始值，采用如下公式：

$\left\{\begin{array}{c}\widehat{\mathrm{\β}}=1+\frac{N{\hat{m}}_4}{N{\hat{m}}_4-(N+2){\hat{m}}_2^2}\\ \widehat{\mathrm{\σ}}=\frac{(\widehat{\mathrm{\β}}-1){\hat{m}}_2}{N}\end{array}\right.$

其中，

定义为k阶样本矩，即有

得到尺度参数估计初始值和形状参数估计初始值

作为本发明的一种优化方法：所述步骤（2）包括如下具体处理：

步骤（21）：令

建立如下函数变量A和G：

$A=\frac{1}{T}\underset{t=1}{\overset{T}{\mathrm{\Σ}}}\frac{N+2{\mathrm{\β}}^{\′}}{y_t^2+2{\mathrm{\σ}}^{\′}},$

$G={\left(\underset{t=1}{\overset{T}{\mathrm{\Π}}}\right[\frac{y_t^2+2{\mathrm{\σ}}^{\′}}{2}\left]\right)}^{1/T}\exp (-\mathrm{\Ψ}\left(\frac{N+2{\mathrm{\β}}^{\′}}{2}\right)),$

其中，Ψ(·)定义为digamma函数；

步骤（22）：令利用步骤（21）中的函数变量A和G与公式

进行迭代，当满足迭代终止条件后，结束迭代过程，令参数估计值

等于最后一次迭代所得的

则σ的估计值由公式

求得，并将迭代终止时的

和

作为最终的参数估计值；若不满足迭代终止条件，则返回步骤（21）。

本发明与现有技术先比具有如下优点：

本发明所设计的基于EM算法的G0分布最大似然参数估计方法的优点为：通过EM算法可以实现对G0分布的参数最大似然估计，该方法比现有的矩估计法和基于Mellin变化的方法具有更高的参数估计精度。

附图说明

图1为本发明所设计的基于EM算法的G0分布最大似然参数估计方法的流程示意图。

具体实施方式

下面结合附图对本发明做进一步的详细说明：

本发明提出一种基于EM算法的G0分布参数最大似然估计方法。由于SAR图像数据具有噪声严重、背景杂波复杂的特点，因此基于统计模型来展开SAR图像解译工作得到了广泛的关注。Frery等人给出了一种新的统计分布模型G分布，G分布的特殊形式G0分布，具有适用范围广、参数估计容易的优点。其中，参数估计是G0分布研究的一个核心问题。矩估计法以及基于Mellin变换的参数估计方法是目前常用的参数估计方法。但是最大似然估计作为统计意义上最优的参数估计方法由于表达式的复杂性一直没有得到应用。为了解决这个问题，本专利提出了基于EM算法的G0分布参数估计方法，该方法是通过迭代的方式来寻找G0分布参数的最大似然估计。

如图1所示，本发明设计了一种基于EM算法的G0分布参数最大似然估计方法，包括以下步骤：

步骤1：采用矩估计法，生成参数估计的初值，具体方法如下：

$\left\{\begin{array}{c}\widehat{\mathrm{\β}}=1+\frac{N{\hat{m}}_4}{N{\hat{m}}_4-(N+2){\hat{m}}_2^2}\\ \widehat{\mathrm{\σ}}=\frac{(\widehat{\mathrm{\β}}-1){\hat{m}}_2}{N}\end{array}\right.,---\left(2\right)$

其中

表示k阶样本矩，即有

${\hat{m}}_k=\frac{1}{T}\underset{t=1}{\overset{T}{\mathrm{\Σ}}}{y}_t^k---\left(3\right)$

步骤2：根据步骤1的矩估计法给出的初始值，采用EM算法以迭代求解的方式估计G0分布参数，具体方法如下:

a)令

根据下面的（4）式和（5）式求取函数变量A和G，

$A=\frac{1}{T}\underset{t=1}{\overset{T}{\mathrm{\Σ}}}\frac{N+2{\mathrm{\β}}^{\′}}{y_t^2+2{\mathrm{\σ}}^{\′}},---\left(4\right)$

$G={\left(\underset{t=1}{\overset{T}{\mathrm{\Π}}}\right[\frac{y_t^2+2{\mathrm{\σ}}^{\′}}{2}\left]\right)}^{1/T}\exp (-\mathrm{\Ψ}\left(\frac{N+2{\mathrm{\β}}^{\′}}{2}\right)),---\left(5\right)$

其中，Ψ(·)表示digamma函数；

b)令

利用（6）式进行迭代，当满足迭代终止条件后，结束迭代过程，令β的估计值

等于最后一次迭代所得的

${\widehat{\mathrm{\β}}}_{(k+1)}=\frac{[\ln {\widehat{\mathrm{\β}}}_{\left(k\right)}-\mathrm{\Ψ}\left({\widehat{\mathrm{\β}}}_{\left(k\right)}\right)]}{\ln \left(\mathrm{AG}\right)}{\widehat{\mathrm{\β}}}_{\left(k\right)}---\left(6\right)$

c)由（7）式求取σ的估计值

$\widehat{\mathrm{\σ}}=\frac{\widehat{\mathrm{\β}}}{A}---\left(7\right)$

d)满足迭代终止条件则结束迭代过程，将迭代终止时的

和

作为最终的参数估计值；若不满足迭代终止条件，则回到a)。

为了进一步说明本发明所提出的方法，下面给出本发明方法的理论推导过程。

G0分布模型是由Frery等人根据经典的乘积模型发展而得。就SAR图像幅度数据而言，可以表示为符合单位均值方根Gamma分布Γ^1/2(n,n)的斑点噪声分量与符合逆方根Gamma分布Γ^-1/2(α,γ)的表征地物RCS起伏特性的后向散射幅度的乘积。设观测得到SAR图像幅度数据为y={y_t,t=1,…,T}，

且独立同分布，则y_t的G0分布模型表达式为：

$p\left(y_t\right|\mathrm{\α},\mathrm{\γ})=\frac{{2n}^n\mathrm{\Γ}(n-\mathrm{\α})y_t^{2n-1}}{\mathrm{\Γ}\left(n\right)\mathrm{\Γ}(-\mathrm{\α}){\mathrm{\γ}}^{\mathrm{\α}}{(\mathrm{\γ}+{\mathrm{ny}}_t^2)}^{n-\mathrm{\α}}},-\mathrm{\α},\mathrm{\γ},n,y_t>0---\left(8\right)$

其中，n表示等效视数，α为形状参数反映了被测区域的均匀度，γ为尺度参数与被测区域的平均能量有关，t定义为SAR图像幅度数据的序号，T定义为SAR图像幅度数据的个数，式（8）是当前关于G0分布的一种常用的表达式。

本专利的方法是基于一种新的G0分布表达式及其推导过程。由贝叶斯理论，我们可以将G0分布写成

$p\left(y_t\right|\mathrm{\λ})={\∫}_0^{\∞}p\left(y_t\right|{\mathrm{\ω}}_t)p\left({\mathrm{\ω}}_t\right|\mathrm{\λ})d{\mathrm{\ω}}_t,---\left(9\right)$

假设变量y_t为N个独立同分布的高斯随机变量（均值为零，方差为ω_t）的平方和的平方根，则有p(y_t|ω_t)符合generalized Rayleigh分布，即

$p\left(y_t\right|{\mathrm{\ω}}_t)=\frac{2y_t^{N-1}}{{\left(2{\mathrm{\ω}}_t\right)}^{N/2}\mathrm{\Γ}(N/2)}\exp \left(\frac{-y_t^2}{2{\mathrm{\ω}}_t}\right)---\left(10\right)$

假设参数ω_t符合inverse Gamma分布，即有

$p\left({\mathrm{\ω}}_t\right|\mathrm{\λ})=\frac{{\mathrm{\σ}}^{\mathrm{\β}}}{\mathrm{\Γ}\left(\mathrm{\β}\right)}{\mathrm{\ω}}_t^{-\mathrm{\β}-1}\exp \left(\frac{-\mathrm{\σ}}{{\mathrm{\ω}}_t}\right)---\left(11\right)$

式中λ=(β，σ)为inverse Gamma分布的参数。将（10）和（11）带入（9）式，积分后可以得到

$p\left(y_t\right|\mathrm{\λ})=\frac{2^{1+\mathrm{\β}}{y}_t^{N-1}{\mathrm{\σ}}^{\mathrm{\β}}}{\mathrm{\Γ}(N/2)\mathrm{\Γ}\left(\mathrm{\β}\right)}{(y_t^2+2\mathrm{\σ})}^{-(N/2+\mathrm{\β})}\mathrm{\Γ}(N/2+\mathrm{\β})---\left(12\right)$

将（12）式与传统的G0分布表达式（8）式相比较，得

$\left\{\begin{array}{c}N=2n\\ \mathrm{\β}=-\mathrm{\α}\\ \mathrm{\σ}=\mathrm{\γ}/N\end{array}\right.,---\left(13\right)$

由此可见（12）式也是一种G0分布的表达式，本专利中的参数估计是基于（12）式的。上述的G0分布的推导也相当于提供了一种新的方式来解释符合G0分布的SAR数据。

由最大似然估计法可知，G0分布的最大似然参数估计使得

$\widehat{\mathrm{\λ}}=\underset{\mathrm{\λ}}{\arg \max }p\left(y\right|\mathrm{\λ})$

(14)

$=\underset{\mathrm{\λ}}{\arg \max }\underset{t=1}{\overset{T}{\mathrm{\Π}}}p\left(y_t\right|\mathrm{\λ})$

由于没有（14）式的解析解，我们采用EM算法来求解（14）式的最大似然估计值。EM算法是一种迭代的算法，根据当前的参数估计值产生一个新的参数估计值，并且新的参数估计值相比当前参数估计值具有更大的似然。EM算法的特性保证了，在达到稳定点之前，新的参数估计值总是比当前的参数估计值具有更大的似然。

EM算法由最大化下面的辅助函数得到：

$\widehat{\mathrm{\λ}}=\underset{\mathrm{\λ}}{\arg \max }{\∫}_{\mathrm{\ω}}p\left(\mathrm{\ω}\right|y,{\mathrm{\λ}}^{\′})\log p(y,\mathrm{\ω}|\mathrm{\λ})\mathrm{d\ω}---\left(15\right)$

λ′为当前的参数估计值。由(15)式可以得到：

$\widehat{\mathrm{\λ}}=\underset{\mathrm{\λ}}{\arg \max }\underset{t=1}{\overset{T}{\mathrm{\Σ}}}{\∫}_0^{\∞}p\left(w_t\right|y_t,{\mathrm{\λ}}^{\′})\log p\left({\mathrm{\ω}}_t\right|\mathrm{\λ})d{\mathrm{\ω}}_t---\left(16\right)$

将（11）式代入（16）式,将其等号右部分别对σ和β进行求导并令其为0，可得

$\frac{1}{\widehat{\mathrm{\σ}}}=\frac{1}{\widehat{\mathrm{\β}}T}\underset{t=1}{\overset{T}{\mathrm{\Σ}}}{\∫}_0^{\∞}p\left({\mathrm{\ω}}_t\right|y_t,{\mathrm{\λ}}^{\′})\frac{1}{{\mathrm{\ω}}_t}d{\mathrm{\ω}}_t---\left(17\right)$

$\mathrm{\Ψ}\left(\widehat{\mathrm{\β}}\right)=\ln \left(\widehat{\mathrm{\σ}}\right)-\frac{1}{T}\underset{t=1}{\overset{T}{\mathrm{\Σ}}}{\∫}_0^{\∞}p\left({\mathrm{\ω}}_t\right|y_t,{\mathrm{\λ}}^{\′})\ln {\mathrm{\ω}}_{t}d{\mathrm{\ω}}_t---\left(18\right)$

为了求解（17），（18）两式中的积分，首先求解下式：

${\∫}_0^{\∞}p\left({\mathrm{\ω}}_t\right|y_t,{\mathrm{\λ}}^{\′}){\mathrm{\ω}}_t^{s}d{\mathrm{\ω}}_t=\frac{1}{p\left(y_t\right|{\mathrm{\λ}}^{\′})}{\∫}_0^{\∞}p\left(y_t\right|{\mathrm{\ω}}_t)p\left({\mathrm{\ω}}_t\right|{\mathrm{\λ}}^{\′}){\mathrm{\ω}}_t^{s}d{\mathrm{\ω}}_t---\left(19\right)$

将（10）、（11）、（12）式带入可得

${\∫}_0^{\∞}p\left({\mathrm{\ω}}_t\right|y_t,{\mathrm{\λ}}^{\′}){\mathrm{\ω}}_t^{s}d{\mathrm{\ω}}_t=2^{-s}{(y_t^2+2{\mathrm{\σ}}^{\′})}^s\frac{\mathrm{\Γ}(\frac{N}{2}+{\mathrm{\β}}^{\′}-s)}{\mathrm{\Γ}(\frac{N}{2}+{\mathrm{\β}}^{\′})}---\left(20\right)$

令（20）式中s=-1，则可得（17）式所需的

${\∫}_0^{\∞}p\left({\mathrm{\ω}}_t\right|y_t,{\mathrm{\λ}}^{\′}){\mathrm{\ω}}_t^{-1}d{\mathrm{\ω}}_t=2{(y_t^2+2{\mathrm{\σ}}^{\′})}^{-1}(\frac{N}{2}+{\mathrm{\β}}^{\′})---\left(21\right)$

利用关系式 $\∫\ln \mathrm{\ωf}\left(\mathrm{\ω}\right)\mathrm{d\ω}=(\∂/\∂s){\∫\mathrm{\ω}}^{s}f\left(\mathrm{\ω}\right)\mathrm{d\ω}{|}_{s=0},$ 可以得到（18）式所需的

${\∫}_0^{\∞}p\left({\mathrm{\ω}}_t\right|y_t,{\mathrm{\λ}}^{\′})\ln {\mathrm{\ω}}_{t}d{\mathrm{\ω}}_t=\left(\frac{\∂}{\∂s}\right){\∫}_0^{\∞}p\left({\mathrm{\ω}}_t\right|y_t,{\mathrm{\λ}}^{\′}){\mathrm{\ω}}_t^{s}d{\mathrm{\ω}}_t{|}_{s=0}$ (22)

$=\ln \left(\frac{y_t^2+2{\mathrm{\σ}}^{\′}}{2}\right)-\mathrm{\Ψ}(\frac{N}{2}+{\mathrm{\β}}^{\′})$

将（21），（22）的结果代入式（17）和（18）可得

$\ln \widehat{\mathrm{\β}}-\mathrm{\Ψ}\left(\widehat{\mathrm{\β}}\right)=\ln \left(\mathrm{AG}\right)---\left(23\right)$

$\widehat{\mathrm{\σ}}=\frac{\widehat{\mathrm{\β}}}{A}---\left(24\right)$

其中，

$A=\frac{1}{T}\underset{t=1}{\overset{T}{\mathrm{\Σ}}}\frac{N+2{\mathrm{\β}}^{\′}}{y_t^2+2{\mathrm{\σ}}^{\′}}---\left(25\right)$

$G={\left(\underset{t=1}{\overset{T}{\mathrm{\Π}}}\right[\frac{y_t^2+2{\mathrm{\σ}}^{\′}}{2}\left]\right)}^{1/T}\exp (-\mathrm{\Ψ}\left(\frac{N+2{\mathrm{\β}}^{\′}}{2}\right))---\left(26\right)$

式（23）～（26）构成了EM算法的迭代求解过程来求解G0分布参数，在给出初始的σ′和β′值时（本专利中，初始的σ′和β′值为矩估计法得到参数估计值），迭代求解过程就可以开始。迭代求解过程可以在满足某些设定的条件的情况下停止。方程（23）没有解析解，本专利采用了不动点迭代（fixed-point iteration）的方法来求解方程（23），这种方法较常用的Newton-Raphson方法简单，而且可以避免求解digamma函数的导数。对于一个给定的ln(AG)值，不动点迭代法给出的迭代求解方程为：

${\widehat{\mathrm{\β}}}_{(k+1)}=\frac{[\ln {\widehat{\mathrm{\β}}}_{\left(k\right)}-\mathrm{\Ψ}\left({\widehat{\mathrm{\β}}}_{\left(k\right)}\right)]}{\ln \left(\mathrm{AG}\right)}{\widehat{\mathrm{\β}}}_{\left(k\right)}.---\left(27\right)$

该迭代过程的起始值可以为当前的估计值β′。由（27）的迭代过程求得

后，

可以直接由（24）式求得，然后将求得

和

值代替β′和σ′，完成一次EM迭代。不难看出，整个参数求解过程有两个迭代过程，其中内部的迭代为用式（27）来求解方程（23），而外部的迭代则为EM算法的迭代过程。迭代过程的终止条件可以设定为连续的参数估计值相差不超过某个设定的门限值或者达到设定的最大迭代次数。

Claims (3)

1.一种基于EM算法的G0分布参数最大似然估计方法，其特征在于，包括如下具体步骤：

步骤（1）：采用矩估计法，为观测得到的SAR图像生成参数估计初始值；

步骤（2）：根据步骤（1）中的参数估计初始值，采用迭代算法估计G0分布参数。

2.根据权利要求1所述的基于EM算法的G0分布参数最大似然估计方法，其特征在于：所述步骤（1）还包括如下具体步骤：

步骤（11）：设观测得到的SAR图像幅度数据为：y={y_t,t=1,…,T}，

且独立同分布，则y_t的G0分布表达式为：

$p\left(y_t\right|\mathrm{\λ})=\frac{2^{1+\mathrm{\β}}{y}_t^{N-1}{\mathrm{\σ}}^{\mathrm{\β}}}{\mathrm{\Γ}(N/2)\mathrm{\Γ}\left(\mathrm{\β}\right)}{(y_t^2+2\mathrm{\σ})}^{-(N/2+\mathrm{\β})}\mathrm{\Γ}(N/2+\mathrm{\β}),\mathrm{\β},\mathrm{\σ},N,y_t>0$

其中；N定义为等效视数，σ定义为形状参数，β定义为尺度参数，t定义为SAR图像幅度数据的序号，T定义为SAR图像幅度数据的个数；

步骤（12）：利用矩估计法来估计参数初始值，采用如下公式：

$\left\{\begin{array}{c}\widehat{\mathrm{\β}}=1+\frac{N{\hat{m}}_4}{N{\hat{m}}_4-(N+2){\hat{m}}_2^2}\\ \widehat{\mathrm{\σ}}=\frac{(\widehat{\mathrm{\β}}-1){\hat{m}}_2}{N}\end{array}\right.$

其中，

定义为k阶样本矩，即有

得到尺度参数估计初始值

和形状参数估计初始值

3.根据权利要求2所述的基于EM算法的G0分布参数最大似然估计方法，其特征在于，所述步骤（2）包括如下具体处理：

步骤（21）：令

建立如下函数变量A和G：

$A=\frac{1}{T}\underset{t=1}{\overset{T}{\mathrm{\Σ}}}\frac{N+2{\mathrm{\β}}^{\′}}{y_t^2+2{\mathrm{\σ}}^{\′}},$

$G={\left(\underset{t=1}{\overset{T}{\mathrm{\Π}}}\right[\frac{y_t^2+2{\mathrm{\σ}}^{\′}}{2}\left]\right)}^{1/T}\exp (-\mathrm{\Ψ}\left(\frac{N+2{\mathrm{\β}}^{\′}}{2}\right)),$

其中，Ψ(·)定义为digamma函数；

步骤（22）：令

利用步骤（21）中的函数变量A和G与公式

进行迭代，当满足迭代终止条件后，结束迭代过程，令参数估计值

等于最后一次迭代所得的则σ的估计值

由公式求得，并将迭代终止时的

和

作为最终的参数估计值；若不满足迭代终止条件，则返回步骤（21）。

Images