一种海杂波K分布加噪声模型参数的快速分数阶矩估计方法与流程

本发明涉及信号处理技术领域，尤其涉及一种海杂波k分布加噪声模型参数的快速分数阶矩估计方法。

背景技术：

海杂波是由大量相互独立的海面散射体的后向散射相互叠加形成的回波，受雷达参数，海面无规律运动的影响，海杂波表现出明显的高幅值、非平稳、非高斯性。k分布模型是能够有效描述海杂波特性的经典统计模型，不仅在很宽的条件范围内可以与海杂波幅度分布很好的进行匹配，还可以正确地描述杂波的时间和空间的相关性，被广泛运用于杂波仿真、目标检测等领域。然而，实测海杂波数据中还存在雷达设备本身产生的加性噪声，对k分布模型的拟合效果有一定影响，为此又提出了k分布加噪声模型。

在杂波背景下的自适应检测中，检测门限值往往与分布模型的各项参数有关。为了在一定恒虚警率下获得较高的检测概率，需要准确地估计k分布加噪声模型的形状参数、尺度参数和噪声功率。目前对k分布加噪声模型的形状参数的估计方法主要有矩估计方法。因为k分布加噪声的矩表达式中含有广义超几何函数，所以不是任意两个不同的矩都可以估计出三个参数的值，但可以运用任意三个偶数阶矩求得。因此，使用三个偶数阶矩估计时，由于阶数较大，数据的计算量会非常高，形状参数估计速度比较慢，并且受异常散射单元影响比较大，稳健性不足。

技术实现要素：

本发明要解决的技术问题是：现有k分布加噪声模型形状参数ν估计速度慢、稳健性不足的问题。

为了解决上述技术问题，本发明提供了如下技术方案：

一种海杂波k分布加噪声模型参数的快速分数阶矩估计方法，包括如下步骤：

s1：从1.0开始，间隔0.1取一个值，取到100.0为止，生成一组序列ν1,ν2,···,νi,···νl作为形状参数ν取值的序列，l是序列中数据的个数，取值为991；

s2：从1.0开始，间隔0.1取一个值，取到100.0为止，生成一组序列μ1,μ2,···,μj,···μs作为数据依赖项q1取值的序列，s是序列中数据的个数，取值为991，q1＝m2/δ²，m2是海杂波样本数据的二阶矩，δ是海杂波样本数据的噪声功率；

s3：获得与形状参数ν的序列和数据依赖项q1的序列对应的数据依赖项q2的序列λ1,λ2,···,λ(j-1)l+i,···λl×s，其中，mr是海杂波样本数据的r阶矩，r是分数阶距的阶数，取值0.1；

s4：将两个数据依赖项序列与形状参数的倒数序列进行非线性拟合，得到如下的近似关系式：

式中，p1，p2，p3，p4，p5，p6，p7，p8为常数项；

s5：利用周期图法获得海杂波样本数据噪声功率δ的估计值

s6：利用海杂波数据的样本序列，结合q1、q2的计算公式以及获得q1的估计值和q2的估计值

s7：将取代近似关系式中的q1，将取代近似关系式中的q2，获得形状参数ν的估计值

优选地，p1，p2，p3，p4，p5，p6，p7，p8的取值分别为-0.70536547，-2.43284673，0.44779562，-0.02296508，-1.86849359，1.45263166，1.67059762，0.06325657。

优选地，步骤s3包括如下步骤：

s31：根据全概率公式，得到模型的r阶矩表达式：

式中，γ(·)是gamma函数，2f0(·)是第二类合流超几何函数，σ是模型的尺度参数；

s32：利用模型的二阶矩表达式获得形状参数与尺度参数的关系式：

s33：结合s31获得的公式和s32获得的公式，得到如下公式：

s34：根据q1和q2的公式，结合s33获得的公式，得到：

s35：利用s34得到的计算方法获得与序列ν1,ν2,···,νi,···νl以及序列μ1,μ2,···,μj,···μs对应的序列λ1,λ2,···,λ(j-1)l+i,···λl×s。

优选地，步骤s4包括如下步骤：

s41：将形状参数的序列用第一个行向量s表示；将数据依赖项q1的序列用第二个行向量u表示，将数据依赖项q2的序列用第个三行向量v表示；

s42：将三个行向量作为matlab中提供的工具包cftool的参数，调用工具包非线性拟合功能，得到形状参数与数据依赖项的近似关系式。

优选地，步骤s5包括如下步骤：

s51：雷达发射机发射多个脉冲信号，雷达接收机接收经过海面散射形成的去除含目标单元的纯杂波幅度数据xr，作为海杂波数据的样本序列x1,x2,...,xj,...,xm×n，xj是参考序列的第j个数据，j＝1,2,...,m×n，m是样本序列数据的距离单元数，n是样本序列数据的脉冲数，该数据的距离单元集合为：ui＝{ui(n),n＝1,2,...,n},i＝1,2,...,m，ui表示第i个距离单元的样本集合，ui(n)表示第i个距离单元第n个脉冲的数据；

s52：对每个距离单元计算其功率谱，得到功率谱密度序列p1,p2,...,pi,...,pk，k是傅里叶变换的采样点数，将功率谱密度序列递增排序，得到功率谱密度增序序列p(1),p(2),...,p(i),...,p(k)，利用功率谱密度增序序列的前半部分数据估计出第i个距离单元杂波的噪声功率：

优选地，在步骤s52中，根据如下公式获得所述功率谱：

优选地，步骤s6包括如下步骤：

s61：利用海杂波数据的样本序列，获得二阶矩的估计值

s62：利用海杂波数据的样本序列，获得0.1阶矩的估计值

s63：利用二阶矩、0.1阶矩、噪声功率的估计值获得和

优选地，还包括s8：利用获得的形状参数的估计值获得海杂波样本数据的尺度参数σ的估计值

优选地，在步骤s1中，从1.0开始，间隔0.1取一个值，取到100.0为止，生成序列ν1,ν2,···,νi,···νl，l是序列中数据的个数，取值为991。

优选地，在步骤s2中，从1.0开始，间隔0.1取一个值，取到100.0为止，生成序列μ1,μ2,···,μj,···μs，s是序列中数据的个数，取值为991。

有益效果

本发明的上述技术方案具有如下优点：

由于本发明使用海杂波k分布加噪声模型的分数阶矩、二阶矩与形状参数的近似数学表达式来估计形状参数，与现有的矩估计方法相比，本发明计算更加简单，求解更加快速。

由于本发明使用较低阶数的矩来估计形状参数，可以有效减少参数估计需要的样本数量，增加估计精度；同时减少了样本中较大数值在估计中占的比重，避免了异常单元的影响，得到参数的稳健估计。

另外，利用本发明提供的方法可以快速地获得k分布加噪声模型的尺度参数的估计值。

附图说明

图1是本发明的实现流程图；

图2是用本发明和现有矩估计方法对仿真数据中的k-s曲线的对比图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合本发明实施例，对本发明的技术方案进行清楚、完整地描述。显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

图1给出了根据本发明的一种实施方式的海杂波k分布加噪声模型(该模型为现有理论，如，watts，s.提出的理论模型)参数的快速分数阶矩估计方法的流程示意图。如图1所示，根据本发明的一种实施方式，海杂波k分布加噪声模型参数的快速分数阶矩估计方法的具体步骤为：

步骤1：生成一组序列ν1,ν2,···,νi,···νl作为形状参数ν取值的序列。

从1.0开始，间隔0.1取一个值，取到100.0为止，得到海杂波k分布加噪声模型形状参数ν取值的序列ν1,ν2,···,νi,···νl。νi是序列的第i个数据，i＝1,2,3,···，l。l是序列中数据的个数，取值为991。

步骤2：生成一组序列μ1,μ2,···,μj,···μs作为数据依赖项q1取值的序列。

所述数据依赖项q1＝m2/δ²，该公式中，m2是海杂波样本数据的二阶矩，δ是海杂波样本数据的噪声功率。

从1.0开始，间隔0.1取一个值，取到100.0为止，得到海杂波k分布加噪声模型数据依赖项q1的序列μ1,μ2,···,μj,···μs。μj是序列的第j个数据，j＝1,2,3,···，s。s是序列中数据的个数，取值为991。

步骤3：求与海杂波k分布加噪声模型形状参数ν的序列和数据依赖项q1的序列对应的数据依赖项q2的序列λ1,λ2,···,λ(j-1)l+i,···λl×s，λ(j-1)l+i是数据依赖项q2序列的第(j-1)l+i个数据，是νi与μj对应的值。

所述数据依赖项该公式中，mr是海杂波样本数据的r阶矩，r是分数阶距的阶数，本发明取值0.1。

具体地，该步骤包括：

步骤31：根据全概率公式，计算得到海杂波k分布加噪声模型的r阶矩表达式：

式中，γ(·)是gamma函数，2f0(·)是第二类合流超几何函数，σ是模型的尺度参数；

步骤32：利用海杂波k分布加噪声模型的二阶矩表达式，求k分布加噪声模型形状参数与尺度参数的关系式，得到如下公式：

步骤33：将公式3-2代入公式3-1，整理后得到：

步骤34：将数据依赖项q1和q2代入公式3-3得到：

步骤35：利用公式3-4得到的计算方法求与海杂波k分布加噪声模型形状参数ν的序列ν1,ν2,···,νi,···νl以及数据依赖项q1的序列μ1,μ2,···,μj,···μs对应的数据依赖项q2的序列λ1,λ2,···,λ(j-1)l+i,···λl×s。

步骤4：将两个数据依赖项序列与形状参数的倒数序列进行非线性拟合，得到k分布加噪声模型的数据依赖项q1，q2与形状参数ν的近似关系式。

具体地，该步骤包括：

s41：将形状参数ν的序列ν1,ν2,···,νi,···νl用第一个行向量s表示；将数据依赖项q1的序列μ1,μ2,···,μj,···μs用第二个行向量u表示，将数据依赖项的q2的序列λ1,λ2,···,λ(j-1)l+i,···λl×s用第个三行向量v表示。

s42：将三个行向量s、u和v作为matlab中提供的工具包cftool的参数，调用工具包非线性拟合功能，得到形状参数与两个数据依赖项的近似关系式：

式中，p1，p2，p3，p4，p5，p6，p7，p8为常数项，取值分别为-0.70536547，-2.43284673，0.44779562，-0.02296508，-1.86849359，1.45263166，1.67059762，0.06325657。

步骤5：利用周期图法获得海杂波样本数据噪声功率δ的估计值

具体地，该步骤包括：

步骤51：雷达发射机发射多个脉冲信号，雷达接收机接收经过海面散射形成的去除含目标单元的纯杂波幅度数据xr，作为海杂波数据的样本序列x1,x2,...,xj,...,xm×n，xj是参考序列的第j个数据，j＝1,2,...,m×n，m是样本序列数据的距离单元数，n是样本序列数据的脉冲数，该数据的距离单元集合为：

ui＝{ui(n),n＝1,2,...,n},i＝1,2,...,m，

ui表示第i个距离单元的样本集合，ui(n)表示第i个距离单元第n个脉冲的数据；

s52：对每个距离单元计算其功率谱：

得到功率谱密度序列p1,p2,...,pi,...,pk，k是傅里叶变换的采样点数。

将功率谱密度序列递增排序，得到功率谱密度增序序列p(1),p(2),...,p(i),...,p(k)。

利用功率谱密度增序序列的前半部分数据估计出第i个距离单元杂波的噪声功率：

[k/2]表示对k/2取整。

步骤6：利用海杂波数据的样本序列，结合q1、q2的计算公式以及获得q1的估计值和q2的估计值

具体地，该步骤包括：

步骤61：利用海杂波数据的样本序列，求二阶矩的估计值

步骤62：利用海杂波数据的样本序列，求0.1阶矩的估计值

步骤63：利用二阶矩、0.1阶矩、噪声功率的估计值求和

步骤7：将取代近似关系式4-1中的q1，将取代近似关系式4-1中的q2，获得形状参数ν的估计值

还可以利用获得的形状参数的估计值获得海杂波样本数据的尺度参数σ的估计值

下面结合仿真实验对本发明的效果做进一步说明。

1.仿真参数

仿真实验采用的数据为随机产生的60000个单元的服从形状参数从1到10，步长为0.5，尺度参数为1，加入杂噪比为5db的高斯白噪声的复合k分布的随机序列。每个形状参数值产生5000组数据，共产生95000组数据。

2.仿真实验内容

分别采用本发明方法和现有的矩估计方法得到海杂波k分布模型参数的估计值，通过k-s检验法，计算两种方法计算得到的模型累积分布函数与实测海杂波数据的经验累积分布函数在不同形状参数下的k-s距离；

然后，对所有数据得到的k-s距离按不同形状参数求均值，得到两条与本发明和矩估计方法对应的k-s距离曲线，结果如图2所示。

仿真实验

首先，在ipix雷达实测海杂波数据中文件名为19980223_215110_antstep.cdf的数据中选取20个不含目标的距离单元的纯海杂波数据。

然后，对些数据使用本发明方法和矩估计方法进行参数估计，计算得到本发明和矩估计方法对应的k-s距离。

最后，根据目标信号信杂比的变化，得到两条与本发明和矩估计方法对应的检测性能曲线，结果如图2所示。

从图2可以看出，本发明在提高形状参数估计速度的前提下检测性能方面与矩估计方法比较接近。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。