盲嵌入算法

关键词： 原始 嵌入 算法 水印

盲嵌入算法（精选八篇）

盲嵌入算法 篇1

近年来相关领域提出了自嵌入技术对鲁棒水印算法的研究。参考文献[6]引入自嵌入技术, 提出一种DWT-SVD域全盲鲁棒量化水印算法, 虽具有良好的不可见性和安全性, 但是该算法中水印序列自嵌入对含水印图像提取特征水印序列的最大值有影响, 因此提取的嵌入水印序列与认证水印序列之间存在误差。参考文献[7]将自嵌入思想引入到鲁棒水印领域, 研究表明该算法在抵抗平滑、添加噪声、JPEG压缩、重采样、剪切和几何攻击上表现出很强的鲁棒性, 但是该算法将原始图像每一子块产生的特征水印嵌入到该块的DCT系数中, 使得算法对矢量量化等攻击鲁棒性较差。

本文通过分析与比较, 针对前述算法的不足, 提出了一种基于非负矩阵分解和提升小波变换的自嵌入盲检测鲁棒数字水印算法。详细介绍了该算法的设计与分析, 最后通过实验表明该方案对常见信号处理具有很强的鲁棒性, 实现了鲁棒水印序列完全盲检测。

1 算法的分析与设计

非负矩阵分解 (NMF) 和提升小波是本文算法的设计基础。对NMF的主要介绍可以见参考文献[8-9], NMF在某种程度上体现了智能化数据描述的本质。

提升 (二代) 小波构造方法由Sweldens提出, 与传统小波变换相比, 提升格式小波变换不仅克服了传统小波计算量大、不能精确重构原始信号等缺点, 而且具有包容传统小波、运算速度快、允许完全原位计算、易于实现并行运算等优点。

基于上述的理论基础, 本文提出了一种自嵌入盲检测鲁棒数字水印算法, 算法描述中使用了非负矩阵分解与提升小波的相关术语与符号。该算法包括鲁棒水印序列提取、鲁棒水印序列自嵌入、鲁棒水印序列检测以及认证水印序列提取4个步骤, 下面进行详细说明。

1.1 鲁棒水印序列提取

为便于鲁棒水印序列提取, 本文使用小波低频逼近系数均值作为统计特征。该统计特征从逼近信号的小波系数得到, 代表原始信号感知上最重要的分量, 对于常见信号处理操作具有较强的鲁棒性。而且由于相邻图像像素间具有高度相关性, 随机剪切等操作即使引起个别小波系数发生较大的改变, 也不会使统计平均值发生很大变化。

本文设计算法的鲁棒水印序列W赞提取过程如下所述:

(1) 将原始图像载体信号I分割为L个互不重叠大小为n×n的子块Ii{i=0, 1, …, L}, 分块个数决定鲁棒水印序列W赞的数据元素个数。

(2) 依次对每一图像子块Ii{i=0, 1, …, L}进行H级小波分解, 选取子块图像小波低频近似分量C (i) Hlow。

(3) 将图像子块低频近似分量C (i) Hlow进行NMF分解, 得到其NMF分解基矩阵UA (i) 和系数矩阵VA (i) 。

(4) 系数矩阵VA (i) 块内求和得V (i) , 并计算均值V=Mean (V (i) ) 作为阈值, 根据V (i) 中元素与V的大小关系产生鲁棒水印序列W赞, 即:

1.2 鲁棒水印序列自嵌入

本文设计算法的鲁棒水印序列自嵌入过程描述如下:

(1) 将原始图像载体信号I分割为L个互不重叠大小为n×n的子块Ii{i=0, 1, …, L}。

(5) 对嵌入水印序列后的小波低频近似分量进行逆向LWT变换, 得到含水印图像IW。

1.3 鲁棒水印序列检测

详细水印序列检测过程如下所述:

(3) 根据下式提取鲁棒水印序列:

1.4 认证水印序列提取

本文设计算法的认证水印序列W"提取过程如下所述:

2 实验结果

2.1 实验参数说明

仿真实验选用的标准灰度图像Lena、Goldhill、Baboon等作为原始载体图像, 实验平台使用Windows XP操作系统和Matlab 7.0编程环境, CPU为Intel CPU E65002.9 GHz, 内存为4 GB。原始载体进行以Daubchies5/3小波为基的1级LWT, 小波低频子带系数进行NMF以产生原始特征水印序列。自嵌入水印序列的长度为1 024 bit。Lena、Goldhill、Baboon 3幅图像的量化步长都为21。

2.2 不可见性实验

分别从主观与客观测评两个方面衡量算法的不可见性。图1 (a) 、图1 (c) 、图1 (e) 为原始标准测试图像, 图1 (b) 、图1 (d) 、图1 (f) 为本文算法下的含水印图像。从图1的实验结果可见, 本文设计的算法具有良好的不可见性。客观测评采用峰值信噪比PSNR (Peak Signal to Noise Ratio) , 定义如下:

根据本文算法嵌入水印前后标准测试图像与含水印图像间的PSNR与参考文献[6, 10, 11]的比较如表1所示。从表1可见, 本文算法与同类文献中算法均具有较好的视觉效果。

2.3 抗攻击鲁棒性实验结果

使用归一化相关系数NC (Normalized Correlation) 作为算法稳健性衡量指标。本文算法在各种常见图像处理攻击和组合攻击下的实验结果如表2所示。从表2可见, 本文算法对常见图像处理攻击和组合攻击均具有较强的鲁棒性, 部分性能优于对比参考文献。

2.4 量化步长选择

量化步长η的选择决定着算法的鲁棒性和不可见性。η越小, 算法的不可见性越好, 但鲁棒性越差;η越大, 算法的鲁棒性越好, 但不可见性越差。图2给出了本文算法PSNR与量化步长之间的关系曲线图。从图2可见, 当3幅图像的量化步长η取21时, 本文算法具有良好不可见性和较强鲁棒性。

本文提出了一种自嵌入盲检测鲁棒数字水印算法。实验结果表明, 该方案对常见信号处理具有很强的鲁棒性, 实现了鲁棒水印序列完全盲检测。此外, 本算法还具有计算简单、容易实现等优点, 大大增强了其用于数字图像作品版权保护的实用性, 具有一定的应用价值。今后的研究将从进一步深入探讨NMF特征序列的稳定性以及分析LWT变换系数特性, 设置更合适的量化步长, 进一步提高算法整体性能等方面展开。

摘要：为了实现鲁棒水印算法的完全盲检测, 提出一种自嵌入盲检测鲁棒数字水印算法。该算法对原始图像分块子图进行小波提升变换得到低频近似分量, 对近似分量作非负矩阵分解得到可近似表示子块图像的基矩阵和系数矩阵;在此基础上将系数矩阵量化得到鲁棒水印序列, 最后鲁棒水印序列自适应量化嵌入原始图像低频自带近似分量。实验结果表明, 该方案对常见信号处理具有很强的鲁棒性, 实现了鲁棒水印序列完全盲检测。

关键词：数字水印,自嵌入,完全盲检测,提升小波变换,非负矩阵分解

参考文献

[1]Wang Xiongyang, Hou Limin, Wu Jun.A feature-based robust digital image watermarking against geometric attacks[J].Image and Vision Computing, 2008, 26 (7) :980-989.

[2]Wang Xiongyang, Cui Changying.A novel image watermarking scheme against desynchronization attacks by SVR revision[J].Journal of Visual Communication and Image Representation, 2008, 19 (5) :334-342.

[3]Li Leida, Qian Jiansheng, PAN J S.Characteristic region based watermark embedding with RST invariance and high capacity[J].International Journal of Electronics and Communications, 2011, 65 (5) :435-442.

[4]左姣, 范九伦.自嵌入全盲鲁棒水印量化算法[J].微型机与应用, 2014, (2) :38-40, 43.

[5]冯祥斌, 陈永红.基于纠错编码和RBF神经网络的盲水印新算法[J].微型机与应用, 2013 (7) :30-33.

[6]叶天语.DWT-SVD域全盲自嵌入鲁棒量化水印算法[J].中国图象图形学报, 2012, 17 (6) :644-650.

[7]叶天语.自嵌入完全盲检测顽健数字水印算法[J].通信学报, 2012, 33 (10) :7-15.

[8]刘竞杰, 陶亮.一种基于非负矩阵分解的鲁棒零水印算法[J].计算机工程与应用, 2012, 48 (16) :90-93, 106.

[9]DONOHO D, STODDEN V.When does non-negative matrix factorization give a correct decomposition into parts[C].Proceeding of Advances in Neural Information Processing Systems, Cambridge:MIT Press, 2004:1141-1148.

[10]李旭东.基于奇异值分解的灰度级数字水印算法[J].武汉大学学报:信息科学版, 2010, 35 (11) :1305-1308, 1359.

盲嵌入算法 篇2

分布嵌入式大气数据系统算法的初步研究

针对某型飞机给出了DFADS系统的构型描述,进行了DFADS的算法设计,算法采用了一种非物理映射的方法来建立各测压点压力和基本大气参数之间的关系.先是利用制定的多传感器数据表决规则表决出2～3对相关的`测压点,然后通过建立相关测压点压力和其压力系数值之间的关系,利用这几对测压点的压力对事先建立的Cpi(α,β,Ma)映射关系表格数据库进行查表计算,从而得到当前的飞行状态.本文利用计算流体动力学(CFD)计算手段获得了332个飞行状态下飞机前机身14个测压点的表面压力系数数据,并以此为基础对DFADS的算法进行了仿真验证.结果表明,该DFADS的算法可以根据14个点的压力输入,正确地解算当前的大气参数.

作 者：王岩 郑伟 WANG Yan ZHENG Wei 作者单位：沈阳飞机设计研究所,辽宁,沈阳,110035刊 名：飞机设计英文刊名：AIRCRAFT DESIGN年，卷(期)：28(6)分类号：V241.7关键词：嵌入式大气数据系统 多传感器 数据表决 映射嵌入式

超完备瞬时盲分离算法研究 篇3

独立分量分析[1] (Independent Component Analysis, ICA) 是一种通过最大化多维观测向量的统计独立性来寻找合适的线性变换的统计方法。其目标是在源信号向量和传输信道均未知的情况下, 从观测数据中重建和恢复各独立源。独立分量分析在盲源分离 (Blind Source Separation, BSS) 方面的应用, 使得BSS的研究取得了突破性的进展, 发展了很多优秀的算法。然而传统的ICA是要求观测信号数目大于或等于源信号数目的, 现实中存在很多源信号数目大于观测信号数目的情况, 这种情况下的盲源分离被称为超完备盲分离[1,2,3,4,5,6,7,8,9,10,11,12,13] (over-complete BSS) 。

由于观测信号数目少于源信号数目, 相当于信源在经过混合信道后, 发生了有损压缩, 因此采用传统的ICA通过对混合系统求伪逆的过程已无法恢复出源信号, 这些丢失的信息无法仅通过数学手段来恢复, 只能通过一些先验、假设或限制条件 (如:独立性、稀疏性等) 进行弥补。

本文就超完备瞬时盲分离的问题描述和目前的主要解决方案做了一个详细的介绍, 并依据不同算法的特点, 将之归纳为两类, 总结了各自的优缺点及其适用条件。

1 问题描述

盲源分离的数据模型如下:

$\mathit{x}=\mathit{A}\mathit{s}+\mathit{\epsilon }(1)$

式中:s=[s1, s2, …, sN]T是N维未知独立源信号矢量, 经过线性系统A混合, 再叠加噪声信号ε=[ε1, ε2, …, εM]后, 得到M维的观测信号矢量x=[x1, x2, …, xM]T, 线性系统A为M×N未知混合矩阵, 在超完备情况下, 假设M<N, 即源信号数目大于观测信号数目。

完备情况下的盲分离一般可以同时估计出混合矩阵和源信号, 但在超完备情况下, 同时估计源和混合矩阵十分困难, 所以分离过程一般分为两步, 先估计混合矩阵, 再由估计出的混合矩阵和观测信号去估计源信号。也可以在同一迭代过程中交替地估计源信号和混合矩阵。因此根据各算法的估计特点, 可以将各算法分为两步估计法和交替估计法两类, 下面将分别介绍这两类算法。

2 两步估计法[7,8,9,10]

2002年, Theis F J在文献[10]及其他学者研究成果的基础上, 提出了两步实现超完备盲分离的思想[7], 即BMMR和BSR, 先由观测信号估计混合矩阵, 再在观测信号和估计出的混合矩阵基础上估计源信号。这种方法多基于无噪数据模型, 即x=As, 各变量定义同式 (1) 。

2.1 BMMR

混合矩阵的盲恢复可以采用聚类[10,12], 几何ICA[8,9]等方法来实现。聚类算法一般要求源信号比较稀疏, 可以近似地认为在每一时刻仅有一个信号非零, 而采用几何ICA算法对稀疏性并不做特别要求。

2.1.1 聚类算法

采用无噪数据模型x=As。以M=2为例介绍该算法, 系统模型可以写为:

$\left[\begin{array}{l}x{}_1\\ x{}_2\end{array}\right]=\left[\begin{array}{ccc}a{}_{11}& a{}_{12}& \cdots \\ a{}_{21}& a{}_{21}& \cdots \end{array}\right]\left[\begin{array}{l}s{}_1\\ s{}_2\\ ⋮\\ s{}_{{\rm N}}\end{array}\right](2)$

一般假设s1, s2, …, sN具有稀疏分布, 那么在x1, x2中它们常常只是单独出现。如果某一时刻t只有s1有值, 则有x${}_1^t$=a11s${}_1^t$, x${}_2^t$=a21s${}_1^t$, 因而此时x${}_1^t$/x${}_2^t$=a11/a21。可见, 属于这种情况的观测数据x在x1～x2散点图上将聚集在斜率为a11/a21的斜线上 (见图1) 。同理, 当只有si出现的时刻, x1～x2散点图上将聚集在斜率为a1i/a2i的斜线上。那么在其散点图上将会聚集N条线, 根据这些线的斜率, 再结合混合矩阵A各列的归一化条件|ai|2=1 (i=1, 2, …, N) , 就可以把矩阵A的各列向量估计出来。

推而广之, 当观测信号数目M>2时, 观测数据的散点图将聚集在各ai方向上, 因此求各ai相当于找数据在M维空间的聚集位置 (这可以通过神经网络的C均值聚类/K均值聚类来实现) 。这些ai即为混合矩阵A各列的矢量。

2.1.2 几何ICA算法

该算法是采用了“赢家通吃” (winner-takes-all) 规则, 可以描述如下:在单位球SM-1⊂RM上选择2N个单位根ω1, ω′1, …, ωN, ω′N, 使得ωi=-ω′i (i=1, 2, …, N) , 且每对根彼此线性独立, 这些根称为神经元。取固定的学习率η:N→R, 使得满足$\eta (n)>0‚{\displaystyle \sum _{n\in {\rm N}}\eta }(n)=\infty$和${\displaystyle \sum _{n\in {\rm N}}\eta }(n){}^2<\infty$。然后按照下面的步骤进行迭代, 直到遇到收敛条件。

选择随机变量x中的一个样本x (t) ∈R, 若x (t) =0, 则重新选择一个非零的样本 (当x的联合概率密度函数ρx为连续时将不会出现这种情况) , 将x (t) 投影到单位球, 得到:y (t) =x (t) /|x (t) |。假设ωi或ω′i是所有根中根据欧氏距离算出的离y (t) 最近的点, 则进行如下迭代:

ωi (t+1) =π{ωi (t) +η (t) sgn[y (t) -ωi (t) ]} (3)

式中:π:RM�}→SM-1表示在RM中的M-1单位球SM-1上的投影。另外有:

$\omega {}_i^{´}(t+1)=-\omega {}_i(t+1)(4)$

其他的神经元在本次迭代中则不进行更新。

所有的神经元均按照这种迭代方式进行更新。迭代的最终结果将获得一组ω1, ω′1, …, ωN, ω′N, 将ω1, ω2, …, ωN组成向量矩阵A, 即为待恢复的混合矩阵。

2.1.3 算法比较

聚类算法要求源信号满足一定的稀疏性, 但随着信号维数的增加, 计算代价不会迅速增长。

几何ICA算法从信号的几何特性出发, 便于理解, 所以得到了更广阔的应用。它对信号的稀疏特性没有特别的要求。且Theis F J在文献[8]中给出了矩阵等价的概念和算法收敛的条件, 并严格证明了分离的混合矩阵的惟一性。几何ICA的主要缺点就是随着观测信号维数的增长, 样本和收敛次数呈指数增长, 因此几何ICA一般都限制在低维情况。

2.2 BSR

在超完备基条件下, 就算己知混合矩阵A, 源信号的分离也是不惟一的, 因此常常假设源信号满足一定的稀疏性, 然后利用稀疏分量分析理论, 根据最大稀疏性准则来分离源信号, 最短路径算法和广义FOCUSS算法都是以寻求信号的最大稀疏性表示为目标的源恢复方法, 最短路径算法是基于单次观测样本进行处理的, 需要的样本数据较多, 但单次计算的代价较低, 而FOCUSS则是对信号的所有样本进行批处理的, 效率较高, 计算代价相对较高。

2.2.1 最短路径算法[7,10,12]

最短路径算法是以“L1范数最小”为目标的一种分解方法, 此时要求:

$\mathit{x}{}^t={\displaystyle \sum _{j=1}^{{\rm N}}\mathit{a}}{}_{j}s{}_j^t‚\text{且}{\displaystyle \sum _{j=1}^{{\rm N}}|}s{}_j^t|=\min (5)$

上式前一部分说明, 每个可行解都是各矢量ajs${}_j^t$的矢量和, 而且其总和等于xt, 上式的后一部分说明, 在全部可行解中寻找路径 (矢量长度之和) , 最短的一条就是所求答案。

以M=2为例, 说明最短路径寻优方法, 如图2所示。

(1) 计算在某时刻t时xt的矢量方向θt=arctan (x${}_2^t$/x${}_1^t$) 。

(2) 选择各aj (j=1, 2, …, N) 中方向居于θt两侧, 且与θt最接近的两个矢量作为基矢量 (图中的aA与aB) 。

(3) 以选定的两个基矢量为方向构造平行四边形, 使其合成矢量为xt, 从而确定矢量的长度s${}_A^t$和s${}_B^t$。也就是说此时的分解结果为xt=aAs${}_A^t$+aBs${}_B^t$, 从而得到t时刻的源信号表示:

$(s{}^t){}_i=\{\begin{array}{l}[(\mathit{a}{}_A|\mathit{a}{}_B){}^{-1}\mathit{x}{}^t]{}_A‚i=A\\ [(\mathit{a}{}_A|\mathit{a}{}_B){}^{-1}\mathit{x}{}^t]{}_B‚i=B\\ 0‚\text{o}\text{t}\text{h}\text{e}\text{r}\text{w}\text{i}\text{s}\text{e}\end{array}(6)$

从而把t时刻的观测信号xt分解成只含有两个信源的稀疏组合。

(4) 依次对t=1～T逐点进行, 便将x分解成若干稀疏信源的组合。

aA和aB是A中最靠近xt的单位矢量, xt=aAs${}_A^t$+aBs${}_B^t$, 而且|s${}_A^t$|+|s${}_B^t$|为最小。

2.2.2 似p范数代价函数和广义FOCUSS算法[11,12]

分散性可以看作稀疏性的反义词, 很多文献认为“似p范数”可以作为分散性的度量:

$J{}_p(\mathit{s})={\displaystyle \sum _{i=1}^{{\rm N}}\left|s{}_i\right|}{}^p,0\le p\le 1(7)$

当p=0时, J0 (s) 将等于s中的非零元素的个数;p=1时, J1 (s) 将等于s中所有元素的绝对值之和。Jp (s) 愈大说明能量越分散, 所以最小化Jp (s) 将会使s的能量越集中, 稀疏性越大。因此为了在等式约束x=As下最小化式 (7) , 定义Lagrange函数L (s, λ) 。

$L(\mathit{s},\mathit{\lambda })=J{}_p(\mathit{s})+\mathit{\lambda }(\mathit{x}-\mathit{A}\mathit{s})(8)$

式中:λ∈RN, 为Lagrange乘子矢量, 那么式 (8) 的平衡点可以估计如下:

$\{\begin{array}{l}\nabla {}_{\mathit{s}}L(\mathit{s}{}^{*},\mathit{\lambda }{}^{*})=\nabla {}_{\mathit{s}}J{}_p(\mathit{s})-\mathit{A}{}^{\text{Τ}}\mathit{\lambda }{}^{*}=0\\ \nabla {}_{\mathit{\lambda }}L(\mathit{s}{}^{*},\mathit{\lambda }{}^{*})=\mathit{x}-\mathit{A}\mathit{s}{}^{*}=0\end{array}(9)$

式中:ᐁsJp (s) =pD-1 (s) s, 其中D (s) =diag$[\left|s{}_i\right|{}^{p-2}]$是对角矩阵, 解方程 (9) 可得:

$\begin{array}{l}\mathit{\lambda }{}^{*}=p[\mathit{A}\mathit{D}(\mathit{s}{}^{*})\mathit{A}{}^{\text{Τ}}]{}^{-1}\mathit{x}(10)\\ \mathit{s}{}^{*}=p{}^{-1}\mathit{D}(\mathit{s}{}^{*})\mathit{A}{}^{\text{Τ}}\mathit{\lambda }{}^{*}=\mathit{D}(s{}^{*})\mathit{A}{}^{\text{Τ}}[\mathit{A}\mathit{D}(\mathit{s}{}^{*})\mathit{A}{}^{\text{Τ}}]{}^{-1}\mathit{x}(11)\end{array}$

从而可以得到估计最优矢量s*的广义FOCUSS迭代算法为:

$\mathit{s}(k+1)=\mathit{D}[\mathit{s}(k)]\mathit{A}{}^{\text{Τ}}\{\mathit{A}\mathit{D}[\mathit{s}(k)]\mathit{A}{}^{\text{Τ}}\}{}^{-1}\mathit{x}(12)$

3 交替估计法

这类算法的特点是源信号和混合矩阵的估计过程交替进行, 直到收敛。它通常以最大化似然函数为目标, 通过梯度寻优的方式实现对混合矩阵的估计, 再在某一假设源信号概率分布的情况下采用最大化后验概率的方法估计源信号, 常采用的源信号模型有拉普拉斯分布[3,4,5]、高斯混合分布[6]、广义高斯分布等。

3.1 最大化似然函数法估计混合矩阵

对于式 (1) 定义的模型, 可以得到观测数据x的似然函数为:

$L=p(\mathit{x}|\mathit{A})={\displaystyle \int p}(\mathit{x}|\mathit{A},\mathit{s})p(\mathit{s})\text{d}\mathit{s}(13)$

采用最大化对数似然函数的方法来估计混合矩阵, 那么目标函数就变为:

$\begin{array}{l}\widehat{\mathit{A}}=\arg \underset{\mathit{A}}{{\displaystyle \max }}[\ln p(\mathit{x}|\mathit{A})]\\ =\arg \underset{\mathit{A}}{{\displaystyle \max }}[\ln {\displaystyle \int p}(\mathit{x}|\mathit{A},\mathit{s})p(\mathit{s})\text{d}\mathit{s}](14)\end{array}$

式中:p (x|A, s) 根据系统假设的噪声模型的不同而不同, 当噪声满足高斯分布时, 有:

$p(\mathit{x}|\mathit{A},\mathit{s})=\left(\frac{1}{\sqrt{2\text{π}}}\right){}^{{\rm M}}\frac{1}{\sigma }\exp \left(-\frac{\parallel \mathit{x}-\mathit{A}\mathit{s}\parallel {}^2}{2\sigma {}^2}\right)(15)$

当噪声假设为具有更宽泛分布的广义高斯分布时, 有:

$p(x|\mathit{A},\mathit{s})=\frac{\rho }{2\lambda \Gamma (\frac{1}{\rho })}\exp \left[-\left(\frac{|\mathit{x}-\mathit{A}\mathit{s}|}{\lambda }\right){}^{\rho }\right](16)$

式中:$\lambda =\sqrt{\frac{\Gamma (1/\rho )}{\Gamma (3/\rho )}}$;Γ (·) 表示伽玛分布。

由于通信系统中很多情况下的噪声都可以近似看作高斯噪声, 所以在后面的讨论中, p (x|A, s) 的分布一般采取式 (15) 。然而不管噪声是何种分布, 在过完备的情况下, 由于A不可逆, 使得式 (13) 中的积分计算相当复杂, 只能对上述积分进行数值近似, 常采用的近似方法有最大值近似[2]和高斯积分法则近似[3,4,5]。最大值近似法一般要求被积函数p (x|A, s) p (s) 具有相当尖锐的分布, 近似的过程忽略了大量的概率信息, 且为了保证算法的收敛, 要对A的列矢量的长度加以限制, 不然将会造成较大的近似误差。然而高斯积分法则相对准确, 对被积函数的分布没有特别要求。下面将分别介绍这两种近似方法。

3.1.1 最大值近似

1997年, Olshausen和Field提出了最大值近似方法[1], 它假设式 (13) 的被积函数p (x|A, s) p (s) 具有相当尖锐的分布, 即:有一个最大值, 其他的都很小, 于是可以仅取此最大值作为式 (13) 的积分近似值, 那么混合矩阵的估计式 (14) 就变为:

$\widehat{\mathit{A}}=\arg \underset{\mathit{A}}{{\displaystyle \max }}<\underset{\mathit{s}}{{\displaystyle \max }}\ln [p(\mathit{x}|\mathit{A},\mathit{s})p(\mathit{s})]>(17)$

式中:arg (·) 表示求值。记E (x, s|A) =-ln[p (x|A, s) p (s) ], 则由式 (15) 可知:

$E(\mathit{x},\mathit{s}|\mathit{A})=\frac{1}{2\sigma {}^2}\parallel \mathit{x}-\mathit{A}\mathit{s}\parallel {}^2-{\displaystyle \sum _{i=1}^{{\rm N}}\text{l}}\text{n}p{}_i(s{}_i)+\text{C}\text{o}\text{n}\text{s}\text{t}.(18)$

那么式 (17) 变为:

$\widehat{\mathit{A}}=\arg \underset{\mathit{A}}{{\displaystyle \min }}<\underset{\mathit{s}}{{\displaystyle \min }}E(\mathit{x},\mathit{s}|\mathit{A})>(19)$

最小化E (x, s|A) 的过程相当于最小化式 (18) 的前两项, 最小化第一项相当于减小重建向量的误差, 从而使$\widehat{\mathit{A}}$张成信号子空间, 最小化第二项使s趋于最大稀疏分布, 这两项在最小化过程中相互影响, 从而使$\widehat{\mathit{A}}$的估计达到最优值。

式 (19) 的优化过程一般分两步:首先, 固定A, 最小化E (x, s|A) , 求得$\widehat{\mathit{s}}$, 然后再固定$\widehat{\mathit{s}}$, 用梯度下降法最小化式 (18) 可得估计矩阵$\widehat{\mathit{A}}$。

3.1.2 高斯积分近似

由于高斯积分近似中有:

${\displaystyle \int f}(\mathit{x})\text{d}\mathit{x}⧋\text{c}\text{o}\text{n}\text{s}\text{t}\cdot f(\widehat{\mathit{x}})\cdot \left|-\frac{\text{d}}{\text{d}{}^2\mathit{x}}\log f(\mathit{x})|{}_{\mathit{x}=\widehat{\mathit{x}}}\right|{}^{-\frac{1}{2}}(20)$

所以式 (13) 可以近似为:

式中:H是p (x|A, s) p (s) 在$\mathit{s}=\widehat{\mathit{s}}$的Hessian矩阵:

$\mathit{{\rm H}}=-\frac{\text{d}}{\text{d}\mathit{s}\mathit{s}{}^{\text{Τ}}}\ln p(\mathit{x}|\mathit{A},\mathit{s})p(\mathit{s})|{}_{\mathit{s}=\widehat{\mathit{s}}}(22)$

将式 (15) 代入式 (21) 并对之两边取对数得:

$\begin{array}{l}\ln {\displaystyle \int p}(\mathit{x}|\mathit{A},\mathit{s})p(\mathit{s})\text{d}\mathit{s}\\ \propto \ln p(\widehat{\mathit{s}})-\frac{1}{2\sigma {}^2}|\mathit{x}-\mathit{A}\widehat{\mathit{s}}|{}^2-\frac{1}{2}\ln |\mathit{{\rm H}}|(23)\\ \mathit{{\rm H}}=\frac{1}{\sigma {}^2}\mathit{A}{}^{\text{Τ}}\mathit{A}-\frac{\text{d}}{\text{d}\mathit{s}\mathit{s}{}^{\text{Τ}}}\ln p(\mathit{s})|{}_{\mathit{s}=\widehat{\mathit{s}}}(24)\end{array}$

将式 (23) , 式 (24) 代入式 (14) , 然后通过梯度上升法寻找其最大值点, 从而得到矩阵A自然梯度迭代估计式如下:

$\text{Δ}\mathit{A}\propto \mathit{A}\mathit{A}{}^{\text{Τ}}\frac{\partial }{\partial \mathit{A}}\ln {\rm P}(\mathit{x}|\mathit{A})⧋-\mathit{A}[\phi (\widehat{\mathit{s}})\widehat{\mathit{s}}{}^{\text{Τ}}+\mathit{{\rm I}}](25)$

$\begin{array}{l}\text{式}\text{中}:\phi (\widehat{\mathit{s}})=[\phi (\widehat{s}{}_1),\phi (\widehat{s}{}_2),\cdots ,\phi (\widehat{s}{}_{{\rm N}})]{}^{\text{Τ}}‚\phi (\widehat{s}{}_i)=\\ \frac{\partial \log {\rm P}(\widehat{s}{}_i)}{\partial \widehat{s}{}_i}\text{称}\text{为}\text{激}\text{活}\text{函}\text{数}。\end{array}$

3.2 最大化后验概率法估计源信号[3,4,5,6]

在估计出混合矩阵后, 可以在假设源信号满足一定的概率分布下, 通过最大化信源的后验概率来实现对源信号的估计:

$\widehat{\mathit{s}}=\underset{\mathit{s}}{{\displaystyle \max }}{\rm P}(\mathit{s}|\mathit{x},\mathit{A})=\underset{\mathit{s}}{{\displaystyle \max }}{\rm P}(\mathit{x}|\mathit{A},\mathit{s}){\rm P}(\mathit{s})(26)$

采用对数表示后可得:

$\widehat{\mathit{s}}=\arg \underset{\mathit{s}}{{\displaystyle \max }}\left[\ln p(\mathit{x}|\mathit{A},\mathit{s})+\ln p(\mathit{s})\right](27)$

对式 (27) 的计算可以采用梯度上升法或拟牛顿法求极值, 两者相比, 拟牛顿法计算复杂度大幅度增加, 但却加快了收敛速度。

定义目标函数J=ln p (x|A, s) +ln p (s) , 则梯度法的源信号迭代过程为:

$\begin{array}{l}\mathit{s}(k+1)=\mathit{s}(k)+\eta {}_{\mathit{s}}\frac{\partial J}{\partial \mathit{s}}\\ =\mathit{s}(k)+\eta {}_{\mathit{s}}[-\mathit{A}{}^{\text{Τ}}\Psi (\widehat{\mathit{\epsilon }})+\phi (\widehat{\mathit{s}})](28)\end{array}$

$\begin{array}{l}\text{式}\text{中}:\widehat{\mathit{\epsilon }}=\mathit{x}-\mathit{A}\mathit{s}‚\Psi (\widehat{\mathit{\epsilon }})=[\frac{\partial \ln {\rm P}(\widehat{\epsilon }{}_1)}{\partial \widehat{\epsilon }{}_1},\frac{\partial \ln {\rm P}(\widehat{\epsilon }{}_2)}{\partial \widehat{\epsilon }{}_2},\cdots ,\\ \frac{\partial \ln {\rm P}(\widehat{\epsilon }{}_{{\rm M}})}{\partial \widehat{\epsilon }{}_{{\rm M}}}]{}^{\text{Τ}}。\end{array}$

拟牛顿法的迭代过程[6]为:

$\mathit{s}(k+1)=\mathit{s}(k)+\eta {}_{\mathit{s}}\left\{\frac{\partial {}^{2}J}{\partial \mathit{s}\mathit{s}{}^{\text{Τ}}}\right\}{}^{-1}\left[\frac{\partial J}{\partial \mathit{s}}\right](29)$

3.3 源信号模型

前面提到, 超完备条件下的盲分离, 由于混合过程中丢失的信息无法用传统的数学手段来恢复, 而只能通过一些先验或假设来弥补这些丢失的信息, 而在交替估计算法中, 无论是对混合矩阵的估计, 还是对源信号的估计, 都需要知道源信号的概率密度函数, 因此需要建立源信号模型, 一般可以采用拉普拉斯分布或高斯混合分布。

3.3.1 拉普拉斯分布模型[3,4,5]

对于语音信号, 大多具有超高斯分布, 所以可以假设源信号服从拉普拉斯分布:P (si) ∝exp (-α|si|) , 这时比较适合做语音信号的盲分离。通常还需假设源信号向量相互独立, 则:

${\rm P}(\mathit{s})={\displaystyle \prod _{i=1}^{{\rm N}}{\rm P}}(s{}_i)\propto {\displaystyle \prod _{i=1}^{{\rm N}}\exp }(-\alpha |s{}_i|)(30)$

此即为N维源信号的联合概率分布。

3.3.2 高斯混合分布模型[6]

对于很多通信信号或图像信号, 其概率密度并不满足拉普拉斯分布。于是提出了采用更宽泛的信号概率模型:高斯混合模型, 它既可以描述单峰分布, 又可以描述多峰分布;既可以描述亚高斯信号, 又可以描述超高斯信号。即每一个信号都可以用一个Q维的高斯信号按照一定的比例混合得到:

$\begin{array}{l}p(\widehat{\mathit{s}})={\displaystyle \sum _{q=1}^Q\mathit{\kappa }}{}_{q}p{}_q(\widehat{\mathit{s}}|\lambda {}_q)=\\ {\displaystyle \sum _{q=1}^Q\frac{\mathit{\kappa }{}_q}{(2\text{π}){}^{\frac{{\rm N}}{2}}|\mathit{\Sigma }{}_q|{}^{\frac{1}{2}}}}\exp \left[-\frac{1}{2}(\widehat{\mathit{s}}-\mathit{\mu }{}_q){}^{\text{Τ}}\mathit{\Sigma }{}_q^{-1}(\widehat{\mathit{s}}-\mathit{\mu }{}_q)\right](31)\end{array}$

式中:N是s的维数;κq, μq和Σq分别是混合权重, 均值和协方差矩阵, 它们均可以通过最大期望算法 (Expectation Maximization, EM) [6]或交替条件期望最大 (ACEM) 算法[13]进行估计, 在此不做具体介绍, 可参看相关文献。

3.3.3 信号模型比较

拉普拉斯分布模型, 参数较少, 因此计算代价相对较低, 但是一般只适用于超高斯稀疏分布的语音信号, 一旦源信号与假设不同, 将会造成较大的误差。

高斯混合模型能更好地模拟不同的源信号分布, 适用的范围比较宽, 但是参数很多, 且会随着源信号数目的增加呈指数增长, 加上EM算法本身运算复杂度高, 所以这种模型的运算量大, 计算代价很高。

4 结 语

本文总结了大部分的超完备盲分离算法, 将各算法归纳为两步估计法和交替估计法两类, 这两类算法各有许多不同的处理方法, 文中较详细地介绍了各算法的思想, 且对各算法的适用条件和优缺点进行了对比和总结。

从文中介绍的各算法来看, 目前超完备盲分离的研究基本上还都限于低维情况的研究, 各算法的复杂度都较高, 且都对系统模型加了许多假设和限制, 所以今后超完备盲分离的研究重点是在非稀疏情况下的分离, 在非对称信号的情况下的盲分离, 在高维观测数据情况下的盲分离, 以及计算复杂度的降低等方面的研究。

信噪比盲估计算法性能比较 篇4

文献[2][3]中使用极大似然估计法,将信号的信噪比估计转变为信号幅度值估计,文献[4][5]中利用接收信号的高阶矩估计信噪比。上述算法在各自特定的条件下非常有效,但大多在非协作通信的条件下不再适用,而且有些需要进行较大维数的矩阵运算,计算量大,难以满足实时性要求。

本文提出了一种基于功率谱差分的信噪比盲估计算法,并将其与基于奇异值分解的信噪比盲估计法在性能上作了详细比较,进行了仿真验证。可以发现,文中提出的算法具有明显的计算简单、估计精确等优点。

在高斯白噪声背景下,对于接收到的信号,可表示为:

其中,x(n)是信号分量,y(n)是观测信号,σ(n)为白噪声,N是采样得到的观测数据长度。

1 基于奇异值分解的信噪比盲估计

由(1)式,可以估计接收信号的m阶自协方差矩阵:

其中,

Cy是非负定矩阵,根据信号分量与噪声分量间的不相关性,可对其进行奇异值分解:

其中,U=[u1,u2,…,up,up+1,up+2,…,um]为Cy的特征空间,DP是由自协方差矩阵Cy的特征值构成的对角阵。

由上式可以看出,σ2为Cy的m-p重特征值,且是最小的m-p个特征值,将它们所对应的特征向量所张成的空间Un=[up+1,up+2,…,um]称为噪声子空间,将由剩余的特征值对应的特征向量所张成的空间Us=[u1,u2,…,up]称为信号子空间。此时接收信号的信噪比可以由特征值表示为:

但实际上,由于信号能量与噪声能量的交叉,使得部分信号能量分布在噪声子空间内,从而导致通过奇异值分解得到的噪声子空间的特征值并不是相等的m-p个值,此时由MDL准则[15],可估计信号子空间的维数,过程如下:

(1)首先对由式(4)得到的特征值进行排序,设λ1≥λ2≥…≥λm定义球形检验函数:

(2)定义MDL目标函数:

其中,K为利用信号遍历性计算信号期望值时的数据长度。信号子空间的维数可估计为:

(3)估计出信号子空间维数后,可估计噪声的功率:

(4)由式(5)修正后可求得信噪比的估计值:

2 基于功率谱差分的信噪比盲估计

在进行信号实时监测中,较大维数的矩阵运算给应用带来了困难,本章从频域出发,利用白噪声功率谱的平坦特征,提出一种基于功率谱分析的信噪比盲估计算法。

由式(1),可估计出接收信号的功率谱:

图1是某一含噪信号的功率谱,从图1可以看出,信号分量存在的频率段,功率谱幅度呈现明显的上凸,从几何学的角度,可以利用导数知识将功率谱图变化最快的地方提取出来,即近似对应着功率谱图中信号分量上凸段的端点。对图1中的功率谱进行前向差分:

其中,nfft是功率谱估计中FFT的长度。得到结果如图2所示。

由图2可以看出,在信号分量出现的频率段端点,功率谱差分图出现了明显的峰值,信号分量频段两端,分别对应着正峰和负峰,因此,实施功率谱差分后,可以设置门限,分别进行正负谱峰搜索,得到信号分量对应的频段估计。考虑到信号的拖尾性,为了保证取得信号分量能量的完整性,对于正谱峰,取其左侧第一个过零点为该信号分量频率段的左端点,对于负谱峰,取其右侧第一个过零点为该信号分量频率段的右端点。为了使功率谱差分后的峰值明显,要求估计出的功率谱尽可能平滑,因此,本文选择Welch法估计功率谱。

令估计出的信号分量频率段左右端点分别为ωL和ωH,则该信号分量的带宽可估计为:

其中,fs为采样频率。该信号分量的中心频率可估计为:

考虑到白噪声谱的平坦性,整个频段内的噪声功率可以利用非信号频率段的功率谱进行估计:

信号分量的功率可估计如下:

进而求得信噪比:

3 实验与仿真分析

实验信号为一段样本长度N为32768的QPSK复信号,采样频率为80MHz,其中QPSK信号的载波频率为20MHz,码速率为5MHz,成形滤波器为升余弦滤波,其滚降系数为0.5,信道为高斯白噪声信道,噪声功率为3dB,信噪比变化范围为-10~10dB,信号的协方差矩阵维数m取值为64(一般取值为50~100),计算期望值时的数据长度K取为N/2。利用Welch法估计功率谱时,每段数据长度L取为256,FFT长度与数据长度一致,段与段之间的重叠度取为50%。

3.1 算法复杂度比较

对于奇异值分解法估计信噪比,包括自协方差矩阵、奇异值分解、MDL法估计信号子空间维数和信噪比估计4个主要计算部分,如表1中所列。其中,决定计算复杂度的是奇异值分解,因此,基于奇异值分解的信噪比盲估计法的运算量为max,m是自协方差矩阵的维数,N是样本数据长度。

对于本文提出的基于功率谱差分算法,包括Welch法估计功率谱、功率谱差分、谱峰搜索和信噪比估计4个主要计算部分,如表2中所列。其中,决定计算复杂度的是利用Welch法估计功率谱及信噪比计算,运算量为o(N),N是样本数据长度。

由表1、表2可知,本文提出算法与基于奇异值分解的信噪比盲估计算法相比,具有更小的计算量,易于实现,能够更好地满足实时性要求。

3.2 算法性能比较

Mente Carlo仿真500次,得到两种算法在-10~10dB范围内的盲信噪比估计偏差与均方误差分别如图3、图4所示。

由图3可以看出,与奇异值分解法相比,本文提出算法在给定范围内的信噪比估计偏差一直较小,两者大约相差3dB。由图4可以看出,随着信噪比降低,两种算法估计出的信噪比均方误差均在逐渐增大,但本文算法估计出的信噪比误差波动小,在稳定度上较奇异值分解法有很大提高,算法的稳健性更强。

4 总结

本文介绍了基于奇异值分解的信噪比盲估计算法的基本原理,针对矩阵分解的计算复杂性,提出了一种基于功率谱差分的信噪比盲估计算法。研究表明,该算法与前者相比,不但具有较小的计算量,而且信噪比估计精度及稳定度均有较大提高。另外,该算法还适用于频域不重叠的单通道多信号的信噪比盲估计,具有较强的工程应用性。

摘要：为了实现非协作环境下的通信信号信噪比估计,本文在功率谱分析的基础上,提出了一种基于功率谱差分的信噪比盲估计算法,并将其与传统的基于奇异值分解的信噪比盲估计算法进行了比较。理论与仿真结果表明,与传统方法相比,本文提出的算法不但复杂度低,而且在低信噪比情形下仍具有较高的估计精确度和稳健性。

关键词：信噪比盲估计,异值分解,功率谱差分,Welch法

参考文献

[1]唱亮,汪芙平.非协作通信中的盲信噪比估计算法.通信学报,2008,29(3):76-81

[2]隋丹,葛临东.一种新的基于改进PASTd的中频信号盲信噪比估计算法.电子与信息学报,2007,29(7):1657~1661

[3]孙钢灿,安建平,杨杰,杨凯.非协作通信中的信噪比估计算法.北京理工大学学报,2009,29(8):708-712

[4]彭耿,黄知涛,陆凤波,姜文利.中频通信信号信噪比的快速盲估计.电子与信息学报,2010,32(1):102-106

粒子群算法用于盲信号分离的研究 篇5

盲信号分离 (BSS) 是20世纪末发展起来的新研究领域, 它指在对源信号、传输通道都未知的情况下, 仅从若干个观测到的混合信号中恢复出源始信号[1]。目前, 对BSS问题存在多种算法。随机梯度、自然梯度等梯度算法是一种较常见的算法, 然而该方法的关键是需要计算非线性激活函数, 使算法的复杂度增加, 同时算法需要计算矩阵的产逆, 计算量大, 而且鲁棒性较差, 收敛性有待改进[2]。文献[3,4]提出基于神经网络的盲分离算法, 该方法收敛速度较慢, 且容易陷入局部极值。本文将粒子群算法引入BSS中, 为BSS领域提供一种新的方法, 该方法简单高效, 与传统方法复杂、收敛慢的特点相比, 有较大的优势。

1 盲分离问题描述

设有n个混合信号X (t) 是n个未知的相互统计独立的源信号S (t) 的线性组合, 即:

$\mathit{X}(t)=\mathit{{\rm H}}\mathit{S}(t)(1)$

式中:X (t) =[x1 (t) , x2 (t) , …, xn (t) ]T;S (t) =[s1 (t) , s2 (t) , …, sn (t) ]T;H为n×n阶混合未知矩阵, H=[H1, H2, …, Hn], 其中Hi=[hi1, hi2, …, hin], i=1, 2, …, n。盲分离就是要在H和S未知的情况下, 使用一定的算法, 对观测到的X (t) 进行处理, 得到一个n×n阶分离矩阵ω, 使得到的输出信号:

$\mathit{Y}=\mathit{\omega }\mathit{X}(t)(2)$

是源信号X (t) 的最好估计。

2 基于粒子群算法的盲分离

2.1 粒子群算法

粒子群算法 (PSO) 是一种智能优化算法[5,6]。在PSO中, 有N个微粒和一个适应度函数, 每个微粒在D维搜索空间中以一定的速度飞行。该飞行速度由个体的飞行经验动态调整。各微粒根据自己的当前位置和速度不断更新自己的位置。微粒i的位置和飞行速度分别表示为[Xi1, Xi2, …XiD]和[Vi1, Vi2, …, ViD]。每个微粒存在一个由它的位置决定的与适应度函数相对应的值, 根据微粒的适应度值判断粒子的优劣, 取最优粒子的位置作为最优解。粒子的飞行速度和位置调整的方程如下:

$\begin{array}{l}V{}_{ij}(t+1)=wV{}_{ij}(t)+c{}_1(t)[p{}_{ij}(t)-x{}_{ij}(t)]+\\ c{}_{2}r{}_2(t)[p{}_{gj}(t)-X{}_{ij}(t)]\\ X{}_{ij}(t+1)=X{}_{ij}(t)+V{}_{ij}(t+1)(3)\end{array}$

式中:t为当前迭代次数;i表示第i个微粒;j表示第j维;c1, c2为区间[0, 2]的常数, 其中c1为调节微粒飞向自身最好位置 (pi) 方向的步长, 为认知系数;c2为调节微粒飞向全局最好位置 (pg) 方向的步长, 为社会学习系数;r1～U (0, 1) , r2～U (0, 1) 为两个独立分布的随机数;w为惯性权重, 让w随着迭代次数增加而减少, 可使PSO算法在初期有较强的全局收敛能力, 而晚期有较强的局部收敛能力[7,8]。

2.2 分离原理

为了对混合信号X (t) 进行盲分离, 本文提出使用PSO算法求解式 (2) 中的分离矩阵ω。根据中心极限定理可知, 由多个独立随机变量之和的分布比其中任一变量更接近高斯分布, 因此可将信号的非高斯性作为信号独立性的度量[9,10]。式 (2) 的目标就是要找出分离矩阵ω, 使输出信号Y=ωX (t) 的非高斯性最强, 这时便由混合信号分离出了源信号。这时有ω=H-1, 于是有:

$\mathit{Y}=\mathit{\omega }\mathit{X}(t)=\mathit{{\rm H}}{}^{-1}\mathit{{\rm H}}\mathit{S}(t)=\mathit{S}(t)$

从而实现混合信号的盲分离。

峰度是随机信号非高斯性的一种判据。随机信号y的峰度定义为:

$\text{k}\text{u}\text{r}\text{t}(y)=E\{y{}^4\}-3E\{y{}^2\}{}^2(4)$

当kurt (y) =0时, y为高斯信号;|kurt (y) |越大, y的非高斯性越强。因此, 本文使用PSO算法, 采用峰度的绝对值如式 (5) 作为适应度函数进行寻优:

$\begin{array}{l}F(\mathit{\omega })=\left|\text{k}\text{u}\text{r}\text{t}(Y)\right|=\left|{\displaystyle \sum _{i=1}^n\text{k}}\text{u}\text{r}\text{t}(Y{}_i)\right|=\left|{\displaystyle \sum _{i=1}^n\text{k}}\text{u}\text{r}\text{t}(\omega {}_{i}X)\right|\\ =|{\displaystyle \sum _{i=1}^{n}E}((\omega {}_{i}X){}^4)-3[E((\omega {}_{i}X){}^2){}^2|(5)\end{array}$

2.3 算法步骤

采用PSO算法对混合信号X (t) 进行盲分离的具体步骤如下:

(1) 确定算法的各个参数值:c1, c2, vmax, MaxDT并初始化各微粒。

(2) 根据式 (5) 计算各微粒的适应度值。

(3) 对每个微粒, 判断是否更新粒子的局部极值和全局极值。

(4) 按进化方程式 (3) 进化各微粒。

(5) 判断是否满足终止条件, 否就转步骤 (3) , 是就转步骤 (6) 。

(6) 由pg得到分离矩阵ω, 从而得到分离信号Y。

3 仿真分析

用Matlab对本算法进行实例仿真分析。选取2个图像作为源信号, 如图1所示, 以线性方式混合, 混合矩阵, 混合图像如图2所示。算法参数设置为:群体大小N=30, 参数维数dim=4, 惯性权重w由1.2线性变到0.9, 粒子群中所有粒子位置的取值范围都限制在[-1, +1], 最大速度Vmax=1, 最大迭代次数MaxDT=100。仿真得到分离出来的两个图像如图3所示, 可以看到分离图像与源图像非常接近, 只是顺序不一样, 这是由盲分离问题的不确定性引起的, 不影响结果的正确性。因此本算法成功地分离出了两个图像。

为了评价本算法的优劣, 将混合图像用自然梯度算法进行盲分离, 将两种算法进行对比。表1是两种算法进行盲分离的收敛情况表。可以看到, 收敛的时候, PSO算法所用的迭代次数要少于自然梯度算法, 因此PSO算法的收敛速度快于自然梯度法的收敛速度。比较两种算法的抽取精度, 可以用信噪比来表示。将分离信号与源信号之差看作是噪声, 分离出的第i个信号的信噪比 (单位:dB) 可以表示为:

$\text{S}\text{Ν}\text{R}{}_i=10\log {}_{2}10\{{\displaystyle \sum _{k=1}^{{\rm N}}s}{}_i^2(t)/{\displaystyle \sum _{k=1}^{{\rm N}}[}y{}_i(t)-s{}_i(t)]{}^2\}(6)$

式中:N为信号的取样点数。

表2是两种算法的信噪比比较表, 可以看到PSO算法比自然梯度算法分离的图像有较高的信噪比, 因此本算法有较高的准确度。

4 结 语

将粒子群算法引入盲信号分离中, 并通过仿真实例成功分离出源信号, 与自然梯度盲分离算法对比, 粒子群盲分离算法在分离精度和收敛速度上有一定的优越性, 精度更高, 收敛速度更快, 为盲信号分离领域提供了一种新的研究思路与方法。

参考文献

[1]张贤达, 保铮.盲信号分离[J].电子学报, 2001, 29 (12A) :1766-1771.

[2]CICHOCKI A, AMARI S, CAOJ.Neural network modelsfor blind separation of ti me delayed and convolved signals[J].Japanese IEICE Trans.on Fundamentals, 1997, E-82-A (9) :1595-1600.

[3]COMON P, JUTTEN C, HERANLT J, Blind separationof sources:problemstatement[J].Signal Processing, 1991, 24 (1) :11-20.

[4]HERAULT J, Jutten C.Blind separation of sources:anadaptive algorithm based on neurominetic[J].Signal Pro-cessing, 1991, 24 (1) :1-10.

[5]谢晓峰, 张文俊, 杨之廉.微粒群算法综述[J].控制与决策, 2003, 18 (2) :129-134.

[6]KENNEDYJ, EBERHART R C.Particle swarmopti miza-tion[C]//IEEE International Conference on Neural Net-works.Perth, Australia:IEEE 1995:1942-1948.

[7]SHI Yu-hui, EBERHART R C.Empirical study of particleswarmopti mization[C]//Proceeding of the IEEE Congresson evolutionary Computation.Piscataway, NJ:IEEEPress, 1998:69-73.

[8]曾建潮, 介婧, 崔志华.微粒群算法[M].北京:科学出版社, 2004.

[9]张贤达, 牛奕龙, 陈海洋.盲信号处理[M].北京:国防工业出版社, 2006.

一种雷达信号自适应盲分离算法 篇6

雷达系统在工作状态中的接收信号包括回波信号、干扰信号、杂波信号、内部噪声等,雷达工作环境的恶化和日益复杂的电子战信号环境,使现代雷达系统面临严峻挑战,使得各式雷达都必须具备非常强的干扰环境中检测目标和提取目标参数的能力[1]。对于较简单的干扰模式,现有的雷达信号处理模块通过经典的滤波理论和传统的信号处理方法加以处理,就可以获取较佳的效果。随着干扰信号模式愈来愈复杂多变,现代雷达侦察设备的信号分选任务也十分艰巨而复杂,传统方法已无法胜任这方面的工作。因为传统的“主动雷达”为了探测目标必须发出电波信号,很容易因暴露自身而遭到攻击,所以使用受到很大限制。相反,近年来发展起来的“被动雷达”(暂且称其为雷达),由于只接收信号而不发出任何信号就可以探测到目标,因而受到各国的广泛重视,实际上,这种被动雷达工作的基本原理就是盲源分离技术。我国这几年来对隐形飞机探测研究所取得的成就也正是将盲源分离领域的最新研究成果应用于上述这类被动雷达的一个具体生动的例子。

盲分离技术是近年来信号处理技术的重要发展方向,是指信号分离时,我们无法预知原始信号和传输信道的基本信息,仅根据接收到的混合信号有效地将图像或信息分离[2]。利用盲信号分离技术能很好地解决复杂环境背景下雷达信号分选的问题,它无需考虑过多的雷达信号的环境、条件的可测性。学习样本的选取,只需根据接收设备所获取的雷达辐射源信号进行处理,就可以得到原始信号的形式,即恢复出原始信号,为电子战中对抗和反对抗采取措施提供了重要的依据。

1 雷达信号盲分离的数学模型

雷达对抗信号环境S(t)是指雷达对抗设备在其所在的地域内存在的各种雷达辐射、散射信号的整体,数学表达式为:

雷达信号的干扰信号包括噪声干扰信号和人为干扰信号。噪声干扰信号一般有噪声调幅干扰信号、压制性干扰信号、噪声调幅—调频干扰信号。人为干扰信号形式繁多,本文以典型的高斯噪声和线性调频信号作为干扰信号为例。系统有M路接收通道,相互等距为d,有N路窄带信号,气象雷达接收信号模型可表示为:

其中,X(t)为雷达接收到的观测信号,A=M×N为混叠矩阵,它为天线阵的响应函数和转播过程中的混叠矩阵的乘积,S(t)为M×1维源信号,包括回波信号和干扰信号等,假设各信号之间与噪声之间均是相互独立;n(t)为N×1维噪声信号,包括外部噪声,内部噪声和电噪声等,通常视为高斯白噪声。

雷达接收机处理数据是通过回波信号和本振信号进行混频,并对混频信号进行调节和跟踪,在复杂的信号背景下,利用自适应神经网络分离系统区分干扰和目标信号,关键在于如何求解一个好的分离系统,这里我们采用两步法———白化矩阵U和分离矩阵W来完成分离的目的。

因此,解混矩阵的数学表达式为:

其中,y(t)=(y1(t),y2(t),L,yn(t))表示恢复出来的信号,是对原始信号S(t)的估计。盲分离的过程,就是根据对源信号S(t)的性质和特征的假设,建立解混矩阵B的目标函数,使得通过解混矩阵B的输出信号y(t)尽可能地逼近S(t)。由(2)、(3)式得:

G称为全局矩阵,最理想的分离效果是G=I(I为n×n阶单位矩阵),即B=A-1。所以盲分离算法的核心问题就是如何寻找到分离矩阵的最佳估计值,它包括两个方面:优化判据(或比值函数)和寻优算法,即学习时首先建立一个以已知信息元素为变元的目标函数ρ,其次是寻找一种有效的算法求解W,若能够算出某个能使ρ达到极大(小)的值,该极值点即为W所需的解。

2 雷达信号的自适应不完整自然梯度盲分离算法

如何求得分离矩阵W,是ICA的核心问题。通常这个问题需要基于目标函数的方法进行求解,ICA问题的解B将在这些目标函数的极小或极大值处找到。因此,对于极小化一个多元函数的无约束问题,最经典的方法就是最速下降法,也称为梯度下降法。本文采用不完整自然梯度法,具体如公式:

不完整自然梯度学习算法的特点是当源信号幅度随时间快速变化或在某段时间为零时,仍能很好地工作,甚至在过高估计源信号数目时,算法性能也不会出现较大波动。而标准自然梯度算法则不具备这种能力,它们会放大不包含在源信号内的微小分量,造成算法性能的恶化,这点Amari等人在文献[3,4]中加以论证,而在具体的仿真实验中,不完整自然梯度算法的分离效果和稳定性也强于自然梯度法,因此本文将采取不完整自然梯度法进行图像分离。公式中步长因子μ的作用是控制分离矩阵W在每次迭代过程中的更新幅度,本文中主要验证自适应算法的可行性,故采用固定步长。

最终得到本文不完整自然梯度变步长算法的具体实现步骤如下:

(1)对混合数据X进行中心化,使其均值为零;

(2)白化:估计白化阵U,获得白化数据Z=UX;

(3)设置初始分离矩阵W。根据需要处理的信号选择合适的φ(·)函数。选择初始步长μ(k)

(4)计算Y=WZ;

(5)计算F(Y)=φ(Y)YT,从F(Y)中获得Λ=diag{λ1,L,λn},其中元素λi=F(Y)中正对角线元素;

(6)更新分离矩阵

(7)如果尚未收敛,返回步骤4,继续迭代直至收敛。

3 仿真实现

假设位于接收机阵列远场有4个雷达信号源,且各雷达信号源之间是相互独立的,噪声为零均值高斯噪声,各接收机噪声之间是相互独立的,信号源与噪声之间也是相互独立的,假设各信号源为4个正弦波调制信号,它们的包络和载频各不相同,运用本文的分选算法可以得到无噪声的仿真结果如图1所示。

4 结束语

一直以来,独立分量分析作为成熟的盲源分离算法,在理论和实际应用上都取得了很大的发展,被广泛地应用于各个学科领域,出现了许多热门方向。而自适应神经网络盲源分离技术是近几年来发展起来的一门新型技术学科,以其无需较多先验知识、快速有效的分离特点,得到更多科研人员的关注,本文把盲源分离技术应用到雷达信号的分离上,并通过算法分析和仿真实验证明了算法的可行性,为雷达技术的发展方向提供了较好的参考价值。

参考文献

[1]杨建民,黄晓刚,徐佳龙,等.基于信号处理优化设计的雷达性能改进[J].现代雷达,2010,3(32):21-23.

[2]杨福生,洪波.独立分量分析的原理与应用[M].北京:清华大学出版社,2006.

[3]Amari S,Cichocki A,Yang HH.A new learning algorithm for blind signal separation[J].Advances in neural information processing systems,1996:757-763.

[4]Amari S,Chen T P,Cichocki A.Non-holonomic Orthogonal Learning Algorithms for Blind Source Separation[J].Neural Computation,2000,6(12):1463-1484.

基于均匀线阵的盲分离算法 篇7

盲信号分离是近年来阵列信号处理领域的研究热点,广泛应用于语音分离与识别、雷达信号处理和声纳信号处理、电子对抗等。基于多数通信信号具有的恒模特性发展起来的恒模盲分离算法[1,2,3]是一种比较著名的迭代型算法。恒模算法(Constant Modules Algorithm,CMA)首先由D.Godard 提出[4],随后Agee于1986 年提出了一种具有快速收敛特性的算法一最小二乘恒模算法(Least Square Constant Modules Algorithm,LSCMA)[5,6],取得了很好的效果。恒模算法利用发送信号的幅度统计特性来调整权系数,使得输出信号的幅度保持恒定,具有计算复杂度低、收敛速度快、不受收敛步长的影响和对阵列模型的偏差不敏感等显著的优点,能够方便地用在通信环境中。利用最小二乘恒模盲分离算法并结合均匀线阵的特点对混合矩阵进行修正,可以达到改善算法分离性能的目的。

1 接收阵列模型

如图1所示,考虑均匀线阵接收系统,阵列天线的个数为N,以阵元1为参考阵元,各阵元间距为d,为入射信号最高频率对应波长的一半(为避免相位模糊)。空间有M个中心频率为f1,f2,…,fM的窄带信号源,入射角分别为θ1,θ2,…,θM,则阵列输出信号为:

X(t)=AS(t)+N(t)。 (1)

式中,阵列输出信号X(t)=[x1(t),x2(t),…xN(t)]T;源信号S(t)=[s1(t),s2(t),…sM(t)]T;N(t)为服从正态分布的复高斯白噪声。

混合矩阵A为:

式中,ωm=2πfm,为第m个信号的角频率;τm=dsinθm/c,为该信号到达相邻天线阵元的时延差,θm为相应的入射角,c为光速。

2 最小二乘恒模盲分离算法

2.1 算法分析

最小二乘恒模算法利用扩展的高斯方法最小化恒模代价函数,扩展的高斯方法定义的代价函数F(w)为:

$F(\mathit{w})={\displaystyle \sum _{l=1}^n\text{ϕ}}{}_l^2=\parallel \text{ϕ}(\mathit{w})\parallel {}_2^2$。 (3)

由泰勒级数展开知:

$F\left(\mathit{w}+\mathit{r}\right)\approx \parallel \text{ϕ}\left(\mathit{w}\right)+J{}^{\text{Η}}\left(\mathit{w}\right)\mathit{r}\parallel {}_2^2$。 (4)

式中,r为偏差向量;J(w)为ϕ(w)的Jacobian行列式,表示为:

$J(\mathit{w})=\left[\nabla \text{ϕ}{}_1(\mathit{w}),\nabla \text{ϕ}{}_2(\mathit{w}),\cdots ,\nabla \text{ϕ}{}_n(\mathit{w})\right]$。 (5)

关于r对F(w+r)求梯度得:

$\nabla {}_r(F(\mathit{w}+\mathit{r}))=2\left[J(\mathit{w})\text{ϕ}(\mathit{w})+J(\mathit{w})J{}^{\text{Η}}(\mathit{w})\mathit{r}\right]$。 (6)

令ᐁr(F(w+r))=0,求出使F(w+r)最小的偏差向量为:

$\mathit{r}=-\left[J(\mathit{w})J{}^{\text{Η}}(\mathit{w})\right]{}^{-1}J(\mathit{w})\text{ϕ}(\mathit{w})$。 (7)

将偏差向量与权向量w(k)相加,可以得到权向量的更新公式为:

w(k+1)=w(k)-r(k)。 (8)

将上述方法应用于最小二乘恒模盲分离算法,令代价函数为:

$F(\mathit{w})={\displaystyle \sum _{k=1}^n(\parallel \mathit{y}(k)\parallel {}^l-1)}{}^2={\displaystyle \sum _{k=1}^n(\parallel \mathit{w}{}^{\text{Η}}\mathit{x}(k)\parallel {}^l-1)}{}^2$。 (9)

此时,

ϕ(w)=‖wHx‖-1。 (10)

与上述推导相似,可以得到权向量更新为:

w(k+1)=(XXH)-1Xr*。 (11)

式中,

$\mathit{X}=\left[\mathit{x}(1),\mathit{x}(2),\cdots ,\mathit{x}({\rm K})\right]$;

$\begin{array}{l}\mathit{y}(k)=\left[\mathit{w}{}^{\text{Η}}(k)\mathit{X}\right]{}^{\text{Τ}}‚\\ \mathit{r}(k)=L(y(k))=\left[\frac{y(1)}{\left|y(1)\right|},\frac{y(2)}{\left|y(2)\right|},\cdots ,\frac{y({\rm K})}{\left|y({\rm K})\right|}\right]{}^{\text{Τ}}。\end{array}$

显然,实施式(11)的算法一次可以得到其中的一个源信号,为了得到所有的源信号,需要运行以上算法M次,为了避免重复分离同一个信号,需要对分离向量做式(12)的正交化处理:

$\mathit{w}{}_p\leftarrow \mathit{w}{}_p-{\displaystyle \sum _{j=1}^{p-1}(}\mathit{w}{}_p^{\text{Τ}}\mathit{w}{}_j)\mathit{w}{}_j$。 (12)

2.2 算法评价

分离出来的信号其实是源信号S(t)的一个估计,可将其记为$\mathit{Y}=\widehat{\mathit{S}}(\mathit{t})$,于是分离信号与源信号之间的相关系数为:

$\zeta (y{}_i,s{}_j)=\left|{\displaystyle \sum _{t=1}^{{\rm M}}y}{}_i(t)s{}_j(t)\right|$$\sqrt{{\displaystyle \sum _{t=1}^{{\rm M}}y}{}_i{}^2(t){\displaystyle \sum _{t=1}^{{\rm M}}s}{}_j{}^2(t)}$。

可用于衡量源信号sj(t)和分离信号yj(t)的相似程度,si(t)和yj(t)相关性越强,ζij越大,当si(t)=cyj(t)(允许盲分离效果在幅度上存在差异)时,ζij=1;当yi与sj相互独立时,ζij=0。由相关系数即可构成相关矩阵,如果分离效果良好,则分离信号与对应的源信号之间的相关系数应该接近1或-1;反过来,如果分离效果很差,则相关系数应接近0。即理想情况下的相关系数矩阵每行每列都有且仅有一个元素近似接近于l,其他元素都近似为0。

3 修正算法

上述最小二乘恒模盲分离算法只利用了通信信号的恒模特性,下面结合均匀线阵的特点对混合矩阵进行重构,修正原算法,改善算法性能。

由式(2)可知,均匀线阵中混合矩阵A的任意相邻2行的比值向量为:

$\frac{\mathit{A}(n+1,:)}{\mathit{A}(n,:)}=[\text{e}{}^{\text{j}\omega {}_1\text{Δ}t{}_1},\text{e}{}^{\text{j}\omega {}_2\text{Δ}t{}_2},\cdots \text{e}{}^{\text{j}\omega {}_{{\rm M}}\text{Δ}t{}_{{\rm M}}}]$。 (13)

式中,A(n,:)表示矩阵A中第n行的所有元素;Δtm=dsinθm/c。

可以看出,混合矩阵中任意相邻2行的比值都相同。根据该特点,先利用已有的最小二乘盲分离算法对白化后的观测信号进行分离,得到分离矩阵W,则总的分离矩阵W1=WQ(Q为白化矩阵),对W1求逆(当N=M)或伪逆(当N>M)就能得到混合矩阵A的估计,设为A1,即A=(W1)-1。在理想情况下,A1也应该具有式(2)所示的结构。实际情况下,由于受噪声影响,A1相邻2行的比值可能不完全相等,因此对A1中第m列m=0,1,…M-1所有相邻2行的比值取平均得到am作为ejωmΔtm的估计:

$a{}_m=\text{e}{}^{\text{j}\omega {}_m\text{Δ}t{}_m}=\text{e}{}^{\text{j}2\text{π}f{}_m\frac{d\text{s}\text{i}\text{n}\theta {}_m}{\text{c}}}$。 (14)

由式(14)估计出来的am值,就可以利用混合矩阵的特点重构出混合矩阵,用A2表示为:

$\begin{array}{l}\mathit{A}{}_2=\left[\begin{array}{ccc}1& \cdots & 1\\ \text{e}{}^{\text{j}\omega {}_1\frac{d\text{s}\text{i}\text{n}\theta {}_1}{\text{c}}}& \cdots & \text{e}{}^{\text{j}\omega {}_{{\rm M}}\frac{d\text{s}\text{i}\text{n}\theta {}_{{\rm M}}}{\text{c}}}\\ ⋮& \ddots & \\ \text{e}{}^{\text{j}\omega {}_1({\rm N}-1)\frac{d\text{s}\text{i}\text{n}\theta {}_1}{\text{c}}}& \cdots & \text{e}{}^{\text{j}\omega {}_{{\rm M}}({\rm N}-1)\frac{d\text{s}\text{i}\text{n}\theta {}_{{\rm M}}}{\text{c}}}\end{array}\right]=\\ \left[\begin{array}{cccc}1& 1& \cdots & 1\\ \widehat{a}{}_1& \widehat{a}{}_2& \cdots & \widehat{a}{}_{{\rm M}}\\ ⋮& ⋮& \ddots & ⋮\\ \widehat{a}{}_1^{({\rm N}-1)}& \widehat{a}{}_2^{({\rm N}-1)}& \cdots & \widehat{a}{}_{{\rm M}}^{({\rm N}-1)}\end{array}\right]。(15)\end{array}$

根据A2,可以估计出修正后的分离矩阵W2为:

W2=(A2)-1。 (16)

则分离信号(源信号的修正估计)为:

y2=W2x。 (17)

综上所述,基于均匀线阵的修正最小二乘恒模盲分离算法步骤如下:

① 对观察信号x进行预白化处理,得到z,白化矩阵为Q;

② 设置W的初始值,用最小二乘恒模盲分离算法对分离矩阵W进行迭代和正交归一化,直到收敛。得到分离矩阵W以及分离信号y=Wz=W×Q×x,进而得到总的分离矩阵W1=W×Q;

③ 对矩阵W1求逆(矩阵为方阵时)或伪逆(矩阵不是方阵时)得到混合矩阵的估计A1=(W1)-1;

④ 求A1中第m列的所有相邻2行比值的平均值am(m=1,2,…,M);

⑤ 利用am,根据式(15)重构得修正算法对应的混合矩阵A2;

⑥ 对A2求逆(或伪逆)得到其相应的分离矩阵W2=(A2)-1,故由修正算法得到的分离信号为y2=W2x。

4 基于盲分离的高分辨测向

利用盲分离算法可以得到所有的源信号,进而可以得到混合矩阵的估计$\widehat{\mathit{A}}$:

$\widehat{\mathit{A}}=\underset{A}{{\displaystyle \text{a}\text{r}\text{g}\text{m}\text{i}\text{n}}}\left\{(\mathit{X}-\mathit{A}\widehat{\mathit{S}})(\mathit{X}-\mathit{A}\widehat{\mathit{S}}){}^{\text{Η}}\right\}$。 (18)

式中,$\widehat{\mathit{S}}$为分离得到的源信号的估计,式(18)是一个关于A的二次型方程,A的最优解一定存在且唯一,有

$\widehat{\mathit{A}}=\mathit{X}\widehat{\mathit{S}}{}^{\text{Η}}(\widehat{\mathit{S}}\widehat{\mathit{S}}{}^{\text{Η}}){}^{-1}$。 (19)

利用$\widehat{\mathit{A}}$的列向量$\left\{\widehat{\mathit{a}}(\theta {}_i)\right\}{}_{i=1}^{{\rm M}}$在真实阵列流形上的投影,就可以得到源信号的方向角估计为:

$\theta {}_i=\arg \max \left(\frac{1}{\mathit{a}(\theta {}_i){}^{\text{Η}}(\mathit{{\rm I}}-\widehat{\mathit{A}}\widehat{\mathit{A}}{}^{\text{Η}})\mathit{a}(\theta {}_i)}\right)$。 (20)

5 仿真实验

仿真中接收阵列为9元均匀线列阵,阵元间距为0.5倍信号波长。信号源1为FM信号,入射角为10°信号源2为BPSK信号,入射角为-10°,2个信号源的带宽为2 MHz,信噪比为10 dB,信号采样频率为40 MHz,样本点数为2 048。仿真在中频进行,中频频率为12 MHz。

原算法和修正算法分离源信号对应分离权得到的波束图如图2和图3所示,利用盲分离测向图和MUSIC测向对比图如图4所示。两算法都可实现混合信号的分离,在分离信号的到达角附近都形成高增益主瓣,而在另一信号到达角附近形成零限,修正算法有更深的零限,说明修正算法大大减弱了分离信号间的相关性。

分离信号与源信号的相关系数矩阵同样证明了以上结论,仿真得到原算法分离信号与源信号的相关系数矩阵为ζ1

$=\left[\begin{array}{cc}0.8922& 0.2321\\ 0.2359& 0.9622\end{array}\right]$

,修正算法分离信号与源信号的相关系数矩阵为ζ2=

$\left[\begin{array}{cc}0.9992& 0.0031\\ 0.0154& 1.0000\end{array}\right]$

,修正算法对应的相关系数矩阵更接近单位阵,分离信号与相应源信号相关性增强,分离信号之间的相关性近似为零,有效地降低了分离信号间的相关性。

修正恒模盲分离算法和经典MUSIC算法有相似的测向精度。通过100次独立仿真实验,基于修正盲分离算法测向的平均测向精度为0.211°,证明该算法可以实现对空间信号源的高精度测向。

6 结束语

介绍了一种基于均匀线阵的修正恒模盲信号分离算法,有效地降低了分离信号之间的相关性,提高了分离精度,并且将分离结果应用于测向实现了对空间信号源的高分辨测向。在测向侦察一体化中,利用该算法既可有效恢复辐射源信号,也可以对辐射源进行测向,对于提高系统对复杂信号环境的适应能力有一定的研究意义。 

参考文献

[1]张贤达,保铮.通信信号处理[M].北京:国防工业出版社,2000:390-403.

[2]CARDOSO J F.Equivariant Adaptive Source Separation.IEEE Trans.on Signal Processing[J].1996,44(12):3 017-3 029.

[3]XU Chang-jiang,LI Jian.A Batch Processing ConstantModulus Algorithm[J].IEEE Com.Letters,2004,8(9):582-584.

[4]GODARD D.Self-Recovering Equalization and CarrierTracking in Two-Dimensional Data Commu-nicationSystems[J].IEEE Trans on Commun.,1980,28(11):1 867-1 875.

[5]AGEE B G.The Least Squares CMA:A New Technique forRapid Correction of Constant Signals[C].USA:ProcICASSP,1986:953-956.

盲嵌入算法 篇8

1 数字水印技术

图像水印技术是采用一定的算法将一些标志性的信息嵌入到图像中,而不影响原图像的画质和使用。按水印嵌入时处理方法的不同,可将目前的图像水印处理算法分成空域水印处理算法和变换域水印处理算法两大类[1,2]。空域处理算法是将水印信息以比特流的形式嵌入到原始图像各像素点最不重要位上,来保证嵌入的水印不可见。但是由于使用的是图像不重要的像素位,因此对于滤波、平滑、压缩等图像处理,算法的鲁棒性较差[3]。变换域水印算法是先将原始图像通过某种变换(如离散余弦变换、小波变换等),转换到变换域,然后再在变换域中嵌入水印,最后经反变换输出含水印的空域图像。这类算法的水印嵌入容量相对较小,但可充分利用人类的感知特性实现良好的不可见性和鲁棒性。

2 DCT域数字水印技术的实现

DCT变换域数字水印是目前数字水印中较为常用的一种稳健的算法。在DCT域中,因为直流分量和中、低频分量是图像信号的主要成分,携带较多的信号能量,因此在此区域嵌入的数字水印具有较高的鲁棒性[4~5]。但改变DC系数容易导致块效应,所以算法中选取AC系数中的中、低频分量来嵌入水印信息。

2.1 水印的嵌入

在水印的嵌入过程中首先采用Arnold方法对图像水印进行置乱处理,处理之后的图像水印变得杂乱无章,打乱了像素空间位置间的关联关系,使其能均匀地分布在载体图像的所有空间,一方面提高了算法的鲁棒性,另一方面也可以加强水印的保密性。具体的水印嵌入处理过程如下:

(1)利用Arnold方法将大小为L×L的二值水印图像W进行n次变换。采用的Arnold置乱公式为:

x、y为图像W像素点的坐标,x'、y'为对应像素点变换后的坐标。

(2)将K×K(K=L×8)原始灰度图像分成8×8大小的图像块并对每个图像块进行DCT变换。

(3)选取各块DCT系数(m,n)位置处的数据Block-dcti(m,n)(i=1,2,…,L×L)嵌入水印信息。仿真中(m,n)=(1,3)。嵌入水印的公式为:

ai为二值水印图像按行扫描排列后的第i个像素值。

(4)对各DCT块进行DCT反变换,得到嵌入水印后的图像。

2.2 水印的提取

水印的提取方法采用盲检形式,不需要原始图像参与,具体过程如下:

(1)对已嵌入水印的图像进行8×8分块,各块再进行DCT变换;

(2)判断各块(m,n)位置处的DCT系数是否大于0。若大于0,则令bi=1,否则bi=0;

(3)将一维序列B={bi}重组为二维二值图像Wate Mark;

(4)对Wate Mark再次进行24-n次Arnold变换,生成提取出的水印图像。

3 仿真实验

通过Matlab6.5编程对算法进行仿真,采用256×256的256级灰度图像作为原始图像Image,以32×32的二值图像作为水印W',并对W进行4次Arnold迭代形成置乱后的水印W',分别如图1中(a)、(b)、(c)所示。在图1(a)中嵌入图1(c)水印之后的图像如图2(a)所示,以及从中提取的图像水印如图2(b)所示。由图可知,人眼无法感知是否嵌入了水印,算法有很好的不可见性。

图2(a)中嵌入水印后图像的峰值信噪比PSNR=35.1108d B,图像的改变在视觉上不明显,满足不可感知的要求。在系统没有受到攻击的情况下,图1(b)与图2(b)归一化相似度NC=1。

结束语

提出并实现了一种直接修改DCT域中、低频AC系数实现二值图像水印嵌入和提取的方法。采用该方法嵌入的水印图像具有较好的透明性,并且在提取时可实现盲检测。通过攻击实验表明:该方法对JPEG压缩、滤波、噪声、剪切等攻击具有较好的鲁棒性。

摘要：介绍了数字水印技术的一般处理方法,给出了一种DCT域盲水印的嵌入和提取的具体操作处理过程,并以Matlab为开发工具,对DCT域图像水印的嵌入与盲检测进行仿真。实验结果表明,该算法对一些常见的攻击,如JPEG压缩、滤波、噪声、剪切具有一定的鲁棒性。

关键词：数字水印,DCT变换,盲检测

参考文献

[1]郎锐.数字图像处理学Visual C++实现[M].北京:北京希望电子出版社,2003.

[2]肖力.一种基于DCT域的数字水印的实现[J].鄂州大学学报,2005,12(6):32-34.

[3]杨文泉,陆阳.基于离散余弦变换图象水印算法的研究[J].计算机工程与应用,2003(13):101-105.

[4]黄继武,SHI Y Q,程卫东.DCT域图象水印:嵌入对策和算法[J].电子学报,2000,28(4):386-38.