简介一份基于离散余弦变换DCT的鲁棒水印MATLAB脚本面向数字图像处理、多媒体安全与版权保护方向的学习者可用于理解频域水印的嵌入、提取与抗攻击原理。压缩包仅含1个m文件体积约1KB代码精简适合作为课程实验或毕业设计的参考实现。脚本通常涵盖图像分块DCT变换、水印信号叠加、逆变换重构以及水印检测等典型流程通过在低频或特定系数中嵌入信息来平衡视觉透明性与鲁棒性同时可能包含加密预处理和强度控制逻辑便于读者体会安全性、不可见性与抗干扰能力之间的权衡。借助该脚本可快速完成从DCT原理到实际代码的映射并验证缩放、压缩等攻击后的水印可检测性对深入理解系数选择、强度调节等关键设计环节很有帮助。目前已有147人学习下载对希望入门频域水印算法的初学者或需要快速搭建验证环境的开发者均有参考价值。1. 鲁棒水印脚本的切入点lubgiviyn.m 在频域里做什么打开 lubgiviyn.m第一眼并不容易判断是干什么的文件名像哈希截断没有 readme多半是从论文配套代码或者老仓库里拆出来的。把它丢进 MATLAB 之前要意识到这类脚本背后通常是同一套逻辑读入灰度图8×8 分块做 DCT在选中频系数上叠加一段 ±1 序列再 IDCT 回去得到含水印图鲁棒性测试时用缩放、旋转、JPEG 压缩把含水印图折腾一遍再根据系数差判断水印位是否存在。鲁棒水印的数字图像处理基础是频带可选择性——低频承载视觉主体高频在压缩和滤波时最先丢只有中频既不显眼又足够持久。这就是水印和条形码的区别看不见但经得起各种破坏。适合需要快速验证和实施数字图像版权保护图片溯源、版权标记的工程师和学生也适合用 MATLAB 做频域水印实验的人事后来补课。下面按模块拆解可运行的代码和参数边界。2. 分块 DCT 与频带选择从系数矩阵到嵌入位置在 DCT 域里图像要分成 8×8 小块原因不是顺手而是 JPEG 统一采用 8×8 分块水印与压缩攻击对齐能最大化存活率。DCT 把空间域 8×8 像素块变成同样尺寸的频率系数矩阵左上角 (1,1) 是直流分量代表块灰度平均值向右下角高频成分代表人眼不敏感的细节。2.1 用 blockproc 和 mat2cell 搭建分块变换MATLAB 里最顺手的是 blockproc它自动按块大小切图并对每块跑回调函数img im2double(imread(cameraman.tif)); % 转成双精度范围[0,1] dct_fun (block) dct2(block.data); % 每块做二维DCT dct_blocks blockproc(img, [8 8], dct_fun);代码执行后dct_blocks 仍然是整幅图尺寸实际相当于把所有块 DCT 后按原空间顺序拼贴。im2double 保证输入范围一致避免 uint8 类型引入过大量化误差。blockproc 的 [8 8] 是不可重叠分块适合直接复用 JPEG 的分块边界如果追求频域分辨率可以改成 [16 16]但块越大几何攻击后块间同步越容易错位。需要更细粒度操作每个块时我常用 mat2cell 把图像切成 cell 数组block_size 8; [r, c] size(img); blocks mat2cell(img, ones(1, r/block_size)*block_size, ... ones(1, c/block_size)*block_size); dct_cells cellfun((b) dct2(b), blocks, UniformOutput, false);这段代码将整张图划分为一个个 8×8 cellcellfun 逐个对块做 dct2返回的 dct_cells{i,j} 就是第 i 行第 j 块的 DCT 系数矩阵。这样后续改系数、查位置都比 blockproc 拼出来的一整张大矩阵直观。2.2 频带属性与可嵌入性对照不同位置的 DCT 系数在视觉权重、抗压缩能力上差异很大。以下表作为位置选择的参考频带块内大致位置视觉影响抗 JPEG 压缩抗噪声/滤波水印适用性直流/低频(1,1)、(2,1)、(1,2)强决定块亮度基色最持久较差可改但强度要小否则块效应明显中频(3,3)、(4,2)、(2,4)、(5,1)中等关联边缘方向较好中等综合首选高频(7,7)、(8,8)、(6,8)弱压缩后几乎丢失对噪声敏感仅适合超大 alpha 或高频密钥模式从表中可以看出水印嵌入位置不是越隐蔽越好而是要选攻击下残留信号仍然可读的位置。JPEG 是有损压缩量化表对高频系数量化步长大系数被丢弃概率高低频能量大但改动极易被人眼发现而且经过高通滤波时反而容易丢失。中频区间在工程上是最安全的折中。2.3 中频系数位置的选取与 Zigzag 展开选定中频位置不能拍脑袋。常见做法是把 8×8 系数按 Zigzag 顺序读成一维数组然后取第 6 到第 16 个索引对应的系数作为中频水印区。Zigzag 顺序模拟频率从低到高的过程zigzag_idx [1 2 9 17 10 3 4 11 18 25 ...]; % 实际补全64个索引 mid_positions zeros(10, 2); for n 1:10 idx zigzag_idx(n5); % 跳过直流和极低频 mid_positions(n, :) [ceil(idx/8), mod(idx-1,8)1]; end上面这段代码里 zigzag_idx 需要按标准路径列出 64 个整数第 6 到第 15 位取出来就是一组环绕在中频区的位置。这样取位的好处是位置自然按频率递增排序且不会无意中嵌入到直流系数。实际脚本里也可以直接用固定矩阵pos [3 3; 4 2; 2 4; 3 4; 4 3];来挑 5 个中频点效果类似。3. 水印嵌入与提取从系数修改到相关检测嵌入过程的核心是把二进制水印序列变成对 DCT 系数的小幅扰动。这个扰动必须小于视觉阈值同时大于后续攻击可能造成的噪声幅度。工程实现里水印序列通常经过扩频同一个 bit 被复制到多个块上提取时通过多数投票或相关计算恢复。3.1 生成水印序列并逐块嵌入以下代码使用固定随机种子生成 64 位水印然后按重复嵌入方式把每位水印加到每个块的 (3,3) 中频系数上rng(2024); wm_bits randi([0 1], 64, 1); % 64 位随机水印 wm_seq 2 * wm_bits - 1; % 0→-1, 1→1方便做正负扰动 alpha 0.15; % 嵌入强度需根据图内容调整 for bi 1:size(dct_cells, 1) for bj 1:size(dct_cells, 2) bit_idx mod((bi-1)*size(dct_cells, 2) bj-1, 64) 1; x dct_cells{bi, bj}(3, 3); % 取中频系数 dct_cells{bi, bj}(3, 3) x alpha * wm_seq(bit_idx); end end上面双重循环中第 9 行的取模运算让水印位按照行列顺序循环覆盖到所有块上。如果一个块对应一个 bit那么一个 32×32 分块图像能嵌入 1024 bit但为了鲁棒性通常把同一 bit 重复几十次用空间冗余换抗攻击性。alpha 直接控制扰动幅度取值在 0.05~0.2 附近超过 0.3 时视觉上会出现横竖网格噪声这是因为 8×8 块边界不连续。3.2 重构含水印图像系数改完后对每个块做 IDCT并把 cell 还原成图像矩阵watermarked_cells cellfun((b) idct2(b), dct_cells, UniformOutput, false); watermarked_img cell2mat(watermarked_cells); watermarked_img max(0, min(1, watermarked_img)); % 截断到 [0,1] watermarked_img im2uint8(watermarked_img);im2uint8 会把浮点图像量化回 8bit 整数这是模拟实际存储过程的必要步骤。截断是为了防止 IDCT 出现的轻微超界导致 im2uint8 饱和度异常。这里不直接保存原始浮点矩阵是因为真实使用场景里水印图像要以 JPG/PNG 形式分发JPEG 压缩本身就是最强的攻击。3.3 提取与检测盲提取做符号判断相关检测做存在性判定如果嵌入端保留原始图像或原始 DCT 系数提取端可以用差分方式恢复水印function wm_ext extract_diff(orig_blocks, att_blocks, alpha, wm_len) delta_sum zeros(wm_len, 1); cnt zeros(wm_len, 1); for bi 1:size(att_blocks, 1) for bj 1:size(att_blocks, 2) idx mod((bi-1)*size(att_blocks, 2) bj-1, wm_len) 1; diff att_blocks{bi,bj}(3,3) - orig_blocks{bi,bj}(3,3); delta_sum(idx) delta_sum(idx) diff / alpha; cnt(idx) cnt(idx) 1; end end wm_ext sign(delta_sum ./ cnt); % 平均后判符号1 或 -1 end这种做法的前提是攻击者没有修改原始图像、也没有对含水印图做全局平移剪切等破坏同步的操作。如果块对齐失效比如图像旋转 30°提取端按原网格取块时 (3,3) 已经不是嵌入时的那个系数diff 的平均值会趋向于 0误码率直接到 50%。所以现代鲁棒水印一定带同步模板或基于特征点的几何校正。对于不需要解码出完整水印、只需证明“此图含水印”的场景用相关性检测更稳定把提取到的系数序列和已知水印序列做归一化相关超过阈值就认为存在水印。function [nc, verdict] detect_by_corr(coeff_list, wm_seq) coeff_list coeff_list(:); wm_seq wm_seq(:); if length(coeff_list) ~ length(wm_seq) error(序列长度不一致); end nc abs(wm_seq * coeff_list) / (norm(wm_seq) * norm(coeff_list)); verdict nc 0.35; endcorr 计算的是投影余弦值当水印被 JPEG 压缩削弱但方向仍保持时NC 会从 0.9 掉到 0.4~0.6但阈值 0.35 仍能给出正判定。阈值不能设太低否则随机噪声也会超阈值。3.4 嵌入参数速查表不同应用场景需要调节的参数和影响面如下表可把它当成调试 lubgiviyn.m 时改参数的依据参数取值范围主要影响调试方向alpha0.05 ~ 0.2视觉透明度和鲁棒性先固定 0.1看 JPEG-q50 下的 NC再上下浮动块大小8×8 / 16×16同步粒度与频率分辨率默认 8×8若经常遇到旋转攻击可考虑 16×16水印长度16 ~ 256 bit信息量与冗余度版权只有 ID 时 32 bit 足够长水印要增强嵌入重复次数取决于块数/水印长抗攻击余量建议每个 bit 至少重复 20 次以上嵌入位置中频若干位置频带持久性高频失效时换 (2,3) 等低频位置代价是可见性这个表里的数值对应 DCT 量化表的步长和人类视觉系统的对比度灵敏度不是随便拍的。4. 攻击模拟与鲁棒性度量把含水印图折腾一遍再算数鲁棒性不是概念是测量出来的。攻击模拟必须和嵌入提取代码放在同一个工程目录先攻击后提取算 NC 和 BER再决定是否调整参数。4.1 常见攻击的 MATLAB 模拟代码下面这个函数是一个最小攻击集覆盖几何和信号处理两大类。实际项目中还应加入高斯低通滤波、泊松噪声、密度噪声等。function attacked attack_watermarked(img, attack_name) switch attack_name case scale small imresize(img, 0.5); % 缩小一半 attacked imresize(small, size(img)); % 恢复原尺寸后再提取 case rotate attacked imrotate(img, 30, bicubic, crop); % 旋转30度裁边补零 case crop attacked imcrop(img, [50 50 199 199]); % 剪去边缘区域 attacked imresize(attacked, size(img)); % 拉伸回原尺寸 case gauss attacked imnoise(img, gaussian, 0, 0.01); % 均值0方差0.01 case jpeg imwrite(img, tmp_q50.jpg, Quality, 50); % 质量因子50 attacked imread(tmp_q50.jpg); end end缩放攻击恢复尺寸是提取端惯用处理因为水印检测一般假设图像有尺寸先验如果检测端不知道原尺寸那需要训练尺度盲检测或者嵌入尺度和方向模板。rotate 用 bicubic 插值模拟更严重的网络传输压缩效果crop 参数选择裁剪中心区域再拉伸模拟用户截屏或平台强制缩放。4.2 NC 和 BER两个标准量化指标解码后水印要和原始水印逐位比较不能只凭肉眼。NC 衡量整段序列的相关程度BER 衡量判错的比特比例。function nc normalized_correlation(wm_a, wm_b) a wm_a(:); b wm_b(:); nc sum(a .* b) / sqrt(sum(a.^2) * sum(b.^2)); end function ber bit_error_rate(wm_a, wm_b) ber sum(wm_a(:) ~ wm_b(:)) / numel(wm_a); endNC 在 [0,1] 之间NC1 表示完全一样0 表示不相关。BER 常用百分数表示20% 以内的误码率在水印纠错码存在时可以恢复原文。注意如果水印序列是随机数单纯看 NC 不够还要跑几百次随机序列算基准 NC 分布才能确定 0.35 阈值是否可靠。4.3 批量跑攻击矩阵并输出报告工程上不会只跑一张图会把 LENA 这类标准图和多张实拍图一起放进循环得到统计意义上的鲁棒性曲线attacks {scale, rotate, crop, gauss, jpeg}; fprintf(%-10s %-8s %-8s\n, attack, NC, BER); for k 1:length(attacks) attacked attack_watermarked(watermarked_img, attacks{k}); att_cells cellfun((b) dct2(b), mat2cell(attacked, ... ones(1, size(attacked,1)/8)*8, ones(1, size(attacked,2)/8)*8), ... UniformOutput, false); wm_extracted extract_diff(dct_cells, att_cells, alpha, 64); nc normalized_correlation(wm_bits, wm_extracted); ber bit_error_rate(wm_bits, wm_extracted); fprintf(%-10s %-8.3f %-8.3f\n, attacks{k}, nc, ber); end这里攻击后重新分块做 DCT与原始块的 DCT 对齐后差分。打印结果时可以直接复制进 Excel 做折线图。如果 rotate 攻击的 NC 跌到 0.3 以下说明嵌入块和提取块的坐标没有对齐这不是 alpha 问题而是同步问题如果 jpeg q50 的 NC 略低但仍在 0.6 以上那提升提取性能的空间主要在嵌入强度上。5. 自适应嵌入强度按块纹理动态调整 alpha 的改造技巧直接对全图用同一个 alpha平坦区域容易露出块状痕迹纹理区域又浪费了可用强度。一个实用技巧是计算每个块的局部纹理度量然后用纹理归一化系数缩放 alpha。DCT 块内中高频能量可以作为纹理粗尺度估计alpha_base 0.12; textures zeros(size(dct_cells)); for bi 1:size(dct_cells, 1) for bj 1:size(dct_cells, 2) block dct_cells{bi,bj}; high_energy sum(block(4:8, 4:8).^2, all); % 中高频区能量 textures(bi, bj) high_energy; end end textures textures / max(textures(:)); % 归一化到 [0,1] for bi 1:size(dct_cells, 1) for bj 1:size(dct_cells, 2) bit_idx mod((bi-1)*size(dct_cells,2) bj-1, 64) 1; alpha_local alpha_base * (0.6 0.8 * textures(bi, bj)); dct_cells{bi,bj}(3,3) dct_cells{bi,bj}(3,3) alpha_local * wm_seq(bit_idx); end endhigh_energy 排除左上角低频区把中高频系数能量作为纹理深浅的度量纹理丰富但系数能量高的区域修改系数不易被察觉因此 alpha 放大多给一些。0.60.8*textures 保证最平坦区域的 alpha 也有 alpha_base 的六成不会因为纹理极低导致水印嵌入太弱。改造完之后要回去跑第 4 章的攻击矩阵重点对比固定 alpha0.12 和自适应 alpha 在 JPEG q50 下的 NC如果自适应版本的平均 NC 提升了 0.05 以上且原图 PSNR 没有下降超过 1dB说明强度分配是有效的。调试时在嵌入代码前加上 tic/toc 统计耗时纹理计算在超大图上可能成为瓶颈必要时改用 blockproc 的并行输出重写避免用两层循环逐块扫描。本文还有配套的精品资源点击获取