基于海浪谱模型的海杂波模拟方法

摘要

本发明公开了一种基于海浪谱模型的海杂波模拟方法，主要解决现有技术难以体现海表面浪高细节的问题。其实现方案是：1.建立基于海浪谱模型的海表面；2.确定雷达主波束照射下的海表面范围；3.在海表面范围内选取矩形框，并对框内海表面离散化构造成由许多等面积的三角形散射单元拼接而成的复杂表面；4.计算每个三角形散射单元的属性参数及回波功率；5.根据回波功率建立三角形散射单元的回波模型；6.按照回波模型产生矩形框内海表面三角形散射单元的回波，判断三角形散射单元的雷达可见性，计算并叠加所有可见的三角形散射单元的回波，产生海杂波的仿真数据并分析。本发明提高了海杂波模拟的可信度与实用性，可用于目标检测。

权利要求

1.基于海浪谱模型的海杂波模拟方法，包括：

(1)根据随机海浪理论，用波高方程η(x,y,t)描述基准水平面上任意一点(x,y)在t时刻的海表面高度，建立海表面的海浪模型，其中波高方程中的振幅由海浪谱函数确定；

(2)建立天线高度为H，波束宽度为2ψ，波束瞄准线入射余角为φ的下视雷达，根据天线主波束照射下的海表面上任意一点(x,y,z)与下视雷达的位置关系，计算出该点运动轨迹x2cos2ψ-2yHcosφsinφ+y2(cos2ψ-cos2φ)＝H2(sin2φ-cos2ψ)，确定海表面的范围：当φ＜ψ时，海表面范围为双曲线及其内部海表面；当φ＝ψ时，海表面范围为抛物线及其内部海表面；当φ＞ψ时，海表面范围为椭圆及其内部海表面；

(3)在确定的海表面范围，选取大小为p×q的一个矩形框，在该矩形框内横向采样间隔取 以该矩形框内海表面任意点(x0,y0,z0)的运动速度v，方位角α0，俯仰角β0，在保持多普勒分辨率与满足多普勒不模糊的条件下计算出纵向采样间隔sδ，并以横纵间隔分别为与sδ在矩形框内的海表面上取采样点，再以每3个采样点构成一个三角形散射单元，其中c为光速，fs为雷达信号采样频率，s为波门起始在y轴投影的位置，δ为方位角分割间隔，N为一个相参处理间隔CPI内的脉冲数，p为波门的起止在y轴投影的长度，q为确定的海表面范围的纵向宽度；

(4)计算任意一个三角形散射单元的参数：

(4a)根据任意一个三角形散射单元与雷达的位置关系，计算该三角形散射单元与雷达之间的距离R，方位角α与俯仰角β，来波方向 法向 面积A0，其中三角形散射单元的面积A0，由 计算得出；

(4b)根据雷达方程，计算该三角形散射单元的回波功率Pr；

(5)建立步骤(4)中三角形散射单元的回波模型 其中

st(·)为雷达发射的线性调频LFM信号，ζ(t)为相位延迟抖动，表示虚数，fd为该三角形散射单元的多普勒频率；

(6)根据射线追踪原则，判断三角形散射单元之间的遮挡关系和三角形散射单元自身的遮挡关系：设某三角形散射单元顶点相对于雷达的仰角为α1，遮挡单元顶点相对于雷达的仰角为α2，当α1≥α2时，该三角形散射单元被遮挡，不计算其回波；当α1＜α2时，该三角形散射单元没被遮挡单元遮挡，则计算该三角形散射单元的局部入射余角αL，再判断该三角形散射单元自身的遮挡关系：当αL＞π/2时，该三角形散射单元被自身遮挡，不计算其回波；

当αL≤π/2时，该三角形散射单元未被遮挡，计算并叠加所有未被遮挡的三角形散射单元的回波，产生一个相参处理间隔CPI回波的海杂波仿真数据，完成对海杂波的模拟。

2.根据权利要求1所述的方法，其中步骤(1)中建立海表面的海浪模型，按如下步骤进行：(1a)根据随机海浪理论，得到三维随机海浪波高方程的数学描述

其中m代表角频率分割数，M代表角频率

的最大分割数，l为方向分割数，L代表方向的最大分割数；km为海面波浪的波数，ωm为海面波浪的角频率，εm,l为随机相位角，Φl为海浪相对于x轴的传播方向；Am,l为不同角频率和方向的简单正弦波的振幅， Δωm和ΔΦl分别表示ωm和Φl的增量，S(ωm,Φl)为不同频率和方向的海浪能量分布，S(ωm,Φl)＝SΩ(ωm)D(Φl)，SΩ(ωm)为海浪频率谱函数，这里选择Pierson-Moskowitz海浪谱模型： g为重力加速度，Uw为高于海表面19.5m处的风速；D(Φl)为方向函数，采用国际船模试验水池会议ITTC提出的建议性形式：(1b)将海表面上任意一点(x,y,z)投影于基准水平面(x,y)点，用波高方程η(x,y,t)描述点(x,y)在t时刻的海表面高度z，建立海表面的海浪模型为：z＝η(x,y,t)。

3.根据权利要求1所述的方法，其中步骤(4)计算该三角形散射单元的回波功率Pr，按如下公式计算：其中Pt为雷达发射功率，λ为信号波长，Ls为雷达收发综合损耗，G(αF,βF；α,β)为波束瞄准线指向(αF,βF)时(α,β)方向上的天线增益，σ为三角形散射单元的雷达散射截面积，其由σ＝ΓHHA0或σ＝ΓVVA0计算，ΓHH为水平极化下的海面散射系数，ΓVV为垂直极化下的海面散射系数， ε为复介

电常数，αL为三角形散射单元的局部入射余角，其由该三角形散射单元来波方向 与法向确定：

说明书

基于海浪谱模型的海杂波模拟方法

技术领域

[0001] 本发明属于信号处理技术领域，具体地说是一种海杂波模拟方法，可用于目标检测。

背景技术

[0002] 雷达海杂波的仿真技术，经历了由浅到深、由简到繁、又粗到细的过程。最初是将海杂波看成具有一定概率密度分布函数PDF的随机序列，采用Monte-Carlo法进行仿真。经典的海杂波幅度统计模型有瑞利分布、对数正太分布、韦布尔分布以及K分布。而实际遇到的海杂波不仅具有幅度上的统计分布特性，同时又满足某种相关性，这就使得仿真的数据要同时满足某种特定的功率谱特性或相关特性和幅度分布特性，比较有代表性且成熟的方法有球不变随机过程法SIRP和零记忆非线性变换法ZMNL，还有一种极少用的随机微分方程法SDE。然而这些从统计理论出发的海杂波幅度统计模型是一种经验模型，未能考虑到雷达的实际工作条件及环境条件，不能直接运用于实际雷达系统的杂波分析。

[0003] 所以需要针对具体的雷达工作场景进行杂波模拟来研究模型参数与雷达参数及环境参数的关系。文献《任意姿态记载PD雷达三维地杂波算法研究》将地海面按距离-多普勒进行了散射单元划分，对机载脉冲多普勒雷达采用下视工作方式时的地海面杂波进行了建模与模拟；文献《PD雷达导引头海杂波谱的建模与仿真》结合实测数据修正后的海面散射系数等模型，将海面按等距离环进行了散射单元划分，对工作在脉冲多普勒体制下的雷达导引头海杂波进行了建模与仿真。前者根据距离-多普勒的散射单元划分方法在理论上很精确，但是采用数学公式表示这些散射单元比较复杂，且计算杂波功率需计算双重积分，实现起来比较困难；后者根据等距离环的散射单元划分方法比较简单，但散射单元面积随着距离的增加而增加，海表面细节难以体现，降低了海杂波模拟的可信度。

发明内容

[0004] 本发明的目的在于针对上述已有技术的不足，提出一种基于海浪谱模型的海杂波模拟方法，提高海杂波模拟的可信度与实用性。

[0005] 为实现上述目的，本发明技术方案包括如下：

[0006] (1)根据随机海浪理论，用波高方程η(x,y,t)描述基准水平面上任意一点(x,y)在t时刻的海表面高度，建立海表面的海浪模型，其中波高方程中的振幅由海浪谱函数确定；

[0007] (2)建立天线高度为H，波束宽度为2ψ，波束瞄准线入射余角为φ的下视雷达，根据天线主波束照射下的海表面上任意一点(x,y,z)与下视雷达的位置关系，计算出该点运动轨迹x2cos2ψ-2yHcosφsinφ+y2(cos2ψ-cos2φ)＝H2(sin2φ-cos2ψ)，确定海表面的范围：当φ＜ψ时，海表面范围为双曲线及其内部海表面；当φ＝ψ时，海表面范围为抛物线及其内部海表面；当φ＞ψ时，海表面范围为椭圆及其内部海表面；

[0008] (3)在确定的海表面范围，选取大小为p×q的一个矩形框，在该矩形框内横向采样间隔取 以该矩形框内海表面任意点(x0,y0,z0)的运动速度v，方位角α0，俯仰角β0，在保持多普勒分辨率与满足多普勒不模糊的条件下计算出纵向采样间隔sδ，并以横纵间隔分别为 与sδ在矩形框内的海表面上取采样点，再以每3个采样点构成一个三角形散射单元，其中c为光速，fs为雷达信号采样频率，s为波门起始在y轴投影的位置，δ为方位角分割间隔， N为一个相参处理间隔CPI内的脉冲数，p为波门的起止在y轴投影的长度，q为确定的海表面范围的纵向宽度；

[0009] (4)计算任意一个三角形散射单元的参数：

[0010] (4a)根据任意一个三角形散射单元与雷达的位置关系，计算该三角形散射单元与雷达之间的距离R，方位角α与俯仰角β，来波方向 法向 面积A0，其中三角形散射单元的面积A0，由 计算得出；

[0011] (4b)根据雷达方程，计算该三角形散射单元的回波功率Pr；

[0012] (5)建立步骤(4)中三角形散射单元的回波模型 其中st(·)为雷达发射的线性调频LFM信号，ζ(t)为相位延迟抖动，表示虚数，fd为该三角形散射单元的多普勒频率；

[0013] (6)根据射线追踪原则，判断三角形散射单元之间的遮挡关系和三角形散射单元自身的遮挡关系：

[0014] 设某三角形散射单元顶点相对于雷达的仰角为α1，遮挡单元顶点相对于雷达的仰角为α2，当α1≥α2时，该三角形散射单元被遮挡，不计算其回波；当α1＜α2时，该三角形散射单元没被遮挡单元遮挡，则计算该三角形散射单元的局部入射角αL，再判断该三角形散射单元自身的遮挡关系：

[0015] 当αL＞π/2时，该三角形散射单元被自身遮挡，不计算其回波；

[0016] 当αL≤π/2时，该三角形散射单元未被遮挡，计算并叠加所有未被遮挡的三角形散射单元的回波，产生一个相参处理间隔CPI回波的海杂波仿真数据，完成对海杂波的模拟。

[0017] 本发明与现有技术相比具有如下优点：

[0018] 1.本发明运用海浪谱模型与雷达工作参数、环境参数结合的方法仿真的海表面增强了模拟的海杂波与海况环境的联系，更加具有实用性；

[0019] 2.由于离散海表面可看成由许多三角形散射单元拼接而成的复杂表面，且三角形散射单元在水平面的投影为直角三角形面元，本发明对海平面采用三角形散射单元进行划分，相比以往基于等距离环的分辨单元划分方法，三角形面元属于等面积划分，划分的面元面积不会随距离的增加而增加，因此模拟的海表面更能体现出浪高细节。同时该划分方法在数学表示上比距离-多普勒分辨单元划分方法简单，不需要复杂的双重积分运算；

[0020] 实验仿真表明：1)本发明单次脉冲回波的海杂波仿真数据的幅度曲线图与海杂波实测数据的幅度曲线图幅度趋势基本一致，说明本发明产生的海杂波具有可信度；2)本发明同一距离单元的海杂波仿真数据的幅度统计分布与传统统计分布模型K分布拟合优度检验最佳，表明本发明的海杂波模拟方法具有实用性，提供研究统计分布模型参数与雷达工作参数及环境参数之间关系的一种途径。

附图说明

[0021] 图1是本发明的操作流程图；

[0022] 图2是本发明中使用的基于海浪的功率密度能量等分法示意图；

[0023] 图3是本发明中使用的在雷达主波束照射下确定海表面范围示意图；

[0024] 图4是本发明中使用的在矩形框内的海表面取采样点示意图；

[0025] 图5是本发明中使用的构造三角形散射单元示意图；

[0026] 图6是本发明中使用的计算三角形散射单元具体参数的三角形散射单元与雷达之间位置关系图；

[0027] 图7是本发明仿真的不同时刻的海表面二维图；

[0028] 图8是本发明的海杂波仿真数据幅度曲线图与海杂波实测数据幅度曲线图的对比图；

[0029] 图9是本发明仿真的海杂波回波数据分析图；具体实施方式：

[0030] 本发明结合基于风速的Pierson-Moskowitz海浪谱模型仿真海表面，根据各类对海平面划分散射单元方法的优缺点，提出等面积的三角形散射单元划分方法，计算并叠加所有散射单元的回波仿真海杂波，仿真时还考虑到了雷达工作参数与环境参数，如天线高度、波束宽度与波束瞄准线入射余角等因素。

[0031] 下面结合附图和具体实施例对本发明作进一步详细描述。

[0032] 参照图1，本发明的基于海浪谱模型的海杂波模拟方法，包括如下步骤：

[0033] 步骤1，建立海表面的海浪模型。

[0034] 由于风速等影响，海面会发生高低起伏的变化，故可知海表面不可能一直都是平整的表面，因此可根据随机海浪理论，将实际海浪看成为不同角频率、不同传播方向、不同波高和不同初始相位的正弦波叠加的结果。

[0035] 1a)选择海浪频率谱函数SΩ(ω)和方向函数D(Φ)：

[0036] 采用国际船模试验水池会议ITTC提出的建议性形式为：

[0037] 海浪频率谱函数SΩ(ω)选择Pierson-Moskowitz海浪谱模型，Pierson-Moskowitz谱是以风速为参量的充分成长状态的海浪频谱，其表达式为：

[0038]

[0039] 其中ω为海面波浪的角频率，Φ为海浪相对于x轴的传播方向，g为重力加速度，Uw为高于海表面19.5m处的风速；

[0040] 1b)根据海浪频率谱函数SΩ(ω)和方向函数D(Φ)，计算海浪谱函数S(ω,Φ)，[0041] S(ω,Φ)＝SΩ(ω)D(Φ)；

[0042] 1c)计算不同角频率和方向的简单正弦波的振幅： 其中Δω和ΔΦ分别为ω和Φ的增量；

[0043] 1d)根据上述参数建立如下三维随机海浪波高方程：

[0044]

[0045] 其中Am,l表示不同角频率和方向的正弦波振幅；km为波数，根据深水的色散关系，海面波浪的角频率 且 kmx＝km cosΦl，kmy＝km sinΦl；εm,l为随机相位角，服从[0,2π)内的均匀分布；l代表方向分割数，L代表方向的最大分割数，其取值越大分割精度越高，l＝1,2,...，L；m代表角频率分割数，M代表角频率的最大分割数，m＝1,
2,...，M；

[0046] 对角频率的分割参照图2，由于角频率可由1a)中Pierson-Moskowitz频率谱函数SΩ(ω)描述，故可采用基于海浪功率密度的能量等分法对角频率进行分割：首先对选定的角频率m次分割，使各角频率间隔的能量ΔE相等，然后确定各子面积的中心频率；

[0047] 定义角频率积累谱为：

[0048] 按照能量等分的概念有：

[0049] 其中b为常数，Mb＝E(∞)，mb为各子面积的中心频率，角频率被分割为ωm＝mb；

[0050] 1e)根据海表面某点(x,y,z)在基准水平面上的投影的坐标为(x,y)，则用η(x,y,t)表示t时刻点(x,y)处的海表面高度，因此可得出点(x,y,z)处海表面的海浪模型为：z＝η(x,y,t)。

[0051] 步骤2，确定天线主波束照射下海表面的范围。

[0052] 雷达照射海表面发生后向散射的主要部分取决于天线高度、波束宽度和波束瞄准线入射余角这3个因素。

[0053] 2a)建立一个建立天线高度为H，波束宽度为2ψ，波束瞄准线入射余角为φ的下视雷达，天线主波束指向y轴方向，波束瞄准线与y轴的交点为点C，则雷达波束瞄准线方向为：其中，Ra为雷达位置， i、j、k分别为三个坐标轴上的单位向
量；

[0054] 2b)参照图3，P(x,y,0)为基准水平面上雷达脚印中的任意一点，则雷达指向该点的方向向量为： 设雷达波束为针状波束或锥形波束，得到向量方程：

[0055]

[0056] 其中 表示向量 与向量 的内积，|·|表示计算一个向量的模值；

[0057] 2c)根据上述向量方程，计算点P(x,y,0)的轨迹：

[0058] Lycosφ+LHsinφ＝(x2+y2+H2)1/2Lcosψ，

[0059] 整理后可得到：

[0060] x2cos2ψ-2yHcosφsinφ+y2(cos2ψ-cos2φ)＝H2(sin2φ-cos2ψ)，

[0061] 对于给定的雷达高度H，给定的波束瞄准线入射余角φ和给定的波束宽度2ψ，当φ＞ψ时，该点的轨迹为双曲线；当φ＝ψ时，该点的轨迹为抛物线；当φ＞ψ时，该点的轨迹为椭圆，该点的轨迹及其内部包含的海表面即为确定的海表面范围。

[0062] 步骤3，在确定的海表面范围内，构造三角形散射单元。

[0063] 雷达信号处理机对上述确定的海表面采样时，只接收波门内的回波信号，所以要选取一个矩形框，在矩形框区域内的海表面上取采样点，构造散射单元。

[0064] 3a)选取矩形框：

[0065] 参照图4，选取大小为p×q的一个矩形框，p为波门的起止在y轴投影的长度，q为确定的海表面范围的纵向宽度，然后在该矩形框范围中的海表面上取采样点；其中可能存在一些海表面不在确定的海表面范围内却在矩形框中，由于天线波束会对这部分回波进行加权削弱，这部分海表面不会影响仿真。

[0066] 3b)计算采样点的横向间隔：

[0067] 设雷达的坐标为(xR,yR,zR)，海表面任意一点的坐标为(x0,y0,z0)，天线主波束指向y轴方向，雷达信号处理机在基准水平面上沿y轴方向等间距采样，计算采样点的横向间隔为：

[0068]

[0069] 其中c为光速，fs为采样频率，φ0为该点处的入射余角，φs＜φ0＜φe，φs、φe分别为波门起始处、波门结束处的入射余角，

[0070] 点(x0,y0,z0)离雷达的距离R0为：

[0071]

[0072] 入射余角为：

[0073] 随着y0增大，该点离雷达距离R0增大，该点处入射余角φ0减小，ΔR减小，当y0＞＞H时，采样点的横向间隔可化简为：

[0074] 3c)计算采样点的纵向间隔：

[0075] 3c1)计算点(x0,y0,z0)相对于雷达的径向速度：vr＝vcosα0cosβ0，

[0076] 其中v为散射单元的运动速度，α0为方位角，

[0077] β0为俯仰角，

[0078] 3c2)计算点(x0,y0,z0)的多普勒频率： 其中fc为信号载频；

[0079] 3c3)计算点(x0,y0,z0)的方位角分辨率：

[0080] 将3c2)中的多普勒频率fD对α0微分，即 得到方位角分辨率为： 其中d(·)为微元；

[0081] 3c4)计算方位角分割间隔：

[0082] 雷达信号的多普勒分辨率为 其中N为一个相参处理间隔CPI内的脉冲数，Tr为脉冲重复周期，为了保持多普勒分辨率，应该有d(fD)≤ΔfD，取d(fD)＝ΔfD、[0083] v取使多普勒不模糊范围内的最大值，若不模糊多普勒频率范围为 则v
的最大取值满足： 将上述d(fD)、α0、v参数代入到3c3)中的方位角分辨率d(α0)中，计算出方位角分割间隔为：

[0084] 为达到不模糊多普勒频率 设N的取值范围为：N≥2；为使δ的划分尽可能精细，可以进行r次取点，r＝1,2,...，得到{β1,β2,...,βr}，取β0＝min{β1,β2,...,βr}；

[0085] 3c5)根据3c4)中的方位角分割间隔δ与波门起始距离在y轴上投影的坐标s，计算海表面采样点的纵向间隔为sδ；

[0086] 3d)在矩形框内的海表面上取采样点：

[0087] 参照图5，以横向间隔 和纵向间隔sδ在矩形框内的海表面上取采样点，采样点在基准水平面上的投影为图5(a)所示的黑色圆点，再以每3个采样点构成一个三角形散射单元，将矩形框内的海表面划分为由许多三角形散射单元拼接而成的复杂表面，散射单元在基准水平面上的投影为图5(b)所示的直角三角形。

[0088] 步骤4，计算任意一个三角形散射单元的参数和回波功率。

[0089] 参照图6任意一个三角形散射单元与雷达位置关系，该三角形散射单元在水平面的投影为图5(b)中的某个直角三角形。该三角形散射单元的参数包括：三角形散射单元与雷达之间的距离R，方位角α与俯仰角β，来波方向 法向 面积A0。

[0090] 4a)计算三角形散射单元的参数：

[0091] 将任意一个三角形散射单元的三个顶点坐标表示为(x1,y1,z1)，(x2,y2,z2)，(x3,y3,z3)，其中心位置坐标为：

[0092]

[0093] 4a1)计算该三角形散射单元与雷达之间的距离：

[0094]

[0095] 4a2)计算该三角形散射单元在雷达天线坐标系下的方向(α,β)：

[0096] 方位角α：

[0097] 俯仰角β：

[0098] 4a3)计算该三角形散射单元的来波方向：

[0099] 4a4)计算该三角形散射单元的法向：

[0100] 4a5)计算该三角形散射单元的面积：

[0101] 4b)计算该三角形散射单元的回波功率：

[0102] 4b1)计算该三角形散射单元局部入射余角： 它的有效取值范围为

[0103] 4b2)根据Guinard和Dalcy给出的计算公式，计算在水平极化和垂直极化下的海面散射系数ΓHH与ΓVV：

[0104]

[0105]

[0106] 其中ε为复介电常数：ε＝ε′-jε″，其中 为虚数，ε′表示ε的一阶导，ε″表示ε的二阶导，ε′与ε″的值可通过查表得到，

[0107] 表1给出了一些在不同频率和温度时海水复介电常数的ε′与ε″值：

[0108] 表1海水的复介电常数

[0109]

[0110] 4b3)根据天线的极化方式选择对应的公式计算该三角形散射单元的雷达散射截面积：

[0111] 天线水平极化时，该三角形散射单元的雷达散射截面积：σ＝ΓHHΑ0，

[0112] 天线垂直极化时，该三角形散射单元的雷达散射截面积：σ＝ΓVVΑ0；

[0113] 4b4)根据上述参数计算该三角形散射单元回波功率：

[0114] 其中Pt为雷达发射功率，G(αF,βF；α,β)为波束瞄准线指向(αF,βF)时(α,β)方向上的天线增益，λ为信号波长，Ls为雷达收发综合损耗。

[0115] 步骤5，建立三角形散射单元的回波模型。

[0116] 5a)计算雷达发射波：

[0117] 设雷达发射波为线性调频脉冲信号：

[0118] st(t)＝rect(t/T)exp(j2πfct)exp(jπΒt2)，

[0119] 其中T为脉冲宽度，B为信号调频率，

[0120] 5b)计算三角形散射单元的点散射回波信号：

[0121] 其中fd为三角形散射单元的多普勒频率，可由三角形散射单元在两个脉冲间相对于雷达的距离变化估算出， Δt为两个脉冲的时间差；

[0122] 5c)建立三角形散射单元的回波模型：

[0123] 三角形散射单元属于海面散射单元，不会是平整的三角形散射面，因此要考虑海面的纹理分量对雷达回波的贡献，给3b)中的点散射回波模型相位延迟加上抖动模拟纹理分量的作用，得到三角形散射单元的回波模型：

[0124]

[0125] 其中ζ(t)为相位延迟抖动，是均值为零、方差为σ0的高斯白噪声，代表了抖动的剧烈程度，其与海况有关。

[0126] 步骤6，判断三角形散射单元之间的遮挡关系和三角形散射单元自身的遮挡关系。

[0127] 实际中，海面三角形散射单元之间会存在遮挡，三角形散射单元本身也会存在自身遮挡，采用射线跟踪和如下原则和判据进行判断：

[0128] 设某三角形散射单元顶点相对于雷达的仰角为α1，遮挡单元顶点相对于雷达的仰角为α2，当α1≥α2时，该三角形散射单元被遮挡，不计算其回波；当α1＜α2时，该三角形散射单元没被遮挡单元遮挡，则计算该三角形散射单元的局部入射角αL，再判断该三角形散射单元自身的遮挡关系：

[0129] 当αL＞π/2时，该三角形散射单元被自身遮挡，不计算其回波；

[0130] 当αL≤π/2时，三角形散射单元未被遮挡，计算并叠加所有未被遮挡的三角形散射单元的回波，产生一个相参处理间隔CPI回波的海杂波仿真数据，完成对海杂波的模拟。

[0131] 步骤7，分析海杂波的仿真数据。

[0132] 7a)将步骤6中产生的一个相参处理间隔CPI回波的海杂波仿真数据存储于大小为u×v的矩阵Q中，其中u为一个CPI包含的脉冲数，v为距离单元的个数；

[0133] 7b)改变天线高度参数，产生两组海杂波仿真数据，分别为Q1、Q2，每组均取第u1个脉冲作为单次脉冲回波的海杂波仿真数据，分别为 其中u1表示脉冲数；分别画出和 的仿真幅度-距离单元曲线图，比较两曲线的幅度；

[0134] 7c)根据线性调频脉冲压缩雷达采集的单次脉冲回波的海杂波实测数据Qre画出实测幅度-距离单元曲线图，将 Qre三者的曲线图作比较，分析仿真幅度-距离单元曲线图与实测幅度-距离单元曲线图的幅度趋势分布；

[0135] 7d)产生多组海杂波仿真数据，任选一组Q3，画出二维回波图，选择其中的第v1个距离单元的海杂波仿真数据 画出它的幅度-脉冲曲线图与幅度统计分布曲线图，将其幅度统计分布曲线与传统统计分布模型曲线进行最小均方误差MSD拟合优度检验，其中v1表示距离单元数。

[0136] 下面结合具体的仿真测试结果进一步说明本发明的有益效果。

[0137] 1.仿真条件

[0138] 如表2所列：

[0139] 表2仿真参数

[0140]

[0141]

[0142] 仿真1：在表2条件下，仿真慢时间的海表面海浪，结果如图7。图7显示了4个不同时刻的100m×100m的仿真海表面二维图，表明本发明建立的海表面海浪模型，能仿真出海表面海浪的动态变化。

[0143] 仿真2：在表2条件下加入天线高度分别为H＝30m与H＝500m的条件，在海表面上仿真海杂波，分析比较不同天线高度的海杂波仿真数据以及海杂波的仿真数据与实测数据，结果如图8。其中：

[0144] 图8(a)为天线高度为30m时的单次脉冲回波的海杂波的仿真幅度-距离单元曲线图；

[0145] 图8(b)为天线高度为500m时的单次脉冲回波的海杂波仿真幅度-距离单元曲线图；

[0146] 图8(c)为单次脉冲回波的海杂波实测幅度-距离单元曲线图，其中海杂波的实测数据由线性调频LFM脉冲压缩雷达采集得到，该雷达位于某山顶照射太平洋，天线高度为500m，入射余角约为1°，带宽5MHz，脉宽42μs；

[0147] 比较图8(a)与图8(b)可知，图8(b)回波幅度较图8(a)回波幅度稍大，表明本发明产生的海杂波符合海杂波后向散射系数随波束瞄准线入射余角的增加而增加的结论；

[0148] 图8(a)与8(b)中的仿真幅度-距离单元曲线图从整体看都呈现出中间高两边低的趋势，与8(c)中的实测幅度-距离单元曲线图趋势一致，这是因为波门开得较小，处于波门两边的散射单元回波未完全包含进波门。

[0149] 图8表明本发明产生的海杂波具有可信度。

[0150] 仿真3：在表2条件下再加入固定天线高度为500m，波束瞄准线入射余角为1°的条件，仿真海杂波，对海杂波的幅度统计分布与传统统计分布模型进行最小均方误差MSD拟合优度检验，结果如图9，其中：

[0151] 图9(a)为一个相参处理间隔的海杂波仿真回波二维图；

[0152] 图9(b)为第500个距离单元的海杂波回波仿真幅度-脉冲曲线图；

[0153] 图9(c)为仿真幅度-脉冲曲线图的幅度统计分布与传统统计分布模型的拟合结果；

[0154] 图9(d)为第500个距离单元回波仿真数据的多普勒谱；

[0155] 图9(c)中对仿真采用最小均方误差MSD拟合优度检验得到的结果如表3所示，由表3可见K分布拟合结果最好，瑞利分布拟合结果次之；

[0156] 表3统计分布拟合结果

[0157]

[0158] 图9(d)中的第500个距离单元仿真数据的多普勒谱有一定的频移和谱宽，因为仿真海面上三角形散射单元相对于雷达的方位角不一样，且与雷达之间的距离和距离变化不同；多普勒谱可看作被天线主波束照射海面的三角形散射单元的速度映射，图9(c)和图9(d)与传统统计模型结论一致。

[0159] 图9结果表明：本发明的海杂波模拟方法具有实用性，可以用来研究统计模型参数、雷达工作参数与环境参数之间的关系。