2023年7月29日发(作者：)

压缩感知正交匹配追踪算法重构二维图像

摘要

在传统采样过程中，为了避免信号失真，采样频率不得低于信号最高频率的2倍。然而，对于数字图像、视频的获取，依照香农定理会导致海量的采样数据，大大增加了存储和传输的代价。压缩感知采用非自适应性投影来保持信号的原始结构，能够通过数值最优化问题准确重构原始信号。该理论指出，如果信号是稀疏的或者在某个基下可压缩，那么用少量的观测值就可以保持信号的结构和相关信息。基于该理论，用于精确重构信号的采样需求数量可以远低于观测的维度，这极大地缓解了宽带信号处理的压力。正交匹配追踪算法正是压缩感知信号检测的一种算法。本文将介绍正交匹配追踪算法的原理以并给出了测试效果。

一、压缩感知简介

压缩感知是一种新的信息获取理论，是建立在信号稀疏表示、测量矩阵的非相关性以及逼近理论上的一种信号采集和重建的方法。该理论指出，只要信号是稀疏的或者在某个基下时刻压缩的，就可以通过远低于奈奎斯特采样定理要求的采样率获取信号的结构信息，再通过重构算法完成信号的精确重构。

压缩感知理论只要包括两个部分：将信号在观测向量上投影得到观测值，以及利用重构算法由观测值重构信号。

设x是一个长度为N的信号，其稀疏度为KKN，系数度为K指x本身有K个非零元素，或者在某种变化域内的展开系数有K个非零元素。

信号（假设信号在变换域内K系数）在观测向量上的投影可以表示为：

yii,x,i1,,M,MN

其中，y为压缩感知获取的M个采样值，的观测基与变换基不相关。

iM是一组观测向量，由i1iM组成i1重构信号的关键是找出信号x在域中的稀疏表示，可以通过l0范数优化问题找到具有系数结构的解：

minTxs.t.0yx

由于上式的优化问题是一个难求解的NP-hard问题，所以可以用l1约束取代l0约束：

min1yx

此时，压缩感知获得的采样值已经保持了原信号的结构及相关信息，因此可以不需要重构信号，利用检测算法直接从采样值中提取特征量进行判断，完成信号检测任务。

二、正交匹配追踪算法

1.最小l0范数模型

从数学意义上讲，基于压缩感知理论的信号重建问题就是寻找欠定方程组(程的数量少于待解的未知数)的最简单解的问题，因l0范数刻画得就是信号中非零元素的个数，而能够使得结果尽可能地稀疏。通常我们采用下式描述最小l0范数最优化问题：

YX (3.1)

实际中，允许一定程度的误差存在，因此将原始的最优化问题转化成一个较简单的近似形式求解，其中是一个极小的常量：

YX022 (3.2)

但是这类问题的求解数值计算极不稳定，很难直接求解。

匹配追踪类稀疏重建算法解决的是最小l0范数问题，最早提出的有匹配追踪(MP)算法和正交匹配追踪(OMP)算法。

2.匹配追踪算法

匹配追踪算法的基本思想是在每一次的迭代过程中，从过完备原子库里(即感知矩阵)选择与信号最匹配的原子来进行稀疏逼近并求出余量，然后继续选出与信号余量最为匹配的原子。经过数次迭代，该信号便可以由一些原子线性表示。但是由于信号在己选定原子(感知矩阵的列向量)集合上的投影的非正交性使得每次迭代的结果可能是次最优的，因此为获得较好的收敛效果往往需要经过较多的迭代次数。

匹配追踪类算法通过求余量r与感知矩阵中各个原子之间内积的绝对值，来计算相关系数u：

uuj|ujr,j,j1,2,…，N

并采用最小二乘法进行信号逼近以及余量更新:

ˆargminYXXiR2

ˆ,

rnewYX3.正交匹配追踪算法

正交匹配追踪算法（Orthogonal Matching Pursuit,OMP），是最早的贪婪迭代算法之一。该算法沿用了匹配追踪算法中的原子选择准则，只是通过递归对己选择原子集合进行正交化以保证迭代的最优性，从而减少迭代次数。OMP算法则有效克服了匹配追踪算法为获得较好的收敛效果往往需要经过较多的迭代次数的问题。

OMP算法将所选原子利用Gram-Schmidt正交化方法进行正交处理，再将信号在这些正交原子构成的空间上投影，得到信号在各个已选原子上的分量和余量，然后用相同方法分解余量。在每一步分解中，所选原子均满足一定条件，因此余量随着分解过程迅速减小。通过递归地对已选择原子集合进行正交化保证了迭代的最优性，从而减少了迭代次数。

OMP的重建算法是在给定迭代次数的条件下重建，这种强制迭代过程停止的方法使得OMP需要非常多的线性测量来保证精确重建。总之，它以贪婪迭代的方法选择的列，使得在每次迭代中所选择的列与当前的冗余向量最大程度地相关，从测量向量中减去相关部分并反复迭代，直到迭代次数达到稀疏度K，强制迭代停止。

OMP算法的具体步骤如下：

(1)初始余量r0Y，迭代次数n1，索引值集合，J；

(2)计算相关系数u ,并将u中最大值对应的索引值存入J中；

(3)更新支撑集，其中J0；

ˆ，同时用式(3.4)对余量进行更新； (4)应用式(3.3)得到X(5)若rnewr2，令rrnew，nn1，转步骤(2)；否则，停止迭代。

三、正交匹配追踪算法的Matlab实现

实验用Matlab编写OMP算法程序，并用数字图像处理标准图像“Lena256”进行了检验，实验程序及结果如下。

Matlab程序：

-----------------------------------------------------------------

% 图像读取变换

X=imread('');% 读文件

X=double(X);

[a,b]=size(X);

% 小波变换矩阵生成

ww=DWT(a);

% 小波变换让图像稀疏化

X1=ww*sparse(X)*ww';

X1=full(X1);

M=190;% 随机矩阵生成 R=randn(M,a);

Y=R*X1;% 测量

-----------------------------------------------------------------

% OMP算法

X2=zeros(a,b); % 恢复矩阵

for i=1:b % 列循环

rec=omp(Y(:,i),R,a);

X2(:,i)=rec;

end

----------------------------------------------------------------

%绘图

figure(1);% 原始图像

imshow(uint8(X));

title('原始图像');

figure(2);% 变换图像

imshow(uint8(X1));

title('小波变换后的图像');

figure(3);% 压缩传感恢复的图像

X3=ww'*sparse(X2)*ww; % 小波反变换

X3=full(X3);

imshow(uint8(X3));

title('恢复的图像');

-----------------------------------------------------------------

% 误差(PSNR)

errorx=sum(sum(abs(X3-X).^2)); % MSE误差

psnr=10*log10(255*255/(errorx/a/b)) % PSNR

-----------------------------------------------------------------

% OMP的函数

% s-测量；T-观测矩阵；N-向量大小

function hat_y=omp(s,T,N)

Size=size(T); % 观测矩阵大小

M=Size(1); % 测量

hat_y=zeros(1,N); % 待重构的谱域(变换域)向量

Aug_t=[]; % 增量矩阵(初始值为空矩阵)

r_n=s; % 残差值

for times=1:M/4; % 迭代次数(稀疏度是测量的1/4)

for col=1:N; % 恢复矩阵的所有列向量

product(col)=abs(T(:,col)'*r_n); % 恢复矩阵的列向量和残差的投影end % 系数(内积值)

[val,pos]=max(product); % 最大投影系数对应的位置

Aug_t=[Aug_t,T(:,pos)]; % 矩阵扩充

T(:,pos)=zeros(M,1); % 选中的列置零

aug_y=(Aug_t'*Aug_t)^(-1)*Aug_t'*s; % 最小二乘,使残差最小

r_n=s-Aug_t*aug_y; % 残差

pos_array(times)=pos; % 纪录最大投影系数的位置

if (norm(r_n)<9) % 残差足够小

break;

end

end

hat_y(pos_array)=aug_y; % 重构的向量

----------------------------------------------------------------------------------------------------------------------

% 构造正交小波变换矩阵函数，图像大小N*N，N=2^P，P是整数。

function ww=DWT(N)

[h,g]= wfilters('sym8','d'); % 分解低通和高通滤波器

% N=256; % 矩阵维数(大小为2的整数幂次)

L=length(h); % 滤波器长度

rank_max=log2(N); % 最大层数

rank_min=double(int8(log2(L)))+1; % 最小层数

ww=1; % 预处理矩阵

for jj=rank_min:rank_max% 矩阵构造

nn=2^jj;

p1_0=sparse([h,zeros(1,nn-L)]);% 构造向量

p2_0=sparse([g,zeros(1,nn-L)]);

for ii=1:nn/2 % 向量圆周移位

p1(ii,:)=circshift(p1_0',2*(ii-1))';

p2(ii,:)=circshift(p2_0',2*(ii-1))';

end

w1=[p1;p2];% 构造正交矩阵

mm=2^rank_max-length(w1);

w=sparse([w1,zeros(length(w1),mm);zeros(mm,length(w1)),eye(mm,mm)]);

ww=ww*w;

clear p1;clear p2;

end

图像结果：