基于粒子群优化算法的阈值分割方法

最新推荐文章于 2025-11-24 15:27:16 发布

原创最新推荐文章于 2025-11-24 15:27:16 发布 · 1k 阅读

28 ·

CC 4.0 BY-SA版权

文章标签：

#算法 #人工智能

一、算法原理与数学模型

1. PSO阈值分割框架

$信息熵\text{分割目标}：\min_{\tau} f(\tau) = \text{类间方差} \ \text{或} \ \text{信息熵}$

粒子表示：每个粒子位置 $x_i$ 对应候选阈值 $τ\tau$
速度更新：
$vi(k+1)=w⋅vi(k)+c1r1(pbesti−xi(k))+c2r2(gbest−xi(k))v_{i}^{(k+1)} = w \cdot v_{i}^{(k)} + c_1 r_1 (pbest_i - x_i^{(k)}) + c_2 r_2 (gbest - x_i^{(k)})$
位置更新：
$x_{i}^{(k+1)} = x_i^{(k)} + v_{i}^{(k+1)}$
（ $w$ ：惯性权重， $c_1,c_2$ ：加速因子， $r_1,r_2$ ：[0,1]随机数）

2. 适应度函数设计

(1) Otsu准则（类间方差最大化）
$f(τ)=ω0(τ)ω1(τ)(μ0(τ)−μ1(τ))2f(\tau) = \omega_0(\tau)\omega_1(\tau)(\mu_0(\tau)-\mu_1(\tau))^2$
其中 $ω0,ω1\omega_0,\omega_1$ 为两类的概率， $μ0,μ1\mu_0,\mu_1$ 为均值

(2) 熵准则（信息最大化）
$f(τ)=HL+HR=−∑i=0τpilog⁡pi−∑i=τ+1L−1pilog⁡pif(\tau) = H_L + H_R = -\sum_{i=0}^\tau p_i \log p_i - \sum_{i=\tau+1}^{L-1} p_i \log p_i$
（ $H_L,H_R$ 为左右区域熵， $p_i$ 为灰度级概率）

二、算法实现流程

1. 初始化阶段

% 参数设置
n_particles = 30;    % 粒子数量
max_iter = 100;      % 最大迭代次数
w = 0.729;           % 惯性权重
c1 = 1.494; c2 = 1.494;% 加速因子

% 粒子初始化（阈值范围[0,255]）
particles = randi([0,255],n_particles,1);
velocities = 0.1*randn(n_particles,1);
pbest = particles;   % 个体最优
gbest = particles(1);// 全局最优

2. 迭代优化过程

for iter = 1:max_iter
    % 计算适应度
    for i = 1:n_particles
        tau = particles(i);
        fitness(i) = otsu_fitness(I,tau); % 调用适应度函数
    end
    
    % 更新个体最优
    update_idx = fitness > prev_fitness;
    pbest(update_idx) = particles(update_idx);
    
    % 更新全局最优
    [best_fitness, best_idx] = max(fitness);
    if best_fitness > gbest_fitness
        gbest = particles(best_idx);
        gbest_fitness = best_fitness;
    end
    
    % 速度与位置更新
    velocities = w*velocities + ...
        c1*rand(n_particles,1).*(pbest - particles) + ...
        c2*rand(n_particles,1).*(gbest - particles);
    particles = particles + velocities;
    
    % 边界处理
    particles = max(0,min(255,particles));
end

3. 关键函数实现

function fitness = otsu_fitness(I,tau)
    % 计算类间方差
    [counts,~] = imhist(I);
    total = sum(counts);
    w0 = sum(counts(1:tau+1))/total;
    w1 = sum(counts(tau+2:end))/total;
    
    mu0 = sum((0:tau)' .* counts(1:tau+1)) / (w0*total);
    mu1 = sum((tau+1:255)' .* counts(tau+2:end)) / (w1*total);
    
    fitness = w0*w1*(mu0 - mu1)^2;
end

三、算法优化策略

1. 自适应参数调整

动态惯性权重：
$w=wmax−itermax_iter(wmax−wmin)w = w_{max} - \frac{iter}{max\_iter}(w_{max}-w_{min})$
（初始 $w_{max}=0.9$ , $w_{min}=0.4$ ）
加速因子调整：
$c1=c1_init⋅e−λiterc_1 = c_{1\_init} \cdot e^{-\lambda iter}$
（ $λ=0.05\lambda=0.05$ 控制收敛速度）

2. 混合优化策略

K-Means+PSO初始化：

% K-Means预分割
[idx,centers] = kmeans(im2double(I),3);
% PSO初始化粒子
particles = centers(:,1)*255;

遗传算法交叉：引入SBX交叉算子增强多样性

function child = sbx_crossover(parent1,parent2,pc)
    if rand < pc
        eta = 2; % 分布指数
        u = rand(size(parent1));
        beta = (2*u).^(1/(eta+1));
        child = 0.5*((1+beta).*parent1 + (1-beta).*parent2);
    else
        child = parent1;
    end
end

3. 并行计算加速

% 使用parfor并行计算适应度
parfor i = 1:n_particles
    fitness(i) = otsu_fitness(I,particles(i));
end

六、Matlab代码

1. 完整代码框架

function [threshold,seg] = pso_threshold_segmentation(I)
    % 参数设置
    n_particles = 30;
    max_iter = 100;
    w = 0.7; c1=1.5; c2=1.5;
    
    % 初始化粒子群
    particles = 255*rand(n_particles,1);
    velocities = 0.1*randn(n_particles,1);
    pbest = particles;
    [gbest,gbest_fitness] = update_global_best(I,particles);
    
    % 迭代优化
    for iter = 1:max_iter
        % 更新速度位置
        velocities = w*velocities + ...
            c1*rand(n_particles,1).*(pbest - particles) + ...
            c2*rand(n_particles,1).*(gbest - particles);
        particles = particles + velocities;
        particles = max(0,min(255,particles));
        
        % 更新最优解
        [current_best,fitness] = update_global_best(I,particles);
        if fitness > gbest_fitness
            gbest = current_best;
            gbest_fitness = fitness;
        end
        
        % 显示进度
        fprintf('Iter %d: Best Fitness=%.4f\n',iter,gbest_fitness);
    end
    
    % 输出结果
    threshold = gbest;
    seg = I > threshold;
end

function [best,fitness] = update_global_best(I,particles)
    n = numel(particles);
    fitness = zeros(n,1);
    for i=1:n
        fitness(i) = otsu_fitness(I,particles(i));
    end
    [best_fitness, idx] = max(fitness);
    best = particles(idx);
end