✅作者简介热爱科研的Matlab仿真开发者擅长毕业设计辅导、数学建模、数据处理、建模仿真、程序设计、完整代码获取、论文复现及科研仿真。 往期回顾关注个人主页Matlab科研工作室 关注我领取海量matlab电子书和数学建模资料个人信条格物致知,完整Matlab代码获取及仿真咨询内容私信。 内容介绍一、引言在信号处理领域去除噪声对于获取准确有效的信息至关重要。完备总体经验模态分解自适应噪声ICEEMDAN是一种强大的信号分解方法然而其分解过程可能存在模态混叠等问题。开普勒优化算法KOA是一种新兴的元启发式优化算法受开普勒定律启发具有良好的全局搜索能力。将 KOA 应用于 ICEEMDAN可优化其分解参数实现更高效的信号去噪。二、相关理论基础一完备总体经验模态分解自适应噪声ICEEMDAN基本原理ICEEMDAN 是对总体经验模态分解EEMD的改进。EEMD 通过多次添加不同的白噪声并进行经验模态分解EMD然后对结果进行平均来抑制模态混叠。ICEEMDAN 在此基础上通过自适应地添加特定的噪声进一步提高了分解的准确性和稳定性。它将原始信号x(t)逐步分解为一系列固有模态函数IMFci(t)和一个残余项r(t)即x(t)∑i1nci(t)r(t)。每个 IMF 满足两个条件一是在整个数据长度内极值点的数量与过零点的数量必须相等或最多相差一个二是在任意点由局部极大值点形成的上包络线和由局部极小值点形成的下包络线的均值为零。优势与不足ICEEMDAN 的优势在于能够自适应地分解非平稳、非线性信号对复杂信号具有良好的分解能力。但在实际应用中仍可能出现模态混叠现象即一个 IMF 中包含了不同尺度的信号成分或者不同 IMF 之间存在相似尺度的成分这会影响后续对信号的分析和处理。二开普勒优化算法KOA算法灵感与原理KOA 模拟了天体在太阳系中的运动依据开普勒定律进行搜索。在 KOA 中每个解被视为一个天体其位置和速度分别对应解的参数和搜索方向。通过更新天体的位置和速度算法在解空间中进行搜索。例如根据开普勒第二定律天体在相等时间内扫过相等的面积这在算法中转化为一种自适应的搜索策略使得算法在搜索初期能够快速探索解空间后期则更注重局部搜索以提高解的精度。算法流程初始化随机生成一组天体初始解并确定每个天体的初始位置和速度。更新位置和速度根据开普勒定律通过迭代更新天体的速度和位置。速度更新公式结合了引力和斥力的作用位置更新则基于速度的变化。适应度评估计算每个天体解的适应度值以评估解的优劣。适应度函数通常根据具体问题进行设计在信号去噪问题中可以是去噪后信号与原始纯净信号的误差度量。选择与更新根据适应度值选择较优的天体并更新种群。重复上述步骤直到满足终止条件如达到最大迭代次数或适应度值收敛。三、基于 KOA - ICEEMDAN 的信号去噪方法一参数优化思路ICEEMDAN 在分解信号时一些参数如添加噪声的强度、分解层数等对分解结果有重要影响。使用 KOA 对这些参数进行优化以找到最优的参数组合从而减少模态混叠提高去噪效果。具体来说将 ICEEMDAN 的参数空间作为 KOA 的搜索空间通过 KOA 的迭代搜索找到使去噪后信号质量最优的参数值。二KOA - ICEEMDAN 算法步骤初始化 KOA 参数设置 KOA 的种群规模、最大迭代次数、引力常数等参数。同时确定 ICEEMDAN 待优化的参数范围如噪声强度的取值范围、最大分解层数等。生成初始种群在 ICEEMDAN 参数空间内随机生成一组参数组合作为 KOA 的初始种群每个参数组合代表一个天体的初始位置。适应度计算对于种群中的每个参数组合将其应用于 ICEEMDAN 对含噪信号进行分解。然后计算去噪后信号的适应度值。适应度函数可以定义为去噪后信号与原始纯净信号的均方误差MSE的倒数即FitnessMSE1MSE 越小适应度值越大表示去噪效果越好。KOA 迭代优化按照 KOA 的算法流程更新天体的速度和位置即调整 ICEEMDAN 的参数组合。每次迭代后重新计算适应度值并根据适应度值选择较优的参数组合。确定最优参数与去噪当 KOA 满足终止条件时得到最优的 ICEEMDAN 参数组合。使用该参数组合对含噪信号进行 ICEEMDAN 分解得到一系列 IMF。通过对 IMF 进行筛选和重构去除噪声相关的 IMF重构得到去噪后的信号。⛳️ 运行结果 部分代码%开普勒优化算法function [Sun_Score,Sun_Pos,Convergence_curve]KOA(SearchAgents_no,Tmax,lb,ub,dim,feval)%%%%-------------------Definitions--------------------------%%%%Sun_Poszeros(1,dim); % A vector to include the best-so-far Solution, representing the SunSun_Scoreinf; % A Scalar variable to include the best-so-far scoreConvergence_curvezeros(1,Tmax);if length(lb)1lb lb*ones(1,dim);ub ub*ones(1,dim);end%%-------------------Controlling parameters--------------------------%%%%Tc3;M00.1;lambda15;%% Step 1: Initialization process%%---------------Initialization----------------------%%% Orbital Eccentricity (e)orbitalrand(1,SearchAgents_no); %% Eq.(4)%% Orbital Period (T)Tabs(randn(1,SearchAgents_no)); %% Eq.(5)Positionsinitialization(SearchAgents_no,dim,ub,lb); % Initialize the positions of planetst0; %% Function evaluation counter%%%%---------------------Evaluation-----------------------%%for i1:SearchAgents_no%% Test suites of CEC-2014, CEC-2017, CEC-2020, and CEC-2022PL_Fit(i)feval( Positions(i,:));% Update the best-so-far solutionif PL_Fit(i)Sun_Score % Change this to for maximization problemSun_ScorePL_Fit(i); % Update the best-so-far scoreSun_PosPositions(i,:); % Update te best-so-far solutionendendwhile tTmax %% Termination condition[Order] sort(PL_Fit); % Sorting the fitness values of the solutions in current population%% The worst Fitness value at function evaluation tworstFitness Order(SearchAgents_no); %% Eq.(11)MM0*(exp(-lambda*(t/Tmax))); %% Eq. (12)%% Computer R that represents the Euclidian distance between the best-so-far solution and the ith solutionfor i1:SearchAgents_noR(i)0;for j1:dimR(i)R(i)(Sun_Pos(j)-Positions(i,j))^2; %% Eq.(7)endR(i)sqrt(R(i));end%% The mass of the Sun and object i at time t is computed as follows:for i1:SearchAgents_nosum0;for k1:SearchAgents_nosumsum(PL_Fit(k)-worstFitness);endMS(i)rand*(Sun_Score-worstFitness)/(sum); %% Eq.(8)m(i)(PL_Fit(i)-worstFitness)/(sum); %% Eq.(9)end%% Step 2: Defining the gravitational force (F)% Computing the attraction force of the Sun and the ith planet according to the universal law of gravitation:for i1:SearchAgents_noRnorm(i)(R(i)-min(R))/(max(R)-min(R)); %% The normalized R (Eq.(24))MSnorm(i)(MS(i)-min(MS))/(max(MS)-min(MS)); %% The normalized MSMnorm(i)(m(i)-min(m))/(max(m)-min(m)); %% The normalized mFg(i)orbital(i)*M*((MSnorm(i)*Mnorm(i))/(Rnorm(i)*Rnorm(i)eps))(rand); %% Eq.(6)end%% a1 represents the semimajor axis of the elliptical orbit of object i at time t,for i1:SearchAgents_noa1(i)rand*(T(i)^2*(M*(MS(i)m(i))/(4*pi*pi)))^(1/3); %% Eq.(23)endfor i1:SearchAgents_no%% a2 is a cyclic controlling parameter that is decreasing gradually from -1 to ?2a2-1-1*(rem(t,Tmax/Tc)/(Tmax/Tc)); %% Eq.(29)%% ? is a linearly decreasing factor from 1 to ?2n(a2-1)*rand1; %% Eq.(28)arandi(SearchAgents_no); %% An index of a solution selected at randombrandi(SearchAgents_no); %% An index of a solution selected at randomrdrand(1,dim); %% A vector generated according to the normal distributionrrand; %% r1 is a random number in [0,1]%% A randomly-assigned binary vectorU1rdr; %% Eq.(21)O_PPositions(i,:); %% Storing the current position of the ith solution%% Step 6: Updating distance with the Sunif randrand%% h is an adaptive factor for controlling the distance between the Sun and the current planet at time th(1/(exp(n.*randn))); %% Eq.(27)%% An verage vector based on three solutions: the Current solution, best-so-far solution, and randomly-selected solutionXm(Positions(b,:)Sun_PosPositions(i,:))/3.0;Positions(i,:)Positions(i,:).*U1(Xmh.*(Xm-Positions(a,:))).*(1-U1); %% Eq.(26)else%% Step 3: Calculating an object? velocity% A flag to opposite or leave the search direction of the current planetif rand0.5 %% Eq.(18)f1;elsef-1;endL(M*(MS(i)m(i))*abs((2/(R(i)eps))-(1/(a1(i)eps))))^(0.5); %% Eq.(15)Urdrand(1,dim); %% A binary vectorif Rnorm(i)0.5 %% Eq.(13)M(rand.*(1-r)r); %% Eq.(16)lL*M*U; %% Eq.(14)Mv(rand*(1-rd)rd); %% Eq.(20)l1L.*Mv.*(1-U);%% Eq.(19)V(i,:)l.*(2*rand*Positions(i,:)-Positions(a,:))l1.*(Positions(b,:)-Positions(a,:))(1-Rnorm(i))*f*U1.*rand(1,dim).*(ub-lb); %% Eq.(13a)elseU2randrand; %% Eq. (22)V(i,:)rand.*L.*(Positions(a,:)-Positions(i,:))(1-Rnorm(i))*f*U2*rand(1,dim).*(rand*ub-lb); %% Eq.(13b)end %% End IF%% Step 4: Escaping from the local optimum% Update the flag f to opposite or leave the search direction of the current planetif rand0.5 %% Eq.(18)f1;elsef-1;end%% Step 5Positions(i,:)((Positions(i,:)V(i,:).*f)(Fg(i)abs(randn))*U.*(Sun_Pos-Positions(i,:))); %% Eq.(25)end %% End If%% Return the search agents that exceed the search spaces boundsif randrandfor j1:size(Positions,2)if Positions(i,j)ub(j)Positions(i,j)lb(j)rand*(ub(j)-lb(j));elseif Positions(i,j)lb(j)Positions(i,j)lb(j)rand*(ub(j)-lb(j));end %% End Ifend %% End ForelsePositions(i,:) min(max(Positions(i,:),lb),ub);end %% End If%% Test suites of CEC-2014, CEC-2017, CEC-2020, and CEC-2022% Calculate objective function for each search agentPL_Fit1feval(Positions(i,:));% Step 7: Elitism, Eq.(30)if PL_Fit1PL_Fit(i) % Change this to for maximization problemPL_Fit(i)PL_Fit1; %% Update the best-so-far solutionif PL_Fit(i)Sun_Score % Change this to for maximization problemSun_ScorePL_Fit(i); % Update the best-so-far scoreSun_PosPositions(i,:); % Update te best-so-far solutionendelsePositions(i,:)O_P;end %% End IFtt1; %% Increment the current function evaluationif tTmax %% Checking the termination conditionbreak;end %% End IFConvergence_curve(t)Sun_Score; %% Set the best-so-far fitness value at function evaluation t in the convergence curveend %% End for iend %% End whileConvergence_curve(t-1)Sun_Score;end%% End Function​ 参考文献往期回顾扫扫下方二维码