CN115567353B - 一种面向雷达组网系统的干扰多波束调度和干扰功率联合优化方法 - Google Patents

Abstract

本发明属于干扰资源管理领域，特别针对雷达组网系统，提出了一种面向雷达组网系统的干扰多波束调度和干扰功率联合优化方法。本发明在离散粒子群算法(DPSO)的框架下，联合优化干扰机的波束调度和干扰功率，其中干扰效能综合考虑了组网雷达系统在搜索阶段对目标的检测概率以及跟踪阶段对目标的瞄准概率，二者加权后将其作为干扰性能评估指标。本发明能够根据干扰场景的变化，自适应配置干扰机波束方向和功率资源，在有限的干扰资源约束下实现更有效的干扰。

Description

一种面向雷达组网系统的干扰多波束调度和干扰功率联合优 化方法

技术领域

本发明属于干扰资源管理领域，提出了一种面向雷达组网系统的干扰多波束调度和干扰功率联合优化的压制干扰策略

背景技术

随着信息融合、数据处理和数字信号处理技术的高速发展，现代雷达的性能不断提高，在目标探测，跟踪和制导方面都有了长足的进步，这对于电子对抗系统提出了新的挑战。雷达干扰作为现代电子对抗系统中不可或缺的一种重要方式，为了削弱敌方雷达系统的作战性能，涌现了各种干扰敌方雷达的技术手段，包括不同类型的压制性干扰和不同方式的欺骗干扰。在压制干扰作战场景当中，干扰系统的资源往往是有限的，采用固定资源分配方式对敌方雷达的干扰效果往往不理想，资源未得到合理的充分利用，难以最大程度限制敌方雷达系统带来的威胁。因此如何合理配置有限的干扰资源，获得最大的整体干扰效果成为一个重要的研究课题。

目前对干扰资源分配问题的研究已经取得了一些进展。在文献(Shen Y,Chen YG,Li XH.Research on optimal distribution of radar jamming resource based onZero-one programming.Acta Armamentarii 2007；28(5):528-532.)当中，针对地对空干扰设备对抗来袭目标场景，研究了敌方目标威胁等级判定以及干扰效果评估方法，采用了0-1规划，运用匈牙利解法得到对应的匹配关系，但是其中的干扰机为单波束模式，未涉及到干扰功率的分配问题。在文献(Liu X,Li DS,Hu R.Application of improved geneticalgorithm in cooperative jamming resource assignment.Journal of Detection&Control 2018；40(05):69-75.)中，针对敌方雷达组网的多干扰机协同干扰资源分配问题，以压制概率为目标函数，使用改进的遗传算法解决了不同干扰样式下的雷达与干扰机匹配问题，实现了匹配关系和干扰方式的联合自适应，但遗传算法对初始种群的选择具有一定的依赖性，在搜索过程中容易陷入“过早”，因此遗传算法在探索新解空间的能力上受到限制。此外，遗传算法没有及时利用反馈信息，导致搜索速度较慢。文献(Wang SJ,LiHD.Research on improved gray wolf algorithm in interference resourceallocation.Applied Science and Technology 2021；48(02):54-57+99.)针对多干扰机干扰雷达组网问题，选择雷达定位精度(GDOP)作为评价指标，对协同干扰资源分配问题采用改进的灰狼算法，在迭代速度上有一定的提升，但仍有陷入局部最优解的可能，算法的稳定性并不高。在文献(Han P,Lu JD,Wang XL.Radar active jamming resourceassignment algorithm based on game theory.Modern Defense Technology 2018；46(04):53-59.)当中，利用博弈论分析干扰资源分配问题，证明了纳什均衡的可行性和存在性，设计迭代干扰策略选择算法，求解出最佳的干扰匹配关系，但算法复杂度太高，降低了收敛速度。文献(Qin Q,Dong W,Lin M,et al.Cooperative jamming resourceallocation of uav swarm based on multi-objective DPSO[C]//2018Chinese ControlAnd Decision Conference(CCDC).IEEE,2018.)中，针对无人机群协同干扰资源分配问题，建立了一种包含干扰效益、损失代价和飞行距离代价的多目标函数优化模型，采用离散粒子群优化算法(DPSO),求解无人机与雷达的匹配关系。文献(尧泽昆,王超,施庆展,张少卿,袁乃昌.基于改进离散模拟退火遗传算法的雷达网协同干扰资源分配模型[J/OL].系统工程与电子技术，2022，1-8)针对分布式干扰机掩护目标突防雷达网背景下的干扰资源分配问题，基于雷达组网融合概率构建了干扰效果评估函数，利用IDSA-GA方法对函数进行寻优求解，联合优化了干扰样式和干扰匹配关系。

上述文章的资源管理策略当中，大都是通过优化波束调度，即干扰机与雷达的匹配关系，从而实现最优干扰。每台干扰机只能发射单个波束，采用固定功率进行干扰，随着多波束干扰机的发展，一台干扰设备可以利用发射的多个干扰波束同时干扰多目标。在文献(GAO Xiaoguang,HU Ming,and ZHENG Jingsong.Jamming strategy for single planeto multi-target in task of penetration[J].Systems Engineering andElectronics,2010,32(6):1239–1243.doi:10.3969/j.issn.1001-50506X.2010.06.028.)当中，针对单机多目标的突防任务，通过贝叶斯理论来建立较为完善的威胁等级评估模型，提出了基于威胁等级评估的干扰功率分配方法。在文献(ZHANG Dalin,YI Wei,and KONGLingjiang.Optimal joint allocation of multi-jammer resources for jammingnetted radar system[J].Journal of Radars,2021,10(4):595–606.doi:10.12000/JR21071.)当中，针对多目标突防组网雷达的场景，建立了包含干扰波束方向和发射功率两个优化变量的资源优化模型，并利用粒子群算法对问题进行了求解，但对于两个优化变量，采用序贯优化方式，即利用智能算法依次获得最佳干扰匹配关系和干扰功率分配结果，增大了计算量。文献(陆德江,王星,陈游,胡星.联合多种资源协同干扰组网雷达系统的自适应调度方法[J/OL].系统工程与电子技术，2022，1-15)当中，针对飞机编队协同干扰组网雷达系统的资源分配问题，提出了一种联合目标选择与功率分配的自适应调度方法，基于此给出了一种结合改进布谷鸟搜索算法与KKT优化条件的求解方法。

上述干扰资源分配策略当中，关于“多对多”的干扰资源调度研究相对较少，相关研究更多的把关注度放在如何确定干扰机和组网雷达各节点的干扰对应关系上，可管控的干扰资源类型较为单一，当添加新的可控参数时，又是在已确定匹配关系的条件下重新采用智能算法等进行求解，大大增加了计算复杂度。同时用于评估干扰效果的目标函数，大多采用各台雷达对目标的检测概率或压制概率，并未考虑干扰对雷达跟踪阶段性能的影响,不足以评价真实的干扰效果。

针对上述问题，本发明提出了一种面向雷达组网系统的干扰多波束调度和干扰功率联合优化方法，建立了以检测概率和瞄准概率为干扰效果评价指标的双因子评估函数。采用DPSO算法，同时实现对干扰机与雷达的匹配关系和干扰机辐射功率的优化。

发明内容

假设在t时刻，有M架干扰机对N个雷达节点进行干扰，每台干扰机的总辐射功率上限为P_total，每台干扰机最多发射B个波束，每个雷达节点最多可受到S个干扰机的干扰，一种面向雷达组网系统的干扰多波束调度和干扰功率联合优化方法具体步骤如下：

步骤1：随机初始化共计D个粒子，每个粒子为一个N行M列的矩阵，令k＝1，第d个粒子在第k次迭代时表示为其中第n行m列元素P_n,m表示干扰机m针对雷达节点n的干扰功率。矩阵/>中元素在0-P_total间随机选取，并按照下式进行归一化处理：

步骤2:按照各粒子的对应干扰关系处的干扰功率，计算每个粒子的适应度值如下:

其中，Q表示目标的数目，w_q表示目标q的重要程度，ρ₁，ρ₂为权重，满足ρ₁+ρ₂＝1。/>分别表示t时刻组网系统针对目标q的检测概率和瞄准概率。

(2)式当中组网系统针对目标q的检测概率计算如下：

其中，N为雷达节点数目，H表示当组网内有H部或多于H部的雷达同时发现目标时，组网雷达系统判定为发现目标，d_n为每个节点针对自身获得的信息，做出一个局部判断结果d_n(d_n∈{0,1})，表示了所有雷达节点对于目标q的局部判断结果之和为j的所有可能排列，/>为t时刻雷达节点n对目标q的检测概率，

(3)式中雷达节点n对目标q的检测概率计算如下：

其中，P_fa表示虚警概率，I表示非相参积累脉冲数目，表示雷达节点n接收到的信干比，计算如下：

其中，

P_n＝k₀T_tΔf_rF_n (6)

表示接收机内部噪声，其中k₀为玻尔兹曼常数，T_t表示接收机噪声温度，Δf_r为接收机带宽，F_n为接收机噪声系数。

表示雷达节点n接收到的自身回波功率，其中P_r,n表示雷达节点n的发射功率，G_r,n表示雷达天线增益，σ_q表示目标q的雷达截面积，λ_r表示雷达波长，L_r,n表示损耗，表示t时刻，雷达节点n和目标q之间的距离，其中(x_n,y_n)，(x_q,y_q)分别表示t时刻雷达节点n和目标q的位置。

表示雷达节点n接收到的干扰机m的干扰功率，表示干扰机m针对雷达节点n的干扰机发射功率，G_j.m表示干扰机m发射增益，λ_j表示干扰机发射波长，L_j,m表示损耗，表示t时刻雷达节点n与干扰机m之间的距离，其中(x_m,y_m)表示t时刻干扰机m的位置。

(2)式中组网系统针对目标q的瞄准概率计算如下：

其中，表示雷达节点n针对目标q的瞄准概率，分别表示距离，速度，方位角三个指标的瞄准概率，计算如下：

Δr,Δv,Δθ分别表示上述三个参数的测量误差允许范围。分别表示三个参数的测量误差的概率密度函数形式，

其中，分别表示三个参数的测量误差均方误差，计算如下：

其中c为光速，τ为脉冲宽度，为信干比，如(5)式所示，θ_0.5为雷达天线主瓣宽度，λ_r为雷达工作波长。

初始化每个粒子的自身最优位置为其本身每个粒子的位置/>选择适应度值最小的一个粒子，将其设定为全局最优位置x_gbest。

步骤3：根据粒子的维度，获得粒子运行速度矢量为N行M列的矩阵v。矩阵v中元素为[-v,v]间均匀分布。对每个粒子的速度和位置进行更新，根据自身的历史最优位置和全局最优位置，更新每个粒子的速度和位置。其中，第d个粒子的第k次速度更新公式如下：

其中，w_k为惯性因子，计算如下:

其中，K_max表示粒子群最大迭代次数，与/>为自身认知因子，计算如下：

位置更新公式为：

经过这一步的粒子更新后，矩阵中元素可能超出功率范围[0,P_total]，当某个粒子中元素小于0时，将其改为0。之后，按照(1)式将各干扰机的功率进行归一化，保证各台干扰机的干扰波束功率之和为P_total。

步骤4：评估更新后的粒子的适应度值。经过步骤3，每个粒子的位置得到了更新，由新的可行解按照(2)式计算其适应度值，并与其自身最优值/>的适应度值进行比较，如果当前值比/>更优，则将/>设置为当前值位置/>否则，保持原有的最优位置。将所有粒子进行上述操作，更新每个粒子的自身最优位置/>选择更新后的/>中适应度值最好的，并与x_gbest进行比较，如果优于原来的，则将其设定为当前时刻的全局最优位置x_gbest，否则保持原来的x_gbest。

步骤5：重复步骤3到步骤4，直到k达到了最大迭代次数K_max，输出x_gbest。之后，由x_gbest计算当前时刻的匹配关系。

步骤6：初始化可行解矩阵U_t为N行M列的0矩阵,令x_con＝x_gbest，遍历x_con中的各个元素，找到其最大值所在位置[n_l,m_l]，并用0替换，之后，令U_t中对应位置为1，并判断是否满足U_t每列之和小于等于B,每行之和小于等于S的约束，如果不满足，则该位置保持为0，不再更新为1，重复上述操作，直到遍历完矩阵x_con中的N×M个元素，得到由松弛解x_con对应的可行解矩阵U_t。

步骤7：令P_t＝U_t⊙x_gbest，⊙表示哈达玛积，得到由匹配关系U_t对应的功率矩阵P_t，此时，重新按照(1)式更新P_t。

采用匹配关系U_t,以及更新后的P_t中对应的各波束的功率情况，计算目标函数，作为t时刻的优化结果。

发明原理

本发明所研究的干扰资源分配问题，是在满足一定的系统要求和有限的干扰资源约束下，合理的分配干扰机的波束方向和发射功率。首先，定义一个用于表示t时刻干扰机与雷达节点干扰关系的匹配关系矩阵：

式中

另外，与匹配关系矩阵所对应，定义一个用于表示t时刻干扰机各波束发射功率的功率分配矩阵：

式中，表示了在t时刻，干扰机m针对雷达节点n的干扰功率，满足：

此外，每个干扰机的总功率是一定的，满足

同时，为了符合现实情况，规定每架干扰机在每个时刻可同时干扰的雷达节点数目是有限的，表示为：

式中，B表示每个时刻每架干扰机可同时干扰的最大雷达节点数目。

出于计算复杂度和干扰资源合理分配的考虑，在每个时刻每个雷达节点分配的干扰机数目也是有限的，应满足：

式中，S表示每个时刻每个雷达节点可分配的最大干扰机数目。

检测概率是表征雷达在搜索阶段能否发现目标的重要性能指标，表示目标真实存在时，雷达判定目标存在的概率。假定在某一时段内，雷达节点n接收到一段脉冲序列，包含两种情况：(1)包含有目标q的回波信号以及接收机内部噪声，(2)包含有目标q的回波信号，接收机内部噪声，以及干扰机对该雷达节点的干扰信号。

由雷达方程易得雷达节点n接收到的目标q的回波信号如(7)式所示，接收机内部噪声一般由(6)式表示，而干扰信号由(8)式表示，由(18)式来表示雷达节点是否受到干扰。

由此，雷达节点n接收到的脉冲信号可表示为：

结合(23)式，可推导在t时刻，雷达节点n接收到的关于目标q的信干比为：

假定雷达采用非相参脉冲累积技术，根据雷达检测原理，可计算得到在t时刻，雷达节点n对目标q的检测概率公式为(4)式所示，接下来，讨论组网雷达系统针对目标q的检测概率，此处假定组网系统是根据秩H准则进行的判定，即组网系统中每个节点根据自身的回波处理结果做出局部判断d_n(d_n∈{0,1})，并将该局部结果送往信息融合中心，根据这些局部结果生成全局判断矢量D＝f_D(d₁,d₂,···,d_N)。D共有2^N种可能，即

之后，组网系统在此基础上进行目标判定，判定的规则为：

因此，由秩H判断准则可得到在t时刻，组网系统针对目标q的检测概率模型如(3)式所示。

当组网系统处于跟踪状态时，为了反映雷达在该阶段对目标的跟踪精度，引入瞄准概率这一指标进行度量。此处考虑了距离，速度和方位角三个参数，并假定其相互独立，近似服从0均值的高斯分布，可得上述三个参数的测量误差均方误差由(12)式表示，接下来依据高斯分布的概率模型，可得各参数的测量误差的概率密度函数为(11)式所示，给出三个参数的允许测量误差的最大上限[Δr,Δv,Δθ]，对(11)式进行积分可得相应参数的瞄准概率如(10)式所示。考虑到三个参数中任何一个参数的测量误差超出允许范围都会导致雷达对目标的跟踪失败，因此，t时刻雷达节点n对目标q的瞄准概率模型为与检测概率类似，采用秩H判断准则，可得组网系统对目标q的瞄准概率如(9)式所示。之后对检测概率和瞄准概率进行加权可得算法的判定指标如(2)式所示。

以检测概率和瞄准概率的加权和作为干扰效果的评价指标，可得优化模型如下：

根据最后两个约束条件，可以看到，若的最佳值为0，意味着/>为0，/>非零即1的取值，可以从最佳的/>取值得到。因此，考虑采用单纯优化P^t的算法去同时获得P^t和U^t的最佳值，上述优化问题可转变为：

对于上述优化问题，可采用离散粒子群算法(DPSO)进行求解。首先假定每台干扰机对各个雷达节点均可进行干扰，以此构建D个如(19)式的初始化粒子，为了满足(28)式中约束条件1，需按照(1)式进行归一化处理，并初始化自身最优值和全局最优值x_gbest，如步骤1、2所示。随后对粒子进行更新，包括速度和位置的更新，如步骤3所示，为了满足(28)式中约束条件1，步骤1最后亦要进行归一化处理，基于此更新自身最优解和全局最优解，如步骤4所示。通过迭代处理，得到初步功率优化结果x_gbest，如步骤5所示。

上述DPSO方式得到的功率分配矩阵x_gbest中各个位置上元素的大小实质反映了该元素索引位置上的波束指向优化结果为1的可能性大小，因此，可按照其中元素从大到小的顺序去确定相应的匹配关系，其中为了满足(27)式中的约束条件1，2，在优化U^t的过程中，需要判断是否每列之和小于等于B，每行之和小于等于S，如步骤6所示。

最后，根据最优的匹配关系矩阵和离散粒子群算法获得的解，得到实际的最优发射干扰功率矩阵，如步骤7所示。

附图说明

图1为仿真场景图

图2为目标函数对比图

图3为雷达节点1干扰模式图

图4为雷达节点2干扰模式图

图5为雷达节点3干扰模式图

图6为雷达节点4干扰模式图

图7为雷达节点5干扰模式图

图8为干扰机1功率分配图

图9为干扰机2功率分配图

图10为干扰机3功率分配图

具体实施方式

假设在一个二维干扰场景中，存在一个分布式雷达组网系统，其中包含5部位置不同的雷达，雷达位置坐标分别为(48km，5km)、(100km，5km)、(160km，5km)、(70km，80km)、(130km，80km)。有3台干扰机经过这片区域，3台干扰机的初始位置坐标分别为(2km，105km)、(1.6km，95.5km)、(1.6km，85.5km)仿真场景如图1所示。仿真参数设置为：雷达发射功率P_r,j＝9×10⁹W，虚警概率P_f＝10^-6，雷达天线主瓣宽度θ_0.5＝2,干扰机发射总功率P_total＝120W,雷达发射波长λ_r＝0.1m,雷达截面积RCS＝1m²,雷达发射增益G_r,i＝40dB,干扰机发射波长λ_j＝0.1m,接收机带宽Δf＝10⁷Hz,干扰机发射增益G_j,m＝10dB,干扰机损耗L_j,m＝6dB,非相干累积脉冲数目I＝20,脉冲宽度τ＝10^-5s,雷达损耗L_r,i＝6dB,粒子群最大迭代次数K_max＝50，粒子个数D＝200，粒子更新速度v＝5，每个雷达节点可分配的最大干扰机数目S＝2，每架干扰机可同时干扰的最大雷达节点数目B＝3。

如图2所示，将本发明算法与固定匹配关系-均分干扰功率算法、就近匹配关系-均分干扰功率算法、优化匹配关系-均分干扰功率算法进行对比。本算法具有最低的目标函数值。

如图3-7所示，本算法可以实现在不同的时刻，随着跟踪场景的变化，自适应的选择最合适的干扰匹配关系，从而实现最佳干扰。图中，干扰模式1表示只有干扰机1进行干扰，干扰模式2表示只有干扰机2进行干扰，干扰模式3表示只有干扰机3进行干扰，干扰模式4表示干扰机1和干扰机2同时对其进行干扰，干扰模式5表示干扰机1和干扰机3同时对其进行干扰，干扰模式6表示干扰机2和干扰机3同时对其进行干扰。

如图8-10所示，各台干扰机可以自适应的实现干扰功率的合理分配，干扰机1在整个过程当中，将大部分功率资源用于干扰雷达节点4和5，且初始阶段大部分用于干扰雷达节点4，在干扰后期大部分用于干扰节点5。干扰机3将大部分功率资源用于干扰雷达节点1，2，3，且干扰中间阶段大部分功率资源用于了雷达节点2。而干扰机2，在初始阶段将大部分功率资源分配给节点2，3，4，在后期，则将大部分功率资源分配给节点1，2，4。

综上，本发明所提算法一种面向雷达组网系统的干扰多波束调度和干扰功率联合优化方法，能够根据实际场景，实时调整干扰机与雷达之间的匹配关系，以及对干扰机各波束功率的自适应分配，在有限的干扰资源约束条件下，实现最佳干扰。

Claims (1)

1.一种面向雷达组网系统的干扰多波束调度和干扰功率联合优化方法，具体技术方案如下：

假设在t时刻，有M架干扰机对N个雷达节点进行干扰，每台干扰机的总辐射功率上限为P_total，每台干扰机最多发射B个波束，每个雷达节点最多可受到S个干扰机的干扰，一种面向雷达组网系统的干扰多波束调度和干扰功率联合优化方法具体步骤如下：

步骤1：随机初始化共计D个粒子，每个粒子为一个N行M列的矩阵，令k＝1，第d个粒子在第k次迭代时表示为其中第n行m列元素P_n,m表示干扰机m针对雷达节点n的干扰功率。矩阵/>中元素在0-P_total间随机选取，并按照下式进行归一化处理：

步骤2:按照各粒子的对应干扰关系处的干扰功率，计算每个粒子的适应度值如下:

其中，Q表示目标的数目，w_q表示目标q的重要程度，ρ₁，ρ₂为权重，满足ρ₁+ρ₂＝1；/>分别表示t时刻组网系统针对目标q的检测概率和瞄准概率；

(2)式当中组网系统针对目标q的检测概率计算如下：

其中，N为雷达节点数目，H表示当组网内有H部或多于H部的雷达同时发现目标时，组网雷达系统判定为发现目标，d_n为每个节点针对自身获得的信息，做出一个局部判断结果表示了所有雷达节点对于目标q的局部判断结果之和为j的所有可能排列，/>为t时刻雷达节点n对目标q的检测概率，

(3)式中雷达节点n对目标q的检测概率计算如下：

其中，P_fa表示虚警概率，I表示非相参积累脉冲数目，/>表示雷达节点n接收到的信干比，计算如下：

其中，

P_n＝k₀T_tΔf_rF_n (6)

表示接收机内部噪声，其中k₀为玻尔兹曼常数，T_t表示接收机噪声温度，Δf_r为接收机带宽，F_n为接收机噪声系数；

表示雷达节点n接收到的自身回波功率，其中P_r,n表示雷达节点n的发射功率，G_r,n表示雷达天线增益，σ_q表示目标q的雷达截面积，λ_r表示雷达波长，L_r,n表示损耗，表示t时刻，雷达节点n和目标q之间的距离，其中(x_n,y_n)，(x_q,y_q)分别表示t时刻雷达节点n和目标q的位置；

表示雷达节点n接收到的干扰机m的干扰功率，表示干扰机m针对雷达节点n的干扰机发射功率，G_j.m表示干扰机m发射增益，λ_j表示干扰机发射波长，L_j,m表示损耗，表示t时刻雷达节点n与干扰机m之间的距离，其中(x_m,y_m)表示t时刻干扰机m的位置；

(2)式中组网系统针对目标q的瞄准概率计算如下：

其中，表示雷达节点n针对目标q的瞄准概率，分别表示距离，速度，方位角三个指标的瞄准概率，计算如下：

Δr,Δv,Δθ分别表示上述三个参数的测量误差允许范围；分别表示三个参数的测量误差的概率密度函数形式，

其中，分别表示三个参数的测量误差均方误差，计算如下：

其中c为光速，τ为脉冲宽度，为信干比，如(5)式所示，θ_0.5为雷达天线主瓣宽度，λ_r为雷达工作波长；

初始化每个粒子的自身最优位置为其本身每个粒子的位置/>选择适应度值最小的一个粒子，将其设定为全局最优位置x_gbest；

步骤3：根据粒子的维度，获得粒子运行速度矢量为N行M列的矩阵v。矩阵v中元素为[-v,v]间均匀分布。对每个粒子的速度和位置进行更新，根据自身的历史最优位置和全局最优位置，更新每个粒子的速度和位置。其中，第d个粒子的第k次速度更新公式如下：

其中，w_k为惯性因子，计算如下:

其中，K_max表示粒子群最大迭代次数，与/>为自身认知因子，计算如下：

位置更新公式为：

经过这一步的粒子更新后，矩阵中元素可能超出功率范围[0,P_total]，当某个粒子中元素小于0时，将其改为0；之后，按照(1)式将各干扰机的功率进行归一化，保证各台干扰机的干扰波束功率之和为P_total；

步骤4：评估更新后的粒子的适应度值；经过步骤3，每个粒子的位置得到了更新，由新的可行解按照(2)式计算其适应度值，并与其自身最优值/>的适应度值进行比较，如果当前值比/>更优，则将/>设置为当前值位置/>否则，保持原有的最优位置。将所有粒子进行上述操作，更新每个粒子的自身最优位置/>选择更新后的/>中适应度值最好的，并与x_gbest进行比较，如果优于原来的，则将其设定为当前时刻的全局最优位置x_gbest，否则保持原来的x_gbest；

步骤5：重复步骤3到步骤4，直到k达到了最大迭代次数K_max，输出x_gbest。之后，由x_gbest计算当前时刻的匹配关系；

步骤6：初始化可行解矩阵U_t为N行M列的0矩阵,令x_con＝x_gbest，遍历x_con中的各个元素，找到其最大值所在位置[n_l,m_l]，并用0替换，之后，令U_t中对应位置为1，并判断是否满足U_t每列之和小于等于B,每行之和小于等于S的约束，如果不满足，则该位置保持为0，不再更新为1，重复上述操作，直到遍历完矩阵x_con中的N×M个元素，得到由松弛解x_con对应的可行解矩阵U_t；

步骤7：令P_t＝U_t⊙x_gbest，⊙表示哈达玛积，得到由匹配关系U_t对应的功率矩阵P_t，此时，重新按照(1)式更新P_t；

采用匹配关系U_t,以及更新后的P_t中对应的各波束的功率情况，计算目标函数，作为t时刻的优化结果。