基于分量分离方法的复合高斯模型参数快速估计方法.pdf

基于分量分离方法的复合高斯模型参数快速估计方法，本发明涉及复合高斯模型参数快速估计方法。本发明是要解决目前复合高斯分布模型参数估计搜索效率低计算量大的问题，而提供了基于分量分离方法的复合高斯模型参数快速估计方法。第一步，利用相关时间滤波器分离复合高斯分布的调制分量和散斑分量；第二步，使用SISE方法估计调制分量的尺度参数b1与散斑分量的尺度参数b2；第三步，利用极大似然法估计出调制分量广义Gamma分布的规模参数a1和形状参数v1和散斑分量广义Gamma分布的规模参数a2和形状参数v2；第四步，将调制分量规模参数a1和散斑分量规模参数a2相乘得到复合高斯分布的规模参数a。本发明应用于雷达领域。

1.  基于分量分离方法的复合高斯模型参数快速估计方法，其特征在于它按以下步骤实现：
第一步，利用相关时间滤波器分离复合高斯分布的调制分量和散斑分量；
第二步，使用SISE方法估计调制分量的尺度参数b₁与散斑分量的尺度参数b₂；
第三步，利用极大似然法估计出调制分量广义Gamma分布的规模参数a₁和形状参数v₁和散斑分量广义Gamma分布的规模参数a₂和形状参数v₂；
第四步，将调制分量规模参数a₁和散斑分量规模参数a₂相乘得到复合高斯分布的规模参数a。

2.  根据权利要求1所述的基于分量分离方法的复合高斯模型参数快速估计方法，其特征在于步骤一具体为：
设某随机过程服从复合高斯分布，且散斑分量的自相关长度小于N₁，而调制分量在相邻N₁点区间内取值不变为s_N1i，其中i＝1,2,3,…，N₁个点应服从广义伽马分布，其PDF为：
fGΓ(x|sN1i)=b1sN1iΓ(v1)(xsN1i)b1v1-1·exp(-(xsN1i)b1)---(1)]]>
其中，b₁为调制分量尺度参数，v₁为调制分量形状参数，x为随机变量，其一阶矩为：
E(x)=sN1iΓ(1b1+v1)Γ(v1)---(2)]]>
令其中的常数项则
E(x)=sN1ia1---(3)]]>
由式(3)可知，把原数据列按每N₁个点分段，并对每段取均值，得到新数列：新数列服从如下广义Gamma分布：
fGΓ(s)=b2a2Γ(v2)(sa2)b2v2-1exp(-(sa2)b2)---(4)]]>
f_GΓ(s)是杂波自相关时间较长的调制分量，其中a₂＝a/a₁为散斑分量规模参数；b₂为 散斑分量的尺度参数，v₂为调制分量形状参数，若令每个数据除以其所在段内的均值则得到原数据中相关时间较短的散斑分量，其PDF为：
fGΓ(x)=b1a1Γ(v1)(xa1)b1v1-1exp(-(xa1)b1)---(5).]]>

3.  根据权利要求2所述的基于分量分离方法的复合高斯模型参数快速估计方法，其特征在于步骤二具体为：
广义Gamma分布的随机过程其PDF为：
fGΓ(x)=biaiΓ(v)(xai)bivi-1exp(-(xai)bi)---(6)]]>
x为随机变量，a_i为规模参数，b_i为尺度参数，v_i为形状参数，其中i＝1,2，；则其t阶原点矩为：
E(xit)=aitΓ(t/bi+vi)Γ(vi)---(7)]]>
设X是一个服从广义Gamma分布的随机过程，定义函数f(t),g(t)，其中f’(t)是f(t)的导数
f(t)=E(xit)---(8)]]>
g(t)=f′(t)f(t)---(9)]]>
将表达式(7-8)带入表达式(9)并整理得：
g(t)=ln(ai)+1biψ(tbi+vi)---(10)]]>
其中是digamma函数，由表达式(8)、表达式(10)得：
f(2vi)f(vi)=E[Xi2vi]E2[Xivi]=1+vivi---(11)]]>
g(bi)-g(0)=1vibi---(12)]]>
利用公式(7～12)构造如下方程：
F(z)=f(2z)f2(z)-{1+z[g(z)-g(0)]}=Γ(vi)Γ(2z/bi+vi)Γ2(z/bi)-{1+zbi[ψ(z/bi+v)-ψ(vi)]}=0---(13)]]>
该方程有唯一解z＝b_i，所以该方程被称为形状独立的尺度估计方程；为了实现估值需将实得到的有限的离散数据带入公式(13)中；
设实测数据个数为N，则随机过程原点矩可以用实测数据的原点矩估计，于是有
fN(t)=E(Xit)=1NΣk=1Nxikt---(14)]]>
gN(t)=fN′(t)/fN(t)=Σk=1NxiktlnxikΣk=1Nxikt---(15)]]>
其中x_ik为相应的分量实际得到的观测值，带入表达式(13)得到：
FN(z)=fN(2z)fN2(z)-(1+z[gN(z)-gN(0)])=0---(16)]]>
当N→∞时，F_N(z)→F(z)，所以把方程F_N(z)＝0的根b_i当做方程F(z)＝0的根b_i的估计值；由于F_N(z)在区间z∈(0，∞)上连续且方程F_N(z)＝0具有唯一解，采用数值方法求解该方程，即得到调制分量的尺度参数b₁与散斑分量的尺度参数b₂。

4.  根据权利要求3所述的基于分量分离方法的复合高斯模型参数快速估计方法，其特征在于步骤三具体为：
参数a_i、v_i的极大似然估计为：
ai=(1NviΣk=1Nxikbi)1/bi---(17)]]>
vi={bi[(Σk=1Nxikbi)-1Σk=1Nxikbilnxik-1NΣk=1Nlnxik]}-1---(18)]]>
由第二步求出的调制分量尺度参数b₁和散斑分量尺度参数b₂分别代替上式中算出调制分量规模参数a₁和形状参数v₁的最大似然估计值，散斑分量规模参数a₂和形状参数v₂的最大似然估计值。

5.  根据权利要求4所述的基于分量分离方法的复合高斯模型参数快速估计方法，其特征在于步骤四具体为：
将步骤三得到的调制分量规模参数a₁和散斑分量规模参数a₂的乘积即为整体的规模参 数a。

基于分量分离方法的复合高斯模型参数快速估计方法
技术领域
本发明涉及复合高斯模型参数快速估计方法。
背景技术
对雷达杂波进行建模是雷达杂波抑制和目标检测器设计的必要步骤，对杂波进行精确的建模可以有效地提高雷达的工作效率和探测能力。从雷达出现到如今已经建立了众多经典的分布：瑞利分布、对数正态分布、韦布尔分布等等，海杂波由于功率强、多普勒特性复杂严重影响对海面目标、掠海目标探测，所以如何对其进行有效的估计和抑制一直是雷达领域研究中重点之一。随着技术的进步，雷达的分辨率和频率越来越高，同时雷达对强杂波强干扰下的目标检测性能提出更高的要求，海杂波的统计特性变得更加复杂，经典模型的误差变得越来越大无法满足适用要求，近年来为解决这一问题，国内外研究人员做了大量探索工作，建立了复合高斯分布模型。复合高斯分布模型不仅可以在较宽的的频谱上拟合实际海杂波统计数据，还可以在杂波散射机理上得到很好的解释，但是由于该模型引用参数多，模型结构复杂，除了在一些特殊参数值情况下，这个模型的积分形式没有封闭形式的表达式(closed-form expression)，所以它的参数估计难度和计算量较大。目前复合高斯分布模型主要参数估计方法为基于矩函数估计的数值算法，如：最小二乘估计、粒子群智能搜索算法等，但这些方法存在的问题主要在于需要进行五维搜索计算量大处理时间长难以满足实时处理要求，同时对搜索的初值较为敏感，容易陷入局部最优解。
发明内容
本发明是要解决目前复合高斯分布模型参数估计搜索效率低计算量大的问题，而提供了基于分量分离方法的复合高斯模型参数快速估计方法。
基于分量分离方法的复合高斯模型参数快速估计方法，它按以下步骤实现：
第一步，利用相关时间滤波器分离复合高斯分布的调制分量和散斑分量；
第二步，使用SISE方法估计调制分量的尺度参数b₁与散斑分量的尺度参数b₂；
第三步，利用极大似然法估计出调制分量广义Gamma分布的规模参数a₁和形状参数v₁和散斑分量广义Gamma分布的规模参数a₂和形状参数v₂；
第四步，将调制分量规模参数a₁和散斑分量规模参数a₂相乘得到复合高斯分布的规模参数a。
发明效果：
欲达到的性能将复合高斯分布模型参数估计运算效率较目前已有算法提高十倍以上(即在保证精度优于或不低于目前已有算法的同时，将运算用时降低到现有算法的十分之一以下)。
一个复合高斯分布模型参数的快速估计算法。首先利用相关滤波器分离复合高斯分布的两个分量(调制分量和三斑分量)，然后利用独立形状参数方程(SISE)先计算出两个形状参数，最后把已经估计出的参数当做已知条件，利用最大似然估计估计出其余参数。同时本算法中分离出的两个分量计算流程相同所以可以采用分时复用的方法，分别计算两个分量的参数，最后再综合计算出剩余参数。
本发明结合分量分离和独立形状参数方程提出了一种新的参数估计方法，可以有效降低计算量，实现对复合高斯模型参数快速有效估计。
附图说明
图1是利用四种方法对IPIX雷达‘lo.dat’数据包数据进行估计拟合效果图；
图2是利用四种方法对雷达hi.dat数据包数据进行估计拟合效果图；
图3利用三种方法对IPIX雷达‘lo.dat’数据包数据进行估计拟合效果误差分析和计算用时图；
图4利用三种方法对IPIX雷达‘hi.dat’数据包数据进行估计拟合效果误差分析图；
图5本发明算法流程图；
图6(a)IPIX雷达的‘hi.dat’数据包实测数据的散斑分量的自相关函数图；
图6(b)IPIX雷达的‘hi.dat’数据包实测数据的调制分量的自相关函数图；
图7分量分离方法简要流程图；
图8按最右侧从上向下依次为v取0.2、0.4、…、2时构造出的函数F(x/b)的曲线图；
图9利用SISE方程二分法估计广义伽玛分布的参数b的流程图。
具体实施方式
具体实施方式一：本实施方式的基于分量分离方法的复合高斯模型参数快速估计方法，它按以下步骤实现：
第一步，利用相关时间滤波器分离复合高斯分布的调制分量和散斑分量；
第二步，使用SISE方法估计调制分量的尺度参数b₁与散斑分量的尺度参数b₂；
第三步，利用极大似然法估计出调制分量广义Gamma分布的规模参数a₁和形状参数v₁和散斑分量广义Gamma分布的规模参数a₂和形状参数v₂；
第四步，将调制分量规模参数a₁和散斑分量规模参数a₂相乘得到复合高斯分布的规 模参数a。
复合高斯分布模型把海杂波分为相关时间较短(微秒级)的散斑分量和相关时间比较长(秒级)的调制分量。并把两种分量都用广义Gamma分布模型来表示。
其散斑分量概率密度函数(Probability density function，PDF)分布为：
fGΓ(x|s)=b1sΓ(v1)[xs]b1v1-1exp[-(xs)b1]---(1)]]>
调制分量概率密度函数为：
fGΓ(s)=b2aΓ(v2)[sa]b2v2-1exp[-(sa)b2]---(2)]]>
复合高斯分布模型概率密度函数(CG-PDF)为：
fGC(x)=∫0∞fGΓ(x|s)fGΓ(s)ds---(3)]]>
把表达式(1)、表达式(2)代入表达式(3)并化简可得：
fGC(x;a,b1,b2,v1,v2)=b1b2Γ(v1)Γ(v2)xb1v1-1ab2v2∫0∞sb2v2-b1v1-1exp[-(sa)b2-(xs)b1]ds---(4)]]>
其中a为规模参数，b₁、b₂为功率参数，v_1、v₂为形状参数。
表达式(4)没有一般形式的解析表达式，但其各阶矩函数是有解析表达式的。其k阶矩函数可表示为：
E(xk)=akΓ(kb1+v1)Γ(kb2+v2)Γ(v1)Γ(v2)---(5)]]>
复合高斯分布模型的参数v₁、b₁由散斑分量单独决定，而b₂、v₂两个参数则由调制分量单独决定。由于调制分量的相关时间远大于散斑分量相关时间，所以可以通过相关滤波器将他们分开，然后分别对两个分量的参数进行估计，最后在综合估计出复合高斯分布模型的参数a。
具体实施方式二：本实施方式与具体实施方式一不同的是：步骤一具体为：
设某随机过程服从复合高斯分布，且散斑分量的自相关长度小于N₁，而调制分量在相邻N₁点区间内取值不变为s_N1i，其中i＝1,2,3,…，N₁个点应服从广义伽马分布，其PDF为：
fGΓ(x|sN1i)=b1sN1iΓ(v1)(xsN1i)b1v1-1·exp(-(xsN1i)b1)---(6)]]>
其中，b₁为调制分量尺度参数，v₁为调制分量形状参数，x为随机变量，其一阶矩为：
E(x)=sN1Γ(1b1+v1)Γ(v1)---(7)]]>
令其中的常数项则
E(x)=sN1a1---(8)]]>
由式(3)可知，对原数据做每相邻N₁点取均值的结果就是这N₁个数据内调制分量的值除以常数项，因此，把(原服从复合高斯分布的)原数据列按每N₁个点分段，并对每段取均值，得到新数列：新数列服从如下广义Gamma分布：
fGΓ(s)=b2a2Γ(v2)(sa2)b2v2-1exp(-(sa2)b2)---(9)]]>
f_GΓ(s)是杂波自相关时间较长的调制分量，其中a₂＝a/a₁为调制分量规模参数；b₂为调制分量的尺度参数，v₂为调制分量形状参数，若令每组N₁个点同时除以其组内的均值E(x)则得到的新的过程应服从分布：
fGΓ(x)=b1a1Γ(v1)(xa1)b1v1-1exp(-(xa1)b1)---(10)]]>
这个过程f_GΓ(x)就是某随机中相关时间较短的散斑分量。
可以看出要求解上面的几个公式难度依然较大，但是由表达式(10～12)可知，若是能够利用其它的方法先行估计出b的值，则可以较快地求出a和v的值。本文选用了基于SISE的方法实现对参数b的快速估计，具体过程将在第4节中推导。对散斑分量、调制分量分别利用上面的方法进行估算，在两者联合的复合高斯分布中，分量的参数b₁、b₂、v₁、v₂就是复合高斯分布的功率参数和形状参数，两个分量的尺度参数a₁、a₂的乘积即为整体的尺度参数的估计值。
其它步骤及参数与具体实施方式一相同。
具体实施方式三：本实施方式与具体实施方式一或二不同的是：步骤一中分离出的调制分量与散斑分量均为广义伽马分布，所以使用同样的方法求解参数b₁和b₂，系统实现时既可以用两个通道处理两个分量也可以通过分时复用处理他们。下面用b_i代指参数b₁和b₂，表述其实现方法。
设X是一个服从广义Gamma分布的随机过程其型概率密度函数为：
fGΓ(x)=biaiΓ(v)(xai)bivi-1exp(-(xai)bi)---(11)]]>
x为随机变量，a_i为规模参数，b_i为尺度参数，v_i为形状参数，其中i＝1,2，；则其t阶原点矩为：
E(Xit)=aitΓ(t/bi+vi)Γ(vi)---(12)]]>
定义函数f(t),g(t)，其中f’(t)是f(t)的导数
f(t)＝E(X^t)   (13)
g(t)=f′(t)f(t)---(14)]]>
将表达式(7-8)带入表达式(9)并整理得：
g(t)=ln(ai)+1biΨ(tbi+vi)---(15)]]>
其中是digamma函数，由表达式(8)、表达式(10)得：
f(2vi)f(vi)=E[Xi2vi]E2[Xivi]=1+vivi---(16)]]>
g(bi)-g(0)=1vibi---(17)]]>
利用公式(7～12)构造如下方程：
F(z)=f(2z)f2(z)-{1+z[g(z)-g(0)]}=Γ(vi)Γ(2z/bi+vi)Γ2(z/bi)-{1+zbi[Ψ(z/bi+v)-Ψ(vi)]}=0---(18)]]>
该方程有唯一解z＝b_i，而参数v的取值仅影响函数F(z)曲线的弯曲程度，与函数F(z)＝0的解无关。所以该方程被称为形状独立的尺度估计方程(SISE方程)；为了实现估值 需将实得到的有限的离散数据带入公式(13)中；
如图8，按最右侧从上向下依次为_i取0.2、0.4、…、2时F(x/b_i)的曲线，可以看出F(x/b_i)＝0方程的解总为1，参数v_i的取值仅影响函数F(z)曲线的弯曲程度
设实测数据个数为N，有
fN(t)=E(Xit)=1NΣk=1Nxikt---(19)]]>
gN(t)=fN′(t)/fN(t)=Σk=1NxiktlnxikΣk=1Nxikt---(20)]]>
带入表达式(13)得到：
FN(z)=fN(2z)fN2(z)-(1+z[gN(z)-gN(0)])=0---(21)]]>
当N→∞时，F_N(z)→F(z)，所以把方程F_N(z)＝0的根b_i当做方程F(z)＝0的根b_i的估计值；求解该方程的解析解较为复杂，但是由于F_N(z)在区间z∈(0，∞)上连续且方程F_N(z)＝0具有唯一解，采用数值方法求解该方程，即得到调制分量的尺度参数b₁与散斑分量的尺度参数b₂。本专利使用的是二分法求解，利用SISE方程求估计广义伽玛分布的参数b的流程图如图9。
其它步骤及参数与具体实施方式一或二相同。
具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：步骤三具体为：
广义Gamma分布的参数a、v的极大似然估计为：
ai=(1NviΣk=1Nxikbi)1/bi---(22)]]>
vi={bi[(Σk=1Nxikbi)-1Σk=1Nxikbilnxik-1NΣk=1Nlnxik]}-1---(23)]]>
由第二步求出的调制分量尺度参数b₁和散斑分量尺度参数b₂分别代替上式中的b算出调制分量规模参数a₁和形状参数v₁的最大似然估计值，散斑分量规模参数a₂和形状参数v₂的最大似然估计值。
其它步骤及参数与具体实施方式一至三之一相同。
具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：步骤四具体为：
将步骤三得到的调制分量规模参数a₁和散斑分量规模参数a₂的乘积即为整体的规模 参数a。
其它步骤及参数与具体实施方式一至四之一相同。
实施例：
实验中采用的两组实验数据是IPIX雷达1993年11月在加拿大Nova Scotia的Dartmouth进行实验所得，数据包的名字为‘hi.dat’和‘lo.dat’。这两组数据已经经过了预处理，长度都是131082*2，第一列为同相分量，第二列为正交分量。这两组数据分别为IPIX雷达在高海况、低海况情况下采集到的实测数据，它们已经被很多研究海杂波的论文采用，本发明也将利用它们来验证算法的性能。采集数据时的环境参数如表1所示：
表格1环境参数

本实验采用的计算机配置如下：cpu：Inter i3-m350@2.27GHz，内存：4G，操作系统：Windous7，计算软件：Matlab2013a。
结果如图1至图4、表2、表3所示。本实施例采用了三种常用的误差检验方法：均方差检验、K-S检验和CHI检验检验对比这本文提到的几种参数估计方法的无偏性，结果如表2、表3所示，可以看出采用本文的参数分离的估计方法估计效果优于复合高斯模型最小二乘估计算法和粒子群估计算法估计。分量分离的参数估计方法计算用时小于于最小二乘参数估计算法和粒子群估计算法用时的十分之一，即可以达到最初设定的目标性能。
表格2地海况lo.dat(误差小、用时短则性能好)

表格3地海况hi.dat(误差小、用时短则性能好)

图1中利用四种方法对IPIX雷达‘lo.dat’数据包数据进行估计拟合效果图,分量分离算法拟合效果较好,最小二乘算法估计波峰峰值偏大且偏右,而粒子群算法分布太分散；
图2.利用四种方法对雷达hi.dat数据包数据进行估计拟合效果图,此时几种方法估计结果差距较小，但仍可以看出分两分离算法同原数据PDF曲线接近程度最高
图3利用三种方法对IPIX雷达‘lo.dat’数据包数据进行估计拟合效果误差分析和计算用时图,为方便查看图中结果为化粒子群算法归一化后的相对误差，误差数值越小，效果越好；计算用时越小，性能越好。可以看出分量分离法误差最小。
图4利用三种方法对IPIX雷达‘hi.dat’数据包数据进行估计拟合效果误差分析图,方便查看图中结果为化粒子群算法归一化后的相对误差，数值越小，效果越好。粒子群算法的均方误差略小于分两分离法，K-S检验两者结果相当，而CHI检验中分两分离法误差小于粒子群算法，此时精度上两个方法相当，而CHI检验中分两分离法误差小于粒子群算法，此时精度上两个方法相当，但考虑到分量分离法计算用时远小于粒子群算法，而且在其他组数据中分量分离算法估计精度要优于粒子群算法，仍可认为分量分离算法性能较好。
图6IPIX雷达的‘hi.dat’数据包实测数据的散斑分量、调制分量的自相关函数。可以看出散斑分量自相关长度只有几十毫米，且随时间差变大迅速下降；而调制分量的自相关时间长达数十秒，且呈现震荡衰减特性。
图8按最右侧从上向下依次为v_i取0.2、0.4、…、2时F(x/b_i)的曲线，可以看出F(x/b_i)＝0方程的解总为1，参数v_i的取值仅影响函数F(z)曲线的弯曲程度。
图9利用SISE方程二分法估计广义伽玛分布的参数b_i的流程图。