企业官网建设流程全解析

简介：面向无线通信与物理层安全研究，这份代码针对智能反射表面辅助的多输入多输出系统，利用坐标下降算法最大化保密率，适用于课程设计、期末大作业与毕业设计。资源共12个文件，其中6个.m脚本分别承担主程序调度、智能反射表面信道建模、坐标下降核心迭代逻辑、固定相位偏移下容量计算、一般特征值分解容量求解以及穷举搜索对比，另有3个.zbak备份文件、1个zip案例数据包、1份README说明文档和1个gitattributes配置，整体仅9KB，结构轻量清晰。代码支持MATLAB 2014、2019a与2021a，采用参数化编程，所有参数均可按需调整，注释详尽，附赠案例数据可直接运行以观察算法收敛效果。已有37人学习，适合计算机、电子信息、数学等专业学生通过修改参数对比不同配置下的保密率表现，结合README中的说明理解坐标下降在IRS-MIMO物理层安全优化中的实际价值。

1. 为什么 IRS 辅助 MIMO 的保密率优化，默认解法是坐标下降

IRS 辅助 MIMO 系统的保密率最大化，表面看是个标准的多变量非凸优化问题，真正在 MATLAB 里跑起来才会意识到，发射协方差、IRS 相移、接收合并三个变量块互相锁死，直接联合优化连个像样的初始点都难找。坐标下降算法之所以成为这个方向绝大多数论文的默认基线，是因为它把问题拆成“固定相移更新协方差、固定协方差更新相移”两个子问题轮流解，每个子问题要么有闭式解，要么只需要在流形上做一步梯度。这篇笔记从系统模型和保密率表达式讲起，一路落到可复现的 MATLAB 实现：注水功率分配、相移数值梯度、收敛判定、参数设置，以及四个高频踩坑。适合正在复现 IRS 物理层安全论文、或者准备把保密率优化写进研究方案和开题报告的从业者。

2. 系统模型与保密率目标：先把三个变量块的耦合关系拆清楚

2.1 三跳信道与 IRS 相移矩阵的建模约定

先把信道模型定死，后面所有代码都基于这套约定。发射端是一个配置 (N_t) 根天线的基站，IRS 上有 (N) 个反射单元，合法用户配 (N_u) 根天线，窃听者配 (N_e) 根天线。窄带平坦衰落假设下，需要生成五个信道矩阵：基站到 IRS 的 (G)（尺寸 (N \times N_t)），IRS 到合法用户的 (H_{ru})（(N_u \times N)），IRS 到窃听者的 (H_{re})（(N_e \times N)），基站到合法用户的直达链路 (H_{du})（(N_u \times N_t)），以及基站到窃听者的直达链路 (H_{de})（(N_e \times N_t)）。这个模型是 IRS 辅助物理层安全文献里最标准的框架，区别只在信道是否加莱斯因子、是否把 IRS 单元换成离散相位。

IRS 的作用全部浓缩在一个对角矩阵里：(\Theta = \mathrm{diag}(\theta_1, \dots, \theta_N))，每个 (\theta_n) 是单位模复数，代表第 n 个反射单元的相移。无源 IRS 只调相位，不放大信号，所以恒模约束 (|\theta_n| = 1) 是硬约束。把三条链路叠加起来，合法用户和窃听者的等效信道分别是

[ H_{eq,u} = H_{du} + H_{ru} \Theta G, \qquad H_{eq,e} = H_{de} + H_{re} \Theta G ]

这个式子值得多看两眼：IRS 路径的信号先经 (G) 到达反射面，被 (\Theta) 逐单元移相，再经 (H_{ru}) 或 (H_{re}) 进入接收端，与直达信号相干叠加。所谓“智能反射”，本质就是调节 (\theta) 的相位组合，让有用信号在合法用户方向同相叠加、在窃听者方向反相抵消。

2.2 从互信息到保密率：目标函数的非凸性来源

在发射协方差矩阵 (Q)（(N_t \times N_t)，半正定，迹不超过总功率 (P)）给定后，合法用户的可达速率是

[ R_u(Q,\theta) = \log_2 \det\left(I + \frac{H_{eq,u} Q H_{eq,u}^H}{\sigma_u^2}\right) ]

窃听者同理，把下标换成 e。这其实就是 MIMO 信道容量公式在具体等效信道上的展开，对应接收端做 MMSE 合并、发送端高斯输入的情形。保密率定义为两个速率之差取非负：

[ R_s = \max{0,\ R_u - R_e} ]

单位是 bit/s/Hz。这个 (\max{0,\cdot}) 不是个平滑函数，它在 (R_u = R_e) 处存在不可导点，后文讲梯度更新时会回来踩这个坑。

为什么说这个优化问题难？三个原因叠在一起。第一，(R_u - R_e) 是凹函数减凹函数，整体非凹，直接最大化没有凸优化工具可用。第二，(\theta) 的恒模约束是复平面单位圆上的流形，不是凸集，投影操作也不直观。第三，(Q) 除了半正定约束还要满足迹约束，和相移变量耦合在一起时，联合可行域的形状非常差。坐标下降能成为默认基线，不是因为它能保证全局最优，而是因为拆出来的每个子问题都足够“乖”：固定 (\theta) 时，(Q) 子问题可以落回注水框架；固定 (Q) 时，(\theta) 子问题只需要梯度信息。

2.3 变量耦合分析：为什么要拆成子问题而不是联合梯度

一个更直接的问题是：既然 (R_s) 对 (Q) 和 (\theta) 都有梯度，为什么不用一个联合梯度下降把所有变量一起更新？答案是两个变量块的“节奏”完全不同。(Q) 是连续半正定矩阵，有效更新方式是 SVD 域上的功率分配；(\theta) 生活在单位模流形上，有效更新方式是相位旋转。把它们硬拼成一个长向量，步长稍大就同时违反功率约束和恒模约束，步长稍小就慢到没法收敛。坐标下降把两个约束解耦，每个子问题只维护自己的结构，这是它能稳定工作的根本原因。

还有一个工程上的现实：IRS 单元数通常远大于天线数。(N) 从 16 加到 256，(\theta) 的维度比 (Q) 的维度高一个数量级，联合梯度里 (Q) 的更新很容易被相移梯度的噪声淹没。分开更新之后，每次迭代的收益来源很清楚，调试时能直接判断是相移在贡献增益还是功率分配在贡献增益。这个可解释性，在写论文、改参数、跟审稿人解释收敛曲线时都非常值钱。

3. 坐标下降框架与 MATLAB 函数实现：注水、流形梯度与收敛判定

3.1 整体循环结构与三个子问题的分工

整个坐标下降循环只有两步：固定 (\theta) 更新 (Q)，固定 (Q) 更新 (\theta)。选代顺序也有讲究，我一般先更新 (Q) 再更新 (\theta)。原因是 (Q) 的注水解会立刻把相移更新的收益反映到功率分配上，如果先更新 (\theta)，这一轮算出来的保密率里还带着旧 (Q)，迭代初期的收敛速度会慢一截。两步各有各的工具：(Q) 子问题用 SVD 加注水，(\theta) 子问题用流形梯度上升。下面三个小节给出可以直接存成.m文件的实现。

3.2 固定相移更新发射协方差：SVD 域注水

固定 (\theta) 后，先把 IRS 的效果合进用户等效信道 (H_{eq,u})，然后对 (H_{eq,u}) 做奇异值分解 (H_{eq,u} = U S V^H)，令 (Q = V P V^H)，其中 (P) 是对角功率分配矩阵。为什么限定这个结构？因为固定 (\theta) 后，用户速率中与 (Q) 有关的项在 (Q) 与用户等效信道右奇异向量对齐时，退化成一组平行子信道的功率分配问题。每个子信道的增益是奇异值平方 (s_i^2)，功率分配用经典注水：

function p = water_filling_power(sigma_sq, Ptot, sigma2) % 输入 % sigma_sq : 用户等效信道各子信道增益（奇异值平方，列向量） % Ptot : 总发射功率（W） % sigma2 : 合法用户噪声功率（W） % 输出 % p : 各子信道分配的功率 sigma_sq = max(sigma_sq, 1e-12); % 防止奇异值为 0 时除零告警 mu_max = max(sigma2 ./ sigma_sq) + Ptot; % 水位上界，保证功率超限 mu_min = 0; for k = 1:50 mu = (mu_min + mu_max) / 2; p = max(0, mu - sigma2 ./ sigma_sq); if sum(p) > Ptot mu_max = mu; else mu_min = mu; end end p = p * (Ptot / sum(p)); % 归一化到总功率 end

这个函数把注水问题用二分法求水位 (\mu)，返回的各子信道功率满足 (\sum p_i = P_{\mathrm{tot}})。二分 50 次足够收敛到浮点精度，再多没有意义。归一化那一行是为了消掉二分残余误差，实际功率偏差通常在 1e-10 量级，不归一化也不影响结果，但归一化能避免后续 (Q) 的迹出现微小超差。注意这里对窃听者链路的处理是近似的，只优化了用户信道方向，在窃听者信道明显弱于用户、或者 IRS 还没收敛到窃听者方向时足够用；等遇到强窃听信道时这个近似会让算法翻车，第五章专门讲。

3.3 固定协方差更新相移：流形梯度上升

固定 (Q) 后，保密率是 (\theta) 的连续函数，但直接做“(\theta = \theta + \mu g) 再投影回单位圆”非常容易把相位弄乱。更干净的做法是直接在流形上移动：(\theta_n \leftarrow \theta_n \exp(j \mu g_n))。这样更新后的 (\theta) 自动满足恒模约束，不需要额外的投影步骤。梯度怎么算？在小规模 IRS 上最省事的是数值梯度，用流形上的中心差分：

function g = numerical_grad_theta(theta, ch, Q, sigma2_u, sigma2_e) % 输入 % theta : 当前 IRS 相移向量（N x 1，单位模复数） % ch : 打包好的信道结构体，包含 G/H_ru/H_re/H_du/H_de % Q : 当前发射协方差矩阵 % 输出 % g : 保密率对每个相移的梯度（实数） d = length(theta); delta = 1e-6; % 有限差分步长，单位弧度 g = zeros(d, 1); for n = 1:d theta_plus = theta; theta_plus(n) = theta(n) * exp(1i * delta); theta_minus = theta; theta_minus(n) = theta(n) * exp(-1i * delta); R_plus = secrecy_rate(theta_plus, ch, Q, sigma2_u, sigma2_e); R_minus = secrecy_rate(theta_minus, ch, Q, sigma2_u, sigma2_e); g(n) = (R_plus - R_minus) / (2 * delta); end end function R = secrecy_rate(theta, ch, Q, sigma2_u, sigma2_e) Theta = diag(theta); H_eq_u = ch.H_du + ch.H_ru * Theta * ch.G; H_eq_e = ch.H_de + ch.H_re * Theta * ch.G; R_u = real(log2(det(eye(size(H_eq_u,1)) + (H_eq_u * Q * H_eq_u') / sigma2_u))); R_e = real(log2(det(eye(size(H_eq_e,1)) + (H_eq_e * Q * H_eq_e') / sigma2_e))); R = max(0, R_u - R_e); end

数值梯度的原理是对每个反射单元的相位分别加减小量 (\delta)，用两次保密率的差值近似导数。这里 (\delta) 取 1e-6 弧度，足够接近真导数。(\delta) 取 1e-3 以上时差分误差会明显变大，梯度带噪声，迭代曲线开始震荡——这个坑后面还会展开。由于保密率计算只涉及两个 (2 \times 2) 矩阵的行列式，(N=32) 时一轮 65 次速率计算在 MATLAB 里是毫秒级；等 (N) 上到 256 再考虑换解析梯度，第五章给出具体策略。

提示：数值梯度里的有限差分步长 (\delta) 和后面迭代用的优化步长 (\mu) 是两个完全不同的量，前者是 1e-6 量级的微扰，后者是 0.01 到 0.05 量级的移动，别混用。

3.4 迭代终止与收敛判定的代码实现

有了注水函数和梯度函数，坐标下降主循环就非常短了：

function [R_hist, Q, theta] = cd_solver(ch, Ptot, sigma2_u, sigma2_e, N_irs, Ntx, max_iter, tol) % 坐标下降主循环 % 输入参数在调用方设置，max_iter 建议 100，tol 建议 1e-4 theta = exp(1i * 2*pi*rand(N_irs, 1)); % 随机相移初值 Q = (Ptot / Ntx) * eye(Ntx); % 均匀功率初值 R_hist = zeros(max_iter, 1); for iter = 1:max_iter % 固定 theta，更新 Q：用户等效信道做 SVD，注水 H_eq_u = ch.H_du + ch.H_ru * diag(theta) * ch.G; [~, S, V] = svd(H_eq_u, 'econ'); p = water_filling_power(diag(S).^2, Ptot, sigma2_u); Q = V * diag(p) * V'; % 固定 Q，更新 theta：流形梯度上升 g = numerical_grad_theta(theta, ch, Q, sigma2_u, sigma2_e); step = 0.05 / (max(abs(g)) + eps); % 把最大相位变化压到 0.05 弧度内 theta = theta .* exp(1i * step * g); R_hist(iter) = secrecy_rate(theta, ch, Q, sigma2_u, sigma2_e); if iter > 3 && abs(R_hist(iter) - R_hist(iter-1)) < tol break; end end R_hist = R_hist(1:iter); end

step 归一化是关键：把最大相位变化压到 0.05 弧度以内，收敛曲线平滑，不会因为某一步梯度特别大而越过峰值。0.05 这个数不是拍脑袋，相移对保密率的敏感度接近正弦函数形状，单步相位超过 0.1 弧度容易跳过极值点；保守一点换来的是更稳定的收敛行为。如果觉得收敛太慢，可以调到 0.1，但要多配几次随机初始化。收敛判定用相邻两轮保密率差值小于 tol，tol 取 1e-4 就够，再小对浮点计算没有意义。把上面几个函数分别存成water_filling_power.m、secrecy_rate.m、numerical_grad_theta.m和cd_solver.m，就是一份能跑的最小实现。

4. 用 MATLAB 主脚本跑通最小可复现示例：参数、代码与结果解读

4.1 仿真参数怎么设

参数设置直接影响收敛速度和保密率的绝对数值。我常用的这套参数贴近论文里的典型配置：基站 4 天线，用户和窃听者各 2 天线，IRS 从 16 到 64 单元扫。发射功率 30 dBm，噪声功率 -90 dBm，信噪比足够高但不会高到注水退化成均匀功率分配。

这里用固定衰减系数占位，实际工程场景应该换成路径损耗公式，比如 (L = L_0 (d/d_0)^{-\alpha})。窃听者的 IRS 链路衰减设得比用户低，是为了模拟窃听者位置更远、反射路径损耗更大的情形，保密率算出来才有区分度。直达链路衰减设 0.1，是为了让 IRS 反射路径成为主要信号来源，否则 IRS 的作用不明显。

4.2 主脚本完整实现

%% 1. 参数区 clear; clc; rng(42); % 固定随机种子，保证可复现 Ntx = 4; Nrx = 2; Nev = 2; N_irs = 32; Ptot = 1e-3 * 10^(30/10); % 30 dBm 转成 W sigma2_u = 1e-9; sigma2_e = 1e-9; % -90 dBm beta_bs_irs = 0.9; beta_irs_user = 0.6; beta_irs_eve = 0.3; beta_bu = 0.1; beta_be = 0.1; %% 2. 信道生成（瑞利平坦衰落） G = sqrt(beta_bs_irs) * (randn(N_irs,Ntx) + 1i*randn(N_irs,Ntx)) / sqrt(2); H_ru = sqrt(beta_irs_user) * (randn(Nrx,N_irs) + 1i*randn(Nrx,N_irs)) / sqrt(2); H_re = sqrt(beta_irs_eve) * (randn(Nev,N_irs) + 1i*randn(Nev,N_irs)) / sqrt(2); H_du = sqrt(beta_bu) * (randn(Nrx,Ntx) + 1i*randn(Nrx,Ntx)) / sqrt(2); H_de = sqrt(beta_be) * (randn(Nev,Ntx) + 1i*randn(Nev,Ntx)) / sqrt(2); ch.G = G; ch.H_ru = H_ru; ch.H_re = H_re; ch.H_du = H_du; ch.H_de = H_de; %% 3. 调用坐标下降求解 max_iter = 100; tol = 1e-4; [R_hist, Q_opt, theta_opt] = cd_solver(ch, Ptot, sigma2_u, sigma2_e, N_irs, Ntx, max_iter, tol); %% 4. 对照组：随机相移 theta_rand = exp(1i * 2*pi*rand(N_irs, 1)); R_rand = secrecy_rate(theta_rand, ch, Q_opt, sigma2_u, sigma2_e); fprintf('随机相移保密率: %.4f bit/s/Hz\n', R_rand); fprintf('坐标下降保密率: %.4f bit/s/Hz\n', R_hist(end)); fprintf('实际迭代次数: %d\n', length(R_hist)); %% 5. 画收敛曲线 figure; plot(1:length(R_hist), R_hist, 'o-', 'LineWidth', 1.5); xlabel('迭代次数'); ylabel('保密率 (bit/s/Hz)'); grid on; title('坐标下降收敛曲线');

信道生成里所有矩阵都除以 (\sqrt{2})，是为了让每个信道系数的能量期望为 1，配合衰减系数后链路增益量级合理。结构体 ch 打包五个信道矩阵，避免后续函数传参写成长长的参数列表，这个习惯在信道矩阵多的时候非常省事。rng(42) 固定种子是关键，否则每次跑出来的保密率数字都不同，没法判断改参数到底有没有效。

4.3 运行结果长什么样：收敛曲线与典型数字

跑完会看到一条先快速上升、然后进入平台的收敛曲线。前 5 轮提升最明显，之后每轮只有零点零几的增益，这是坐标下降的典型行为——前几步把最明显的方向性增益吃干净，后面在细调相位。一次随机信道实现下，结果大致长这样：

换随机种子数字会浮动，但趋势是固定的：IRS 单元数从 16 加到 64，保密率能拉开 2 到 3 倍的距离，继续加到 128、256，边际增益明显放缓。这个表格也可以当自检标准——如果 32 单元下坐标下降结果跟随机相移差不多，多半是数值梯度或注水环节出了 bug，而不是算法本身有问题。还可以顺手把 (R_u) 和 (R_e) 各自画出来，观察差值曲线。用户信道和窃听者信道接近时，(R_s) 会出现零值平台，这不是 bug，是 (\max{0,\cdot}) 截断在起作用。写论文画图时我习惯把两条速率曲线和保密率曲线叠在一起，否则审稿人会问中间那段为什么是零。

5. 保密率坐标下降的四个经典踩坑：从现象到排查

5.1 迭代几轮后保密率不升反降

现象：收敛曲线前几轮正常上升，第四五轮突然掉头向下，之后虽然反弹但整体收敛质量很差。

原因：固定 (\theta) 更新 (Q) 时，注水只看了用户等效信道。如果窃听者在空间上恰好落在用户主特征方向附近，注水把功率集中到用户信道主方向的同时，也把窃听者的接收功率抬上去了，(R_e) 的增速超过 (R_u)，保密率反而下降。

解决：把 (R_u) 和 (R_e) 的逐轮记录打出来，确认是不是窃听者速率在涨。如果是，说明当前信道实现下“对用户信道注水”这个近似不够好。最简单的处理是这一步固定 (Q) 不更新，只让 IRS 相移去调整空间方向；更彻底的方案是用 WMMSE 更新 (Q)，把窃听者项也纳入目标，代价是每轮多几次矩阵迭代。我在代码里保留注水版本作为基线，但遇到强窃听场景会直接切 WMMSE。

5.2 迭代曲线震荡：有限差分步长和优化步长别混用

现象：保密率曲线锯齿状震荡，相位向量在连续几轮里跳来跳去，最后收敛到明显低于正常水平的解。

原因：两个参数同时出错。一是数值梯度的有限差分步长 (\delta) 设成了 1e-3 这类偏大的值，中心差分结果带着明显误差；二是迭代步长 step 没有按梯度模归一化，某个相位的相位变化超过 0.1 弧度，直接越过峰值。

解决：(\delta) 固定为 1e-6，这是差分近似精度和浮点误差的平衡点。step 用0.05 / max(abs(g))做归一化，保证单次最大相位变化不超过 0.05 弧度。调试时先把梯度的模打出来看量级，如果某个单元的梯度比其他单元大两个数量级以上，多半不是物理问题，而是差分公式写错了。

5.3 迭代全程保密率等于零：初始点选错方向

现象：跑完 100 轮，R_hist 全是 0，坐标下降像没启动一样。

原因：初始 (\theta) 是随机相位，如果恰好把 IRS 的反射能量聚焦到窃听者方向，(R_e) 会远大于 (R_u)，此时 (R_s = \max{0, R_u - R_e} = 0)。问题在于这个不可导点在梯度算法看来是个“死区”：(R_u < R_e) 时速率差为负，梯度近似为零，坐标下降感知不到任何改进方向。

解决：换初始化策略。我一般同时跑三个初始点：全 1 相移、随机相移、以及把 IRS 相位对准合法用户信道方向的匹配相移（取 (H_{ru}) 和 (G) 的主特征向量做相位对齐），最后取保密率最高的结果。多初始点在坐标下降里不是可选项，是必选项，因为算法本身没有全局收敛保证。每次跑完把三个初始点的收敛曲线画在一起，还能顺便判断问题是不是凸的、有多少个局部极值。

5.4 N 上到 256，数值梯度慢到怀疑人生：换解析梯度

现象：32 单元时一轮迭代不到 1 秒，改成 256 单元后每轮要算 513 次保密率，每次里面有两个行列式，一轮十几秒，整体慢到没法调参。

原因：数值梯度的复杂度是 (O(N)) 次速率计算，而每次速率计算里有两个矩阵行列式。IRS 单元数线性增长，耗时也线性增长，但 (N=256) 时线性系数已经很大了。

解决：换成解析梯度。核心思路是利用 (\log\det(\cdot)) 的导数公式：(\partial \log_2\det(X) = \mathrm{tr}(X^{-1} \partial X))。把 (H_{eq,u} = H_{du} + \sum_n \theta_n a_n b_n^H) 代入，其中 (a_n) 是 (H_{ru}) 的第 n 列、(b_n) 是 (G) 的第 n 行转置，(\partial \Phi_u / \partial \theta_n) 会化简成秩一结构，最终变成一个向量内积，整个梯度的代价从 (O(N)) 次行列式降到 (O(N)) 次矩阵乘向量。这个推导很容易在常数因子和共轭转置上出错，我的做法是先保留数值梯度做对照，解析梯度写完跑一次数据对拍，两者相对误差小于 1e-6 才正式替换。这个对拍脚本留下来，后面改信道模型、加莱斯因子、换路径损耗配置时还能复用。

6. 从基线走向实战：离散相移、WMMSE 与多用户扩展

坐标下降跑通只是第一步，论文和工程实现里还有三件绕不开的事。第一件是 IRS 离散相移。实际反射单元只能取有限个相位状态，常见是 1 bit 到 3 bit。把连续解量化到 2 bit 的写法就一行：

phase_q = pi/2 * round(angle(theta_opt) / (pi/2)); theta_q = exp(1i * mod(phase_q, 2*pi)); R_q = secrecy_rate(theta_q, ch, Q_opt, sigma2_u, sigma2_e);

量化后保密率通常损失 10% 到 20%，具体看信道实现。评估量化损失时，拿连续解做上界、量化解做可达点，两个数字一起写进论文，这是审稿人最想看到的对比。第二件是把注水更新 (Q) 升级成 WMMSE。注水只在窃听者链路弱时可靠，强窃听场景下 WMMSE 把窃听者速率纳进目标，通过 (Q) 和接收合并矩阵交替更新逼近局部最优，代价是每轮多几次矩阵求逆，但能解决第五章第一个坑。多用户场景的思路也是这条线：把系统速率换成加权和速率，IRS 相移和所有用户的协方差交替更新，复杂度可控。第三件是保留一个验证习惯——每次改信道模型、加参数、换初始点，先跑数值梯度对拍，再跑随机相移对照，确认每轮收益的来源。我到现在写新模块还是先过这两关，不然参数一多，黑匣子就出来了。希望帮到你。

本文还有配套的精品资源，点击获取