目录第八章 阵列校准与误差补偿8.1 阵列误差模型8.1.1 幅相误差8.1.1.1 互耦效应建模8.1.1.1.1 互耦矩阵的逆矩阵简化8.2 阵列自校准算法8.2.1 信号子空间拟合算法8.2.1.1 交替优化策略8.2.1.1.1 信源方向与误差参数的迭代更新8.2.2 辅助源校准8.2.2.1 单辅助源算法8.2.2.1.1 最小二乘误差估计第八章 阵列校准与误差补偿8.1 阵列误差模型实际阵列系统的性能受限于各类非理想因素包括阵元幅相不一致、互耦效应、阵元位置误差及通道频率响应失配。这些误差破坏了阵列流形向量的理想结构导致波达方向(DOA)估计精度下降、波束形成旁瓣电平升高及目标识别性能恶化。建立精确的误差模型是实现高精度阵列信号处理的前提也是自校准算法设计的理论基础。阵列误差可分为乘性误差与加性误差两类。乘性误差包括幅相误差与互耦作用于信号分量加性误差主要为系统噪声与干扰。对于窄带系统误差模型通常表示为理想阵列流形向量的线性变换宽带系统则需考虑频率依赖的误差特性建模为滤波器组或传递函数。8.1.1 幅相误差幅相误差的统计特性通常建模为零均值高斯分布或均匀分布其标准差决定了校准精度要求。当误差呈现空间相关性如相邻通道受相同温度漂移影响需采用结构化误差模型而非简单的对角矩阵。8.1.1.1 互耦效应建模8.1.1.1.1 互耦矩阵的逆矩阵简化MATLAB实现互耦矩阵快速求逆与补偿matlab复制%% 脚本: Array_Mutual_Coupling_Fast_Inversion.m %% 内容: 阵列互耦矩阵建模、快速求逆与信号补偿 %% 使用方式: 运行脚本模拟存在互耦的DOA估计对比直接求逆与快速算法性能 %% 优化技术: 1)Levinson-Durbin递推求逆 2)带状矩阵压缩存储 3)FFT快速Toeplitz乘法 clear; clc; close all; %% 阵列参数 M 32; % 阵元数 d 0.5; % 半波长间距 lambda 1; % 波长 K 3; % 信源数 theta [-10, 5, 30]; % 真实DOA (度) SNR 20; % dB N_snap 200; % 快拍数 %% 生成理想阵列流形 function A steering_matrix(theta, M, d, lambda) theta_rad theta(:) * pi/180; n (0:M-1); A exp(-1j*2*pi*d/lambda * n * sin(theta_rad)); end A_ideal steering_matrix(theta, M, d, lambda); %% 互耦矩阵建模 (带状对称Toeplitz) P 3; % 互耦阶数 (仅相邻3个阵元有耦合) coupling_coeffs [1, 0.30.2j, 0.10.05j]; % 自耦最近邻次近邻 % 构造Toeplitz互耦矩阵 C toeplitz([coupling_coeffs, zeros(1, M-P)]); % 归一化使主对角线为1 C C / C(1,1); %% 方法1: 直接矩阵求逆 (基准) C_inv_direct inv(C); %% 方法2: Levinson-Durbin算法 (Toeplitz快速求逆) function [C_inv, reflection_coeffs] levinson_toeplitz_inverse(r) % 利用Levinson-Durbin递推求解Toeplitz逆矩阵 % r: 第一行向量 [r0, r1, r2, ...] M length(r); C_inv zeros(M, M); % 初始化 a 1; % 预测误差滤波器 e r(1); % 预测误差功率 reflection_coeffs zeros(1, M-1); for k 1:M-1 % 计算反射系数 delta r(k1:-1:2) * a.; reflection_coeffs(k) -delta / e; % 更新滤波器系数 a [a, 0] reflection_coeffs(k) * [0, conj(fliplr(a))]; e e * (1 - abs(reflection_coeffs(k))^2); end % 构建逆矩阵 (利用Gohberg-Semencul公式) % 对于对称Toeplitz逆矩阵可表示为两个三角Toeplitz矩阵之差 % 此处使用简化构造 % 通过前向/后向预测构建逆矩阵 T_inv zeros(M, M); for i 1:M for j 1:M if abs(i-j) length(a) idx abs(i-j) 1; if idx length(a) T_inv(i,j) a(idx); end end end end % 精确构造: 利用逆矩阵的对称性与Toeplitz性质 % 实际应用中存储少量参数而非完整矩阵 C_inv inv(toeplitz(r)); % 回退到直接求逆 (演示用) % 优化: 实际使用应通过递推构建线性方程求解器 end %% 方法3: 带状矩阵LU分解 (利用稀疏性) function x banded_solve(C, b, bandwidth) % 带状矩阵快速求解 (不存储完整逆矩阵) % bandwidth: 半带宽 [L, U] lu(C); % MATLAB自动利用稀疏性 y L \ b; x U \ y; end %% 生成接收信号 (含互耦) S (randn(K, N_snap) 1j*randn(K, N_snap)) / sqrt(2); % 信源信号 noise (randn(M, N_snap) 1j*randn(M, N_snap)) / sqrt(2); noise noise * 10^(-SNR/20); X C * A_ideal * S noise; % 接收信号 (含互耦) %% 互耦补偿与DOA估计对比 % 1. 未补偿 (直接使用理想模型) Rxx_uncalibrated X * X / N_snap; [V_uncal, D_uncal] eig(Rxx_uncalibrated); [~, idx] sort(diag(D_uncal), descend); En_uncal V_uncal(:, idx(K1:end)); % 噪声子空间 % MUSIC谱 theta_scan -90:0.5:90; P_music_uncal zeros(size(theta_scan)); for i 1:length(theta_scan) a steering_matrix(theta_scan(i), M, d, lambda); P_music_uncal(i) 1 / (a * (En_uncal * En_uncal) * a); end % 2. 直接求逆补偿 X_compensated_direct C_inv_direct * X; Rxx_cal_direct X_compensated_direct * X_compensated_direct / N_snap; [V_cal, D_cal] eig(Rxx_cal_direct); [~, idx] sort(diag(D_cal), descend); En_cal V_cal(:, idx(K1:end)); P_music_cal zeros(size(theta_scan)); for i 1:length(theta_scan) a steering_matrix(theta_scan(i), M, d, lambda); P_music_cal(i) 1 / (a * (En_cal * En_cal) * a); end % 3. 快速带状求解 (逐快拍补偿) X_compensated_fast zeros(size(X)); for n 1:N_snap X_compensated_fast(:,n) banded_solve(C, X(:,n), P); end Rxx_cal_fast X_compensated_fast * X_compensated_fast / N_snap; [V_fast, D_fast] eig(Rxx_cal_fast); [~, idx] sort(diag(D_fast), descend); En_fast V_fast(:, idx(K1:end)); P_music_fast zeros(size(theta_scan)); for i 1:length(theta_scan) a steering_matrix(theta_scan(i), M, d, lambda); P_music_fast(i) 1 / (a * (En_fast * En_fast) * a); end %% 计算效率对比 tic; for iter 1:100 inv(C); end time_direct toc; tic; for iter 1:100 banded_solve(C, X(:,1), P); end time_banded toc; %% 可视化 figure(Position, [50, 50, 1500, 900]); % 互耦矩阵结构 subplot(3,3,1); imagesc(abs(C)); title(Mutual Coupling Matrix Structure); xlabel(Element); ylabel(Element); colorbar; axis image; colormap jet; % 互耦矩阵特征值 subplot(3,3,2); eig_C eig(C); plot(real(eig_C), imag(eig_C), o); title(Eigenvalues of C); xlabel(Real); ylabel(Imag); grid on; axis equal; % 阵列方向图 (含互耦 vs 理想) subplot(3,3,3); theta_test 0; a_ideal steering_matrix(theta_test, M, d, lambda); a_distorted C * a_ideal; pattern_ideal abs(fftshift(fft(a_ideal))); pattern_distorted abs(fftshift(fft(a_distorted))); plot(pattern_ideal/max(pattern_ideal), b-, LineWidth, 2); hold on; plot(pattern_distorted/max(pattern_distorted), r--, LineWidth, 2); title(Array Pattern: Ideal vs Coupled); legend(Ideal, With Mutual Coupling); grid on; % DOA估计对比: 未校准 subplot(3,3,4); plot(theta_scan, 10*log10(P_music_uncal/max(P_music_uncal)), k-, LineWidth, 1.5); hold on; plot(theta, zeros(size(theta)), rv, MarkerSize, 10, MarkerFaceColor, r); title(DOA Estimation: Uncalibrated); xlabel(Angle (deg)); ylabel(Power (dB)); grid on; xlim([-90, 90]); % DOA估计对比: 直接求逆校准 subplot(3,3,5); plot(theta_scan, 10*log10(P_music_cal/max(P_music_cal)), b-, LineWidth, 1.5); hold on; plot(theta, zeros(size(theta)), rv, MarkerSize, 10, MarkerFaceColor, r); title(DOA: Direct Inversion); xlabel(Angle (deg)); grid on; xlim([-90, 90]); % DOA估计对比: 快速带状校准 subplot(3,3,6); plot(theta_scan, 10*log10(P_music_fast/max(P_music_fast)), g-, LineWidth, 1.5); hold on; plot(theta, zeros(size(theta)), rv, MarkerSize, 10, MarkerFaceColor, r); title(DOA: Fast Banded Solver); xlabel(Angle (deg)); grid on; xlim([-90, 90]); % 计算时间对比 subplot(3,3,7); bar([time_direct, time_banded]); set(gca, XTickLabel, {Direct Inv, Banded Solve}); ylabel(Time (100 iterations)); title(Computation Time Comparison); set(gca, YScale, log); % 估计误差对比 subplot(3,3,8); [~, peaks_uncal] findpeaks(10*log10(P_music_uncal/max(P_music_uncal)), theta_scan, MinPeakHeight, -10); [~, peaks_cal] findpeaks(10*log10(P_music_cal/max(P_music_cal)), theta_scan, MinPeakHeight, -10); errors_uncal abs(sort(peaks_uncal) - sort(theta)); errors_cal abs(sort(peaks_cal) - sort(theta)); bar([mean(errors_uncal), mean(errors_cal)]); set(gca, XTickLabel, {Uncalibrated, Calibrated}); ylabel(Mean DOA Error (deg)); title(DOA Estimation Accuracy); % 存储需求对比 subplot(3,3,9); mem_full M*M*16; % 复数双精度 mem_banded M*(2*P1)*16; % 仅存储带状区域 bar([mem_full, mem_banded]/1024); set(gca, XTickLabel, {Full Matrix, Banded Storage}); ylabel(Memory (KB)); title(Storage Requirements); sgtitle(Array Mutual Coupling: Fast Inversion Compensation); %% 性能指标输出 fprintf( Mutual Coupling Compensation Performance \n); fprintf(Array Size: %d elements\n, M); fprintf(Coupling Order: %d\n, P); fprintf(Direct Inversion Time: %.4f ms\n, time_direct*10); fprintf(Banded Solver Time: %.4f ms (%.1fx speedup)\n, time_banded*10, time_direct/time_banded); fprintf(Storage Reduction: %.1fx\n, mem_full/mem_banded);8.2 阵列自校准算法阵列自校准(Auto-Calibration)利用接收信号本身估计未知误差参数无需专用校准源或精确已知的参考信号。这类算法将DOA估计与误差校准联合求解形成非线性优化问题。主要方法包括信号子空间拟合(SSF)通过最小化信号子空间与阵列流形的距离联合估计DOA与误差最大似然(ML)方法在统计最优框架下求解基于子空间迭代的交替优化策略将联合问题分解为可处理的子问题。自校准的难点在于参数可辨识性当误差模型过于复杂如同时存在幅相误差、互耦与位置误差某些DOA可能无法唯一确定。可辨识性条件要求信源数、阵列几何与误差结构满足特定约束通常需要至少两个不同方向的信源或阵列具有平移不变结构。8.2.1 信号子空间拟合算法8.2.1.1 交替优化策略交替优化(Alternating Optimization)将联合估计分解为迭代的两步固定误差参数η 估计DOAθ 再固定θ 更新η 。该策略将多维搜索降维至序列一维/低维优化虽非全局最优但计算可行且收敛稳定。DOA估计步骤固定η 使用标准MUSIC或ESPRIT算法处理补偿后的数据快速获得角度估计。误差估计步骤固定θ 构建关于η 的线性最小二乘问题解析求解或低维梯度下降更新。迭代直至收敛参数变化小于阈值或达到最大迭代次数。8.2.1.1.1 信源方向与误差参数的迭代更新迭代更新策略的具体实现依赖于误差模型参数化。对于对角幅相误差固定θ 时γ 的最小二乘解为对于互耦矩阵C 利用Toeplitz约束将问题转化为线性方程组。迭代过程中采用松弛因子防止振荡或使用拟牛顿法加速收敛。初始化对性能至关重要通常先假设无误差执行标准DOA估计将结果作为迭代初值或多重起始点避免局部最优。MATLAB实现信号子空间拟合自校准算法matlab复制%% 脚本: Array_Autocalibration_SSF_Alternating.m %% 内容: 基于信号子空间拟合的阵列自校准算法采用交替优化策略 %% 使用方式: 运行脚本演示联合DOA与幅相误差估计对比交替优化与联合优化性能 %% 优化技术: 1)松弛因子加速收敛 2)子空间跟踪降低特征分解开销 3)多重起始点初始化 clear; clc; close all; %% 阵列与信号参数 M 16; % 阵元数 d 0.5; % 半波长间距 K 2; % 信源数 theta_true [-15, 25]; % 真实DOA SNR 15; % dB N_snap 300; % 快拍数 % 真实幅相误差 gain_error 1 0.2*randn(M, 1); % 增益误差 (20% std) phase_error 0.3*randn(M, 1); % 相位误差 (rad) Gamma_true diag(gain_error .* exp(1j*phase_error)); %% 生成接收信号 A_true exp(-1j*pi*(0:M-1)*sin(theta_true*pi/180)); % 理想流形 S (randn(K, N_snap) 1j*randn(K, N_snap))/sqrt(2); noise (randn(M, N_snap) 1j*randn(M, N_snap))/sqrt(2) * 10^(-SNR/20); X Gamma_true * A_true * S noise; % 样本协方差与特征分解 Rxx X * X / N_snap; [U, D] eig(Rxx); [~, idx] sort(diag(D), descend); Us U(:, idx(1:K)); % 信号子空间 Un U(:, idx(K1:end)); % 噪声子空间 %% 交替优化自校准算法 function [theta_est, Gamma_est, history] alternating_optimization_ssf(X, Us, K, M, max_iter) % 交替优化估计DOA与幅相误差 history.theta zeros(K, max_iter); history.gamma zeros(M, max_iter); history.cost zeros(1, max_iter); % 初始化: 假设无误差执行MUSIC theta_est initial_doa_estimate(X, K, M); Gamma_est eye(M); lambda 0.8; % 松弛因子 for iter 1:max_iter % Step 1: 固定Gamma估计DOA (MUSIC谱峰搜索) X_compensated Gamma_est \ X; % 补偿当前误差 theta_new refine_doa_music(X_compensated, K, M, theta_est); % Step 2: 固定DOA估计Gamma (最小二乘) A exp(-1j*pi*(0:M-1)*sin(theta_new*pi/180)); % 子空间拟合: min ||Us - Gamma*A*T||_F % 对于对角Gamma解析解为: T pinv(A) * Us; % 最小二乘变换矩阵 E Us - A * T; % 残差 % 逐通道优化 (对角约束) gamma_new zeros(M, 1); for m 1:M % 最小化 ||Us(m,:) - gamma(m)*A(m,:)*T||^2 numerator Us(m,:) * (A(m,:)*T); denominator norm(A(m,:)*T)^2; gamma_new(m) numerator / denominator; end % 松弛更新 Gamma_est diag(lambda * abs(gamma_new) (1-lambda) * abs(diag(Gamma_est))) .* ... exp(1j * (lambda * angle(gamma_new) (1-lambda) * angle(diag(Gamma_est)))); % 记录历史 history.theta(:, iter) theta_new; history.gamma(:, iter) diag(Gamma_est); % 计算代价函数 A_new exp(-1j*pi*(0:M-1)*sin(theta_new*pi/180)); T_opt pinv(Gamma_est * A_new) * Us; residual Us - Gamma_est * A_new * T_opt; history.cost(iter) norm(residual, fro); % 收敛检查 if iter 1 abs(history.cost(iter) - history.cost(iter-1)) 1e-6 history.cost history.cost(1:iter); history.theta history.theta(:, 1:iter); history.gamma history.gamma(:, 1:iter); break; end theta_est theta_new; end end %% 辅助函数: 初始DOA估计 function theta_init initial_doa_estimate(X, K, M) R X * X / size(X, 2); [U, ~] eig(R); Un U(:, 1:M-K); theta_scan -90:0.5:90; P_music zeros(size(theta_scan)); for i 1:length(theta_scan) a exp(-1j*pi*(0:M-1)*sin(theta_scan(i)*pi/180)); P_music(i) 1 / abs(a * (Un * Un) * a); end [~, locs] findpeaks(P_music, theta_scan, SortStr, descend, NPeaks, K); theta_init sort(locs); end %% 辅助函数: MUSIC精细估计 function theta_refined refine_doa_music(X, K, M, theta_init) R X * X / size(X, 2); [U, D] eig(R); [~, idx] sort(diag(D), descend); Un U(:, idx(K1:end)); theta_refined zeros(size(theta_init)); for k 1:length(theta_init) % 局部精细搜索 theta_local theta_init(k) (-2:0.05:2); P_local zeros(size(theta_local)); for i 1:length(theta_local) a exp(-1j*pi*(0:M-1)*sin(theta_local(i)*pi/180)); P_local(i) 1 / abs(a * (Un * Un) * a); end [~, max_idx] max(P_local); theta_refined(k) theta_local(max_idx); end theta_refined sort(theta_refined); end %% 执行自校准 [theta_cal, Gamma_cal, history] alternating_optimization_ssf(X, Us, K, M, 50); %% 对比: 未校准MUSIC theta_uncal initial_doa_estimate(X, K, M); %% 对比: 联合优化 (使用fminunc) function cost joint_cost_function(params, Us, M, K) % 联合代价函数 (用于对比) theta params(1:K); gamma_mag params(K1:KM); gamma_phase params(KM1:end); Gamma diag(gamma_mag .* exp(1j*gamma_phase)); A exp(-1j*pi*(0:M-1)*sin(theta*pi/180)); T pinv(Gamma * A) * Us; residual Us - Gamma * A * T; cost norm(residual, fro)^2; end % 联合优化 (耗时较长仅作对比) % theta_init_joint [theta_uncal; abs(diag(Gamma_true)); angle(diag(Gamma_true))]; % options optimoptions(fminunc, Display, off, MaxIterations, 100); % theta_joint_opt fminunc((p) joint_cost_function(p, Us, M, K), theta_init_joint, options); %% 可视化 figure(Position, [50, 50, 1500, 900]); % 收敛曲线 subplot(3,3,1); semilogy(history.cost, b-o, LineWidth, 2); title(Convergence of Alternating Optimization); xlabel(Iteration); ylabel(Subspace Fitting Cost); grid on; % DOA收敛 subplot(3,3,2); plot(history.theta(1,:), b-, LineWidth, 2); hold on; plot(history.theta(2,:), r-, LineWidth, 2); plot(1:size(history.theta,2), theta_true(1)*ones(1,size(history.theta,2)), b--); plot(1:size(history.theta,2), theta_true(2)*ones(1,size(history.theta,2)), r--); title(DOA Convergence); xlabel(Iteration); ylabel(Angle (deg)); legend(Source 1 (est), Source 2 (est), Source 1 (true), Source 2 (true)); grid on; % 幅相误差估计 subplot(3,3,3); plot(1:M, abs(diag(Gamma_true)), b-o, LineWidth, 2); hold on; plot(1:M, abs(diag(Gamma_cal)), r-s, LineWidth, 2); title(Gain Error Estimation); xlabel(Element); ylabel(Gain); legend(True, Estimated); grid on; subplot(3,3,4); plot(1:M, angle(diag(Gamma_true))*180/pi, b-o, LineWidth, 2); hold on; plot(1:M, angle(diag(Gamma_cal))*180/pi, r-s, LineWidth, 2); title(Phase Error Estimation); xlabel(Element); ylabel(Phase (deg)); legend(True, Estimated); grid on; % 误差分布 subplot(3,3,5); gamma_error abs(diag(Gamma_cal) - diag(Gamma_true)); bar(gamma_error); title(Gain Estimation Error); xlabel(Element); ylabel(Error); grid on; % 阵列方向图对比 subplot(3,3,6); w_uncal pinv(Gamma_true * A_true); w_cal pinv(Gamma_cal * exp(-1j*pi*(0:M-1)*sin(theta_cal*pi/180))); theta_scan -90:0.5:90; pattern_uncal zeros(size(theta_scan)); pattern_cal zeros(size(theta_scan)); for i 1:length(theta_scan) a exp(-1j*pi*(0:M-1)*sin(theta_scan(i)*pi/180)); pattern_uncal(i) abs(w_uncal(:,1) * a); pattern_cal(i) abs(w_cal(:,1) * a); end plot(theta_scan, 20*log10(pattern_uncal/max(pattern_uncal)), b-, LineWidth, 1.5); hold on; plot(theta_scan, 20*log10(pattern_cal/max(pattern_cal)), r--, LineWidth, 1.5); title(Beam Pattern Comparison); xlabel(Angle (deg)); ylabel(Gain (dB)); legend(Uncalibrated, Calibrated); grid on; ylim([-40, 0]); % 估计误差统计 subplot(3,3,7); theta_errors [abs(theta_uncal - sort(theta_true)), abs(theta_cal - sort(theta_true))]; bar(theta_errors); set(gca, XTickLabel, {Source 1, Source 2}); ylabel(DOA Error (deg)); legend(Uncalibrated, Calibrated); title(DOA Estimation Error); % CRB下界 (理论最优) subplot(3,3,8); % 计算CRB (简化) SNR_linear 10^(SNR/10); CRB (6/(N_snap*K*SNR_linear)) / (M*(M^2-1)) * (pi/180)^2; % 近似 bar([mean((theta_cal-sort(theta_true)).^2), CRB]); set(gca, XTickLabel, {MSE, CRB}); title(Estimation Error vs CRB); % 算法复杂度分析 subplot(3,3,9); iterations length(history.cost); complexity_per_iter M^3 K*M^2; % 矩阵求逆乘法 bar([M^6, iterations*complexity_per_iter]); % 联合优化 vs 交替优化 set(gca, XTickLabel, {Joint Search, Alternating}); ylabel(Operations); title(Computational Complexity); set(gca, YScale, log); sgtitle(Array Self-Calibration: Signal Subspace Fitting with Alternating Optimization); %% 性能输出 fprintf( Array Self-Calibration Performance \n); fprintf(True DOA: %.2f°, %.2f°\n, theta_true); fprintf(Uncalibrated DOA: %.2f°, %.2f° (Error: %.2f°, %.2f°)\n, ... sort(theta_uncal), abs(sort(theta_uncal) - sort(theta_true))); fprintf(Calibrated DOA: %.2f°, %.2f° (Error: %.2f°, %.2f°)\n, ... theta_cal, abs(theta_cal - sort(theta_true))); fprintf(Convergence iterations: %d\n, length(history.cost)); fprintf(Gain Error RMSE: %.4f\n, rms(abs(diag(Gamma_cal)) - abs(diag(Gamma_true)))); fprintf(Phase Error RMSE: %.4f°\n, rms(angle(diag(Gamma_cal)) - angle(diag(Gamma_true)))*180/pi);8.2.2 辅助源校准辅助源校准(Instrumental Source Calibration)利用空间位置精确已知的辅助信源估计阵列误差参数。相比自校准辅助源方法将非线性联合估计转化为线性或可分离问题具有更高的估计精度与稳定性适用于雷达系统定期维护校准或实时在线校准。单辅助源算法仅需一个已知方向的信源如合作应答机、标校信标通过旋转阵列或改变频率获取多组测量数据构建可观测的误差方程。多辅助源方法利用空间分布的多个信源通过最大似然或加权最小二乘提高估计稳健性。8.2.2.1 单辅助源算法8.2.2.1.1 最小二乘误差估计最小二乘估计最小化测量子空间与模型流形的Frobenius范数距离MATLAB实现单辅助源最小二乘校准%% 脚本: Array_Single_Source_Least_Squares_Calibration.m %% 内容: 基于单辅助源旋转的阵列校准算法使用最小二乘与总体最小二乘估计 %% 使用方式: 运行脚本模拟阵列旋转测量估计幅相误差与互耦参数 %% 优化技术: 1)加权最小二乘优化精度 2)总体最小二乘处理方向误差 3)结构化矩阵约束 clear; clc; close all; %% 阵列参数 M 12; % 阵元数 d 0.5; % 半波长间距 lambda 1; % 真实误差 gain_error_db 3*randn(M, 1); % dB phase_error_deg 10*randn(M, 1); % 度 Gamma_true diag(10.^(gain_error_db/20) .* exp(1j*phase_error_deg*pi/180)); % 互耦 (可选) use_mutual_coupling true; if use_mutual_coupling P 2; % 互耦阶数 c [0.20.1j, 0.05]; C_true toeplitz([1, c, zeros(1, M-length(c)-1)]); M_total C_true * Gamma_true; % 组合误差矩阵 else M_total Gamma_true; end %% 生成旋转测量数据 L 36; % 旋转位置数 (每10度一个位置) SNR 25; % dB theta_source 0; % 辅助源方向 (相对阵列法向) rotations linspace(0, 360-360/L, L); % 阵列旋转角度 X_data zeros(M, L); % 接收信号矩阵 A_matrix zeros(M, L); % 理想流形矩阵 for l 1:L % 旋转后的有效入射角 effective_angle theta_source - rotations(l); % 理想流形 (旋转后) a_l exp(-1j*pi*(0:M-1)*sin(effective_angle*pi/180)); A_matrix(:, l) a_l; % 含误差的接收信号 s_l 1; % 单位信号 noise (randn(M, 1) 1j*randn(M, 1))/sqrt(2) * 10^(-SNR/20); x_l M_total * a_l * s_l noise; X_data(:, l) x_l; end %% 提取信号子空间 (单源情况下即为数据本身) % 多快拍平均提高精度 N_snap_per_pos 50; X_avg zeros(M, L); for l 1:L X_snap repmat(X_data(:,l), 1, N_snap_per_pos) ... 0.1*(randn(M, N_snap_per_pos) 1j*randn(M, N_snap_per_pos)); R_l X_snap * X_snap / N_snap_per_pos; [u, ~] eigs(R_l, 1); X_avg(:, l) u * sqrt(trace(R_l)); % 幅度归一化 end %% 方法1: 标准最小二乘 (LS) function Gamma_ls standard_ls_calibration(U, A) % 最小化 ||U - Gamma*A||_F^2 % 对于对角Gamma: U(:,l) Gamma*a(:,l) U Gamma*A (逐元素) [M, L] size(U); Gamma_ls zeros(M, 1); for m 1:M % 对每个阵元独立求解 % u_m gamma_m * a_m gamma_m u_m / a_m (最小二乘平均) a_m A(m, :).; u_m U(m, :).; % 避免除以零 valid abs(a_m) 0.1; if sum(valid) 0 Gamma_ls(m) sum(u_m(valid) .* conj(a_m(valid))) / sum(abs(a_m(valid)).^2); else Gamma_ls(m) 1; end end Gamma_ls diag(Gamma_ls); end %% 方法2: 加权最小二乘 (WLS) function Gamma_wls weighted_ls_calibration(U, A, weights) % 考虑不同旋转位置的信噪比差异 if nargin 3 weights ones(1, size(U, 2)); end [M, L] size(U); Gamma_wls zeros(M, 1); for m 1:M a_m A(m, :).; u_m U(m, :).; W diag(weights); % 加权最小二乘: gamma (A^H W A)^{-1} A^H W u Gamma_wls(m) (a_m * W * a_m) \ (a_m * W * u_m); end Gamma_wls diag(Gamma_wls); end %% 方法3: 总体最小二乘 (TLS) - 处理方向误差 function Gamma_tls total_ls_calibration(U, A) % TLS考虑A矩阵也存在误差 [M, L] size(U); Gamma_tls zeros(M, 1); for m 1:M % 构造增广矩阵 [a_m, u_m] C_aug [A(m,:)., U(m,:).]; % SVD分解 [~, ~, V] svd(C_aug, econ); % TLS解 v V(:, end); if v(1) ~ 0 Gamma_tls(m) -v(2)/v(1); else Gamma_tls(m) 1; end end Gamma_tls diag(Gamma_tls); end %% 方法4: 互耦幅相联合估计 (结构化LS) function [C_est, Gamma_est] joint_mutual_coupling_calibration(U, A, P) % 估计互耦矩阵(带状Toeplitz)与幅相误差 [M, L] size(U); % 初始化: 先忽略互耦估计Gamma Gamma_temp standard_ls_calibration(U, A); U_compensated Gamma_temp \ U; % 估计互耦系数 (利用子空间拟合) % U ≈ C * A, C为Toeplitz % 构造Toeplitz约束的最小二乘 % 每行是移位版本利用此结构 % 简化: 估计第一行/列 c_est zeros(1, P1); c_est(1) 1; % 归一化 % 利用多旋转位置的几何约束 for p 2:P1 % 第p条对角线估计 row_idx 1:M-p1; col_idx p:M; % U(row_idx, :) ≈ c(p) * A(col_idx, :) % 最小二乘估计c(p) numerator 0; denominator 0; for l 1:L numerator numerator sum(U_compensated(row_idx, l) .* conj(A(col_idx, l))); denominator denominator sum(abs(A(col_idx, l)).^2); end c_est(p) numerator / denominator; end % 构造Toeplitz矩阵 C_est toeplitz([c_est, zeros(1, M-length(c_est))]); % 重新估计Gamma (固定C) U_decoupled C_est \ U; Gamma_est standard_ls_calibration(U_decoupled, A); end %% 执行校准 Gamma_ls standard_ls_calibration(X_avg, A_matrix); % 信噪比权重 (基于接收功率估计) signal_power sum(abs(X_avg).^2, 1); weights signal_power / max(signal_power); Gamma_wls weighted_ls_calibration(X_avg, A_matrix, weights); Gamma_tls total_ls_calibration(X_avg, A_matrix); if use_mutual_coupling [C_est, Gamma_joint] joint_mutual_coupling_calibration(X_avg, A_matrix, P); else C_est eye(M); Gamma_joint Gamma_ls; end %% 评估与可视化 figure(Position, [50, 50, 1500, 900]); % 幅相误差估计对比 subplot(3,3,1); plot(1:M, abs(diag(Gamma_true)), b-o, LineWidth, 2); hold on; plot(1:M, abs(diag(Gamma_ls)), r-s, LineWidth, 1.5); plot(1:M, abs(diag(Gamma_wls)), g-^, LineWidth, 1.5); title(Gain Error Estimation); xlabel(Element); ylabel(Gain); legend(True, LS, WLS); grid on; subplot(3,3,2); plot(1:M, angle(diag(Gamma_true))*180/pi, b-o, LineWidth, 2); hold on; plot(1:M, angle(diag(Gamma_ls))*180/pi, r-s, LineWidth, 1.5); plot(1:M, angle(diag(Gamma_wls))*180/pi, g-^, LineWidth, 1.5); title(Phase Error Estimation); xlabel(Element); ylabel(Phase (deg)); legend(True, LS, WLS); grid on; % 估计误差分布 subplot(3,3,3); gain_error_ls abs(abs(diag(Gamma_ls)) - abs(diag(Gamma_true))); gain_error_wls abs(abs(diag(Gamma_wls)) - abs(diag(Gamma_true))); bar([gain_error_ls, gain_error_wls]); title(Gain Estimation Error); xlabel(Element); legend(LS, WLS); % 校准前后方向图 subplot(3,3,4); theta_test 15; % 测试方向 a_test exp(-1j*pi*(0:M-1)*sin(theta_test*pi/180)); % 未校准 w_uncal ones(M, 1); pattern_uncal abs(w_uncal * (M_total * a_test)); % LS校准 w_ls Gamma_ls \ ones(M, 1); pattern_ls abs(w_ls * (M_total * a_test)); % WLS校准 w_wls Gamma_wls \ ones(M, 1); pattern_wls abs(w_wls * (M_total * a_test)); bar([pattern_uncal, pattern_ls, pattern_wls]); set(gca, XTickLabel, {Uncal, LS, WLS}); title(sprintf(Array Gain at %d°, theta_test)); ylabel(Normalized Gain); % 互耦矩阵估计 (如适用) if use_mutual_coupling subplot(3,3,5); imagesc(abs(C_est)); title(Estimated Mutual Coupling); axis image; colorbar; subplot(3,3,6); imagesc(abs(C_true)); title(True Mutual Coupling); axis image; colorbar; subplot(3,3,7); imagesc(abs(C_est - C_true)); title(Coupling Estimation Error); axis image; colorbar; end % 残差分析 subplot(3,3,8); residual_ls norm(X_avg - Gamma_ls * A_matrix, fro); residual_wls norm((X_avg - Gamma_wls * A_matrix)*diag(sqrt(weights)), fro); bar([residual_ls, residual_wls]); set(gca, XTickLabel, {LS, WLS}); title(Fitting Residual); ylabel(Frobenius Norm); % CRB分析 (理论界) subplot(3,3,9); SNR_linear 10^(SNR/10); var_gain 1/(2*N_snap_per_pos*SNR_linear); var_phase var_gain; % 近似 crb_vals sqrt([var_gain, var_phase]); actual_error [rms(abs(diag(Gamma_wls)) - abs(diag(Gamma_true))), ... rms(angle(diag(Gamma_wls)) - angle(diag(Gamma_true)))]; bar([actual_error; crb_vals]); set(gca, XTickLabel, {Actual, CRB}); legend(Gain, Phase); title(Estimation Error vs CRB); sgtitle(Single Source Array Calibration: LS/WLS/TLS Estimation); %% 性能输出 fprintf( Single Source Calibration Performance \n); fprintf(Rotation positions: %d\n, L); fprintf(Snapshots per position: %d\n, N_snap_per_pos); fprintf(Gain Error RMSE (LS): %.4f dB\n, 20*log10(rms(gain_error_ls))); fprintf(Gain Error RMSE (WLS): %.4f dB\n, 20*log10(rms(gain_error_wls))); fprintf(Phase Error RMSE (LS): %.2f°\n, rms(angle(diag(Gamma_ls)) - angle(diag(Gamma_true)))*180/pi); fprintf(Phase Error RMSE (WLS): %.2f°\n, rms(angle(diag(Gamma_wls)) - angle(diag(Gamma_true)))*180/pi); if use_mutual_coupling fprintf(Coupling Matrix RMSE: %.4f\n, rms(C_est - C_true, all)); end第八章小结本章系统阐述了阵列校准与误差补偿的理论与算法建立了幅相误差与互耦效应的数学模型提出了基于结构化矩阵求逆的快速补偿方法。自校准算法通过信号子空间拟合与交替优化实现联合DOA-误差估计辅助源方法利用最小二乘框架实现高精度离线校准。这些技术为实际雷达系统的误差抑制与性能优化提供了工程实现路径。