【雷达信号优化】第E章 分布式与网络化雷达（系统级扩展）
·
目录
E.1.1.2 多基地RCS闪烁与NLOS（非视距）路径利用
E.2.1.1 集中式 vs 分布式粒子滤波（Consensus-based Tracking）
E.2.1.2 通信带宽受限下的量化检测与censoring策略
脚本7：集中式 vs 分布式粒子滤波（Consensus Tracking）
E.1 多静态雷达信号模型
E.1.1 双/多基地几何与信号特性
E.1.1.1 椭圆定位与多普勒椭球几何
多静态雷达系统通过分离布置的发射站与接收站构建立体探测网络，其目标定位几何由双基地椭圆方程组描述。对于第$m$个发射站与第$n$个接收站构成的发射-接收对，目标存在时的距离和约束定义了三维空间中的椭球面，其焦点位于发射站与接收站坐标。该椭球面的数学表达为发射-接收距离和等于常数的点的集合，目标位置必须满足所有参与观测的发射-接收对所定义的椭球约束交集。
双基地距离和$R_{mn} = R_{Tm} + R_{Rn}$决定椭圆面的主轴长度，其中$R_{Tm}$为目标到第$m$个发射站的距离，$R_{Rn}$为目标到第$n$个接收站的距离。椭圆偏心率由发射-接收基线长度$L_{mn}$与距离和之比确定，基线越长则椭圆越扁长，距离分辨率在基线方向与垂直方向呈现各向异性。多基地定位通过多个椭球面的三维交集确定目标位置，当发射站与接收站几何配置适当时，交集退化为空间单点实现无模糊定位。 多普勒频移在多基地几何中呈现复杂的方位依赖性。目标速度矢量$\mathbf{v}$在双基地方向上的投影决定多普勒频移大小，其表达式包含发射方向与接收方向的单位矢量组合。对于运动目标，多普勒频移由双基地角$\beta$与速度方向夹角共同决定，形成多普勒椭球面——等频移点的轨迹构成以发射-接收焦点为基准的旋转椭球面。多普勒分辨能力取决于速度矢量在双基地角平分线方向的投影精度，宽双基地角配置增强速度分辨率但降低信号能量积累效率。
椭圆定位的精度由几何精度因子表征，该矩阵依赖于发射-接收站的布局几何与测量误差协方差。当发射站与接收站均匀分布于目标周围时，几何精度因子各向同性，定位误差在三维空间呈球形分布；当站点集中于某一方位时，误差椭球沿该方向拉伸形成定位盲区。优化布站策略通过最小化几何精度因子的迹或行列式，在约束条件下寻求费雪信息矩阵的最大化。
$$(R_{Tm} + R_{Rn})^2 = \|\mathbf{p} - \mathbf{p}_{Tm}\|^2 + \|\mathbf{p} - \mathbf{p}_{Rn}\|^2 + 2\|\mathbf{p} - \mathbf{p}_{Tm}\|\|\mathbf{p} - \mathbf{p}_{Rn}\|\cos\beta_{mn}$$
其中$\mathbf{p}$为目标位置矢量，$\mathbf{p}{Tm}, \mathbf{p}{Rn}$分别为发射与接收站位置，$\beta_{mn}$为双基地角。
E.1.1.2 多基地RCS闪烁与NLOS（非视距）路径利用
多基地雷达目标截面积呈现显著的视角依赖性与闪烁特性。复杂目标在不同双基地角照射下，其雷达截面积随视角快速起伏，闪烁频率与目标姿态变化速率及波长相关。多基地配置通过空间分集平滑雷达截面积起伏，当各通道的雷达截面积统计不相关时，非相干积累获得分集增益，检测概率曲线陡峭化。
双基地雷达截面积建模需考虑目标几何结构与材料特性。光学区目标适用物理光学近似，雷达截面积与表面曲率半径和目标投影面积相关；谐振区目标需采用矩量法或有限元法全波仿真。多基地配置下，前向散射增强效应在双基地角接近180度时显著，此时目标阴影区衍射贡献主导，雷达截面积可达光学区数值的十倍以上，适用于隐身目标的探测与跟踪。
非视距路径在多基地网络中既是干扰源也是潜在信息载体。城市峡谷与室内环境中，直射路径被遮挡，多径反射信号成为主要能量来源。传统雷达将多径视为干扰加以抑制，而NLOS利用技术主动挖掘反射路径中的目标信息。通过构建多径传播图模型，将墙面、地面反射视为虚拟发射源，扩展有效发射-接收对数量。多径参数估计联合优化目标位置与反射面几何，利用到达时间差与到达角度实现镜像定位。
多径利用的精度受限于反射面位置标定误差与散射系数不确定性。当反射面粗糙度满足基尔霍夫近似条件时，镜面反射主导，路径可预测；当粗糙度增加导致漫散射显著时，信号能量分散降低信噪比。多基地网络通过冗余路径测量补偿单一路径的不确定性，利用一致性检验识别并剔除异常多径，提升NLOS场景下的定位鲁棒性。
$$\sigma_{bi}(\beta, \theta) = \sigma_{mono}(\theta_i) \cdot \frac{\cos^2\theta_i}{\cos^4(\beta/2)} \cdot D(\beta)$$
其中$\sigma_{bi}$为双基地雷达截面积，$\theta_i$为入射角，$D(\beta)$为方向性函数，前向散射增强时$D(\pi) \gg 1$。
E.1.1.3 直达波与目标回波的时间同步
多基地雷达的相干处理严格依赖于发射站与接收站间的时间同步精度。直达波信号作为时间参考基准，提供发射时间戳与载波相位基准，其提取质量直接影响目标回波的相关处理增益。直达波路径与目标回波路径的传播时延差异导致时间窗口偏移，需通过预测模型或实时校准消除。
时间同步误差分为系统误差与随机误差两类。系统误差源于时钟频偏与固定传播延迟，可通过双向时间传递与校准报文补偿；随机误差源于相位噪声与传播环境起伏，需通过锁相环跟踪与平滑滤波抑制。在分布式相参阵列中，时间同步精度需达到亚波长尺度，即皮秒级时间精度或毫米级空间精度，以保证多通道信号的相位一致性。
直达波提取算法在强直达波与弱直达波场景采用不同策略。强直达波场景通过波束形成在主瓣内提取参考信号；弱直达波场景采用相关检测与积累提升信噪比，或利用导航卫星信号作为外部时间基准。直达波与目标回波的互相关处理提供双基地距离和估计，其精度由信号带宽与信噪比决定，时延分辨率为带宽倒数的量级。
时钟同步网络采用主从架构或分布式共识架构。主从架构以高精度原子钟为主节点，分发时间基准至从节点，单点失效风险高；分布式共识架构通过节点间双向比对实现时间一致性，鲁棒性强但收敛速度慢。无线传感器网络时间同步协议在能耗与精度间权衡，采用间歇同步与漂移补偿降低通信开销。
$$\Delta t = \frac{\Delta\phi}{2\pi f_c} + nT_s$$
其中$\Delta t$为时间同步误差，$\Delta\phi$为相位测量误差，$f_c$为载波频率，$n$为整数周期模糊数，$T_s$为采样周期。
E.1.2 分布式信号级融合
E.1.2.1 相参积累的相位同步误差补偿
分布式相参雷达通过多接收站信号的相位相干叠加，实现虚拟孔径扩展与信噪比增益。相参积累增益理论上与接收站数量的平方成正比，前提是各通道信号相位精确对齐。相位同步误差导致积累增益退化，均方根相位误差$\sigma_\phi$引起的增益损失近似为$\exp(-\sigma_\phi^2/2)$，当误差超过60度时相参增益急剧下降。 相位误差来源包括时钟相位噪声、传播路径差异与目标视角变化。时钟相位噪声在短基线场景可通过共视参考源抑制；长基线场景需采用独立高稳时钟与事后相关处理。传播路径差异由大气折射率起伏引起，可通过气象数据模型补偿或利用直达波实时校准。目标视角变化导致各接收站回波的多普勒相位差异，需在积累前进行运动补偿。 相位同步补偿算法在频域与时域实现。频域方法基于互相关相位估计，在窄带假设下计算通道间相位差；时域方法基于自适应滤波，通过最小均方误差准则调整通道相位权重。迭代投影算法交替估计目标参数与通道相位误差，联合优化实现自校准。基于最大似然的联合估计将相位误差视为nuisance参数，通过边缘化或联合求解消除其影响。 相参积累的性能边界由克拉美罗界表征，该边界依赖于相位误差统计与信号模型。在随机相位误差场景，贝叶斯相参积累利用相位先验分布优化加权系数；在确定性相位误差场景，校准辅助积累通过导频信号实时跟踪相位漂移。分布式相参的稳健性设计预留相位误差容限，在允许一定增益损失的前提下降低系统复杂度。 $$G_{coh} = \left|\sum_{n=1}^{N} e^{j\phi_n}\right|^2 = N^2 \cdot \left|\frac{1}{N}\sum_{n=1}^{N} e^{j\phi_n}\right|^2 \approx N^2 e^{-\sigma_\phi^2}$$ 其中$G_{coh}$为相参积累增益，$N$为接收站数量，$\phi_n$为第$n$个通道的相位误差。
E.1.2.2 非相参积累的GLRT推导与CFAR检测
非相参积累在相位同步困难或相位误差过大场景提供稳健的性能保障。广义似然比检验框架处理多基地雷达检测问题，假设检验针对各接收站信号的联合分布建立。零假设下各通道仅含噪声，备择假设下各通道含幅度未知的信号分量，相位在$[0,2\pi)$均匀分布或视为未知参数。 对数似然比在各通道统计独立假设下分解为单通道似然比之和。对于高斯噪声背景，平方律检波器输出服从指数分布（零假设）与非中心卡方分布（备择假设）。非相参积累对各通道检波输出求和或求极大值，求和策略适用于高信噪比场景，极大值策略适用于低信噪比与信号闪烁场景。 恒虚警率检测在非相参积累中需考虑多通道参考单元的选择策略。集中式CFAR利用所有通道的参考单元联合估计背景杂波水平，提高估计精度但增加通信开销；分布式CFAR各通道独立估计后融合，降低复杂度但损失估计效率。有序统计CFAR在多通道场景对参考单元排序后选取第$k$大值作为阈值，对异常值鲁棒。 检测概率与虚警概率的解析表达式涉及Marcum Q函数或广义卡方分布，多通道场景计算复杂。高信噪比近似采用高斯近似简化分析，低信噪比场景采用级数展开或数值积分。非相参积累的性能损失相对于理想相参积累为$N$倍（功率积累）相对于$N^2$倍（电压积累），但在相位失配严重时实际性能可能优于相参方案。 $$\Lambda_{GLRT} = \sum_{n=1}^{N} \ln I_0\left(\frac{|z_n|}{\sigma_n^2}\sqrt{2P_n}\right) \approx \sum_{n=1}^{N} \frac{P_n}{\sigma_n^2} |z_n|^2$$ 其中$I_0$为修正贝塞尔函数，$P_n$为第$n$个通道的信号功率，$z_n$为匹配滤波输出，$\sigma_n^2$为噪声方差，高信噪比近似下退化为平方律积累。
E.1.2.3 分布式压缩感知与稀疏恢复（多视角联合重构）
分布式压缩感知理论突破单传感器奈奎斯特采样限制，利用多传感器联合稀疏性降低采样率。多基地雷达场景中，目标在三维空间稀疏分布，各接收站观测的是同一目标集合在不同视角下的投影。联合稀疏模型要求各通道信号共享相同的支撑集，但允许不同的系数值。
多视角联合重构通过协作稀疏恢复提升分辨率与鲁棒性。同时正交匹配追踪算法扩展单通道贪婪算法至多通道，在每次迭代中选择对所有通道残差贡献最大的原子。分布式基追踪算法采用交替方向乘子法将全局优化问题分解为本地子问题与一致性约束，通过双上升迭代实现分布式求解。共识优化框架在节点间交换本地估计的软信息，通过平均共识达成全局一致解。
感知矩阵的块对角结构反映分布式观测特性，各块对应一个接收站的观测矩阵。互块相关性影响恢复性能，当各站视角差异大时感知矩阵列间相关性低，受限等距性质易满足。联合恢复的性能边界优于单通道独立恢复，所需测量数随通道数增加而减少，节省率取决于目标场景的相关性结构。
模型失配与阵列误差在分布式场景更为复杂，各通道可能存在不同的阵列流形误差。联合校准算法将通道增益与相位误差纳入稀疏恢复框架，通过扩展字典原子实现自校准。贝叶斯分布式压缩感知引入误差先验分布，通过变分推断近似后验，提供恢复结果的不确定性量化。
$$\min_{\mathbf{X}} \|\mathbf{Y} - \mathbf{A}\mathbf{X}\|_F^2 + \lambda\|\mathbf{X}\|_{2,1}$$
其中$\mathbf{Y}$为多通道观测矩阵，$\mathbf{A}$为块对角感知矩阵，$\mathbf{X}$为行稀疏信号矩阵，$|\cdot|_{2,1}$为行范数诱导的联合稀疏正则化。
E.2 传感器网络协同处理
E.2.1 网络拓扑与数据融合架构
E.2.1.1 集中式 vs 分布式粒子滤波（Consensus-based Tracking）
目标跟踪在多传感器网络中面临数据融合架构的选择。集中式粒子滤波将全部传感器的测量数据汇集至融合中心，利用完整信息计算后验分布，性能最优但通信负担重、鲁棒性差。分布式粒子滤波在本地节点维持粒子集合，通过有限的邻域通信实现全局后验的近似，适用于带宽受限与无中心架构。
共识粒子滤波通过迭代平均实现全局似然的分布式计算。各节点基于本地测量计算权重，通过共识协议与邻域节点交换权重信息，迭代收敛至全局似然的加权平均。gossip算法在随机连接网络中通过异步平均达成共识，收敛速度与网络代数连通度相关。分布式重采样面临粒子枯竭与副本问题，采用全局随机数生成或分布式决议确保各节点粒子集合的一致性。
混合架构结合集中与分布的优势，分层粒子滤波在簇头节点执行局部集中处理，簇间通过分布式共识协调。自适应架构根据网络状态动态切换工作模式，拓扑稳定时采用集中式保证精度，拓扑剧变时切换至分布式保证连续性。网络编码技术将测量信息编码传输，提高带宽效率与容错能力。
粒子滤波的可扩展性通过边际化与分解实现。Rao-Blackwellised粒子滤波将线性状态用卡尔曼滤波解析处理，仅对非线性状态采样，降低粒子数需求。分布式实现中，各节点对局部线性状态边际化，仅交换非线性状态的粒子云，减少通信量。分布式粒子滤波的收敛性由网络连通性与共识步长保证，有限次数共识迭代引入的偏差可通过自适应补偿修正。
$$\mathbf{x}_k^{(i)} \sim \sum_{j=1}^{N} \tilde{w}_{k-1}^{(j)} p(\mathbf{x}_k | \mathbf{x}_{k-1}^{(j)}), \quad \tilde{w}_k^{(i)} \propto p(\mathbf{z}_k | \mathbf{x}_k^{(i)})$$
其中共识步骤确保各节点的归一化权重$\tilde{w}_k^{(i)}$通过邻域通信收敛至全局似然平均。
E.2.1.2 通信带宽受限下的量化检测与censoring策略
传感器网络通信带宽受限要求对传输信息进行压缩与选择。量化检测将本地软信息或硬判决映射至有限比特，降低传输开销但引入量化误差。最优量化器设计在检测性能与通信速率间权衡，基于广义率失真理论或检测误差最小化准则。独立量化各传感器观测忽略传感器间相关性，分布式信源编码利用相关结构实现压缩增益。
Censoring策略仅在本地检测统计量超越阈值时触发传输，沉默时发送端节省能量与带宽。阈值设计平衡检测性能与通信负载，高阈值降低传输率但增加漏警，低阈值提高检测概率但浪费资源。自适应censoring根据信道状态与目标先验动态调整阈值，信道恶劣时提高阈值避免无效传输，目标不确定度高时降低阈值保证跟踪连续性。
压缩感知理论指导稀疏事件的低带宽检测。当目标稀疏存在时，各传感器本地投影观测至随机向量，传输压缩后的测量值。融合中心通过联合稀疏恢复重构目标状态，实现低于奈奎斯特速率的检测。二值量化极端压缩至单比特，符号型压缩感知理论证明在适当条件下单比特测量足以支持稀疏恢复，适用于极低带宽的物联网雷达。
联合优化检测与传输策略在跨层设计中实现。物理层波形设计考虑MAC层碰撞，MAC层协议优化考虑检测性能。博弈论框架建模各传感器的传输竞争，纳什均衡策略指导分布式接入，避免碰撞与拥塞。强化学习动态适应流量模式，在突发检测需求时自动提升传输优先级。
$$\max_{\mathcal{Q}, \mathcal{T}} P_D(\mathcal{Q}, \mathcal{T}) \quad \text{s.t.} \quad R(\mathcal{Q}, \mathcal{T}) \le R_{max}, \quad P_{FA} \le \alpha$$
其中$\mathcal{Q}$为量化策略集合，$\mathcal{T}$为censoring阈值集合，$R$为平均通信速率，$R_{max}$为带宽约束。
E.2.2 网络鲁棒性与抗毁性
E.2.2.1 节点失效下的波形重分配与覆盖补偿
传感器网络的节点失效由硬件故障、能量耗尽或外部攻击引起，破坏网络覆盖与连通性。鲁棒网络设计预留冗余节点与备用链路，通过拓扑重构维持功能。波形重分配在节点失效后调整剩余节点的发射参数，补偿覆盖盲区并维持探测性能。
覆盖补偿算法识别失效节点负责的监视区域，计算覆盖缺口的几何与频率特性。邻近节点通过功率增加与波束扩展接管失效区域，优化问题在功率约束与干扰约束下最小化覆盖损失。多输入多输出波形设计调整发射方向图，在失效节点方向形成波束零陷避免能量浪费，在补偿方向增强增益。移动节点重定位通过势能场或凸优化方法调整网络几何，快速填补覆盖缺口。
频谱重分配在节点失效后避免频谱空洞与干扰加剧。失效节点占用的频段释放后，邻近节点扩展带宽提升分辨率，或新接入节点使用该频段维持容量。图着色算法指导频谱重分配，确保邻域节点正交性同时最大化频谱复用。博弈论频谱接入在节点失效后重新计算纳什均衡，动态调整策略适应新拓扑。
网络性能退化评估量化节点失效的影响。连通度保持概率衡量随机失效下网络保持连通的概率，$k$-连通网络可容忍$k-1$个节点失效。覆盖保持概率衡量区域覆盖比例，与节点密度与失效概率相关。弹性指标定义为性能退化曲线下面积，综合评估网络在级联失效中的承受能力。
$$\min_{\mathbf{w}_i} \sum_{j \in L} \left(P_d(\mathbf{p}_j) - P_d^{target}\right)^2 + \lambda \sum_{i \in A} \|\mathbf{w}_i\|^2$$
其中$L$为覆盖缺口位置集合，$A$为活跃节点集合，$P_d$为检测概率，$\mathbf{w}_i$为第$i$个节点的波形权矢量。
E.2.2.2 对抗网络干扰的分布式一致性算法
对抗环境下传感器网络面临拜占庭故障与恶意干扰，部分节点发送错误信息或故意破坏共识。传统一致性算法在存在恶意节点时收敛至错误值，需采用鲁棒一致性或弹性分布式算法。鲁棒统计量如中值、截断均值替代均值共识，限制异常值的影响。弹性共识算法要求网络连通度超过恶意节点数的两倍，确保正常节点子图保持连通。
拜占庭容错算法通过多数表决与交叉验证识别恶意节点。各节点收集邻域信息，若某节点报告值与本地估计偏差过大则标记为可疑，通过多轮通信确认故障。分布式故障检测利用观测一致性检验，统计异常报告频率，隔离高频异常节点。网络自清洁机制在检测到恶意节点后重构拓扑，切断与恶意节点的连接，维持剩余网络的连通性。
对抗干扰的物理层技术结合波形设计与网络协议。跳频与波形分集使干扰难以瞄准，认知无线电动态避开干扰频段。协作波束形成将多个节点的发射能量相干叠加至目标方向，同时在干扰方向形成零陷。物理层安全利用信道特性实现密钥生成与安全传输，防止窃听与篡改。
分布式优化在对抗环境下采用安全多智能体算法。梯度下降算法结合梯度裁剪与验证，确保更新方向不被恶意梯度误导。交替方向乘子法的变体在增广拉格朗日函数中加入鲁棒惩罚项，抑制异常乘子。区块链技术为分布式雷达网络提供不可篡改的数据记录与共识机制，确保测量历史与融合结果的可信性。
$$x_i(t+1) = x_i(t) + \epsilon \sum_{j \in N_i} a_{ij} \cdot \psi(x_j(t) - x_i(t))$$
其中$\psi(\cdot)$为鲁棒非线性函数（如符号函数、饱和函数），抑制异常差值$x_j(t) - x_i(t)$的影响，$N_i$为节点$i$的邻域集合。
第二部分 代码实现
脚本1：椭圆定位与多普勒椭球几何分析
内容说明：本脚本实现多基地雷达的双基地椭圆定位几何建模，计算并可视化椭球约束面、多普勒频移等值面，以及多站联合定位的精度几何因子分布。支持任意部署的发射-接收站配置。
使用方式：直接运行main_Elliptic_Geometry，配置txPositions（发射站坐标）、rxPositions（接收站坐标）与targetGrid（目标搜索网格），系统将生成三维椭球定位曲面与GDOP分布图。
matlab
复制
%% 脚本1：椭圆定位与多普勒椭球几何分析
% 功能：多基地双基地椭圆定位，多普勒椭球几何，GDOP分析
% 使用：配置txPositions(Ntx x 3)、rxPositions(Nrx x 3)、targetGrid范围
% 输出：椭球约束面、多普勒等值面、定位精度分布、双基地角分布
function main_Elliptic_Geometry()
clear; close all; clc;
%% 网络配置（三维坐标，单位：km）
txPos = [0, 0, 0; % 发射站1
50, 0, 0; % 发射站2
25, 43, 0]; % 发射站3（等边三角形）
rxPos = [100, 0, 0; % 接收站1
50, 86, 0; % 接收站2
0, 86, 0; % 接收站3
50, 43, 30]; % 接收站4（空中）
numTx = size(txPos, 1);
numRx = size(rxPos, 1);
%% 生成观测（模拟）
trueTarget = [45, 35, 10]; % 真实目标位置 [x,y,z] km
c = 3e5; % 光速 km/s
wavelength = 0.03; % 波长 30cm
% 计算真实距离和与多普勒
[trueRanges, trueDopplers, bistaticAngles] = calculateBistaticObservations(...
trueTarget, txPos, rxPos, c, wavelength, [50, 30, 0]); % 目标速度50km/s向北，30向东
fprintf('真实目标位置: [%.1f, %.1f, %.1f] km\n', trueTarget);
fprintf('距离和观测:\n');
disp(trueRanges);
%% 定位求解（解析法与数值法）
% 最小二乘定位
estPosition = ellipticLocalization(txPos, rxPos, trueRanges, [50, 50, 5]);
fprintf('估计位置: [%.2f, %.2f, %.2f] km (误差: %.2f km)\n', ...
estPosition, norm(estPosition - trueTarget));
%% 可视化
figure('Name', '多基地椭圆定位几何', 'Position', [50 50 1400 900]);
% 3D几何配置
subplot(2, 3, 1);
scatter3(txPos(:,1), txPos(:,2), txPos(:,3), 300, 'r^', 'filled'); hold on;
scatter3(rxPos(:,1), rxPos(:,2), rxPos(:,3), 300, 'bo', 'filled');
scatter3(trueTarget(1), trueTarget(2), trueTarget(3), 200, 'g*', 'LineWidth', 3);
% 绘制基线
for i = 1:numTx
for j = 1:numRx
plot3([txPos(i,1), rxPos(j,1)], [txPos(i,2), rxPos(j,2)], ...
[txPos(i,3), rxPos(j,3)], 'k:', 'LineWidth', 0.5);
end
end
title('多基地几何配置');
xlabel('X (km)'); ylabel('Y (km)'); zlabel('Z (km)');
legend('发射站', '接收站', '目标'); grid on; axis equal;
% 椭球约束面（选择第一对Tx-Rx展示）
subplot(2, 3, 2);
plotEllipsoid(txPos(1,:), rxPos(1,:), trueRanges(1), 50);
title(sprintf('双基地椭球面 (Tx1-Rx1, R=%.1f km)', trueRanges(1)));
xlabel('X (km)'); ylabel('Y (km)'); zlabel('Z (km)');
grid on; axis equal;
% 多普勒椭球（等频移面）
subplot(2, 3, 3);
plotDopplerEllipsoid(txPos(1,:), rxPos(1,:), trueDopplers(1), wavelength, 50);
title(sprintf('多普勒等值面 (fd=%.1f Hz)', trueDopplers(1)));
xlabel('X (km)'); ylabel('Y (km)'); zlabel('Z (km)');
grid on; axis equal;
% GDOP分布（水平面切片）
subplot(2, 3, 4);
gdopMap = calculateGDOPMap(txPos, rxPos, [0:2:100], [0:2:100], 10, c);
imagesc([0 100], [0 100], gdopMap);
title('水平面定位精度几何因子 (GDOP)');
xlabel('X (km)'); ylabel('Y (km)');
colorbar; colormap jet; axis xy;
hold on;
scatter3(txPos(:,1), txPos(:,2), txPos(:,3), 100, 'r^', 'filled');
scatter3(rxPos(:,1), rxPos(:,2), rxPos(:,3), 100, 'bo', 'filled');
% 双基地角分布
subplot(2, 3, 5);
bistaticGrid = calculateBistaticAngleMap(txPos, rxPos, [0:5:100], [0:5:100], 10);
imagesc([0 100], [0 100], bistaticGrid);
title('双基地角分布 (度)');
xlabel('X (km)'); ylabel('Y (km)');
colorbar; colormap hsv; axis xy;
% 定位误差椭圆（三维投影）
subplot(2, 3, 6);
errorEllipsoid = calculateErrorEllipsoid(txPos, rxPos, trueTarget, c, 0.1); % 0.1 km测距误差
plotErrorEllipsoid(errorEllipsoid, trueTarget);
title('定位误差椭球 (3-sigma)');
xlabel('X (km)'); ylabel('Y (km)'); zlabel('Z (km)');
grid on; axis equal;
%% 精度分析
fprintf('\n=== 几何精度分析 ===\n');
fprintf('平均双基地角: %.1f°\n', mean(bistaticAngles) * 180/pi);
fprintf('GDOP at target: %.2f\n', calculateGDOP(txPos, rxPos, trueTarget, c));
fprintf('水平定位精度: ±%.2f km (假设测距误差0.1km)\n', ...
calculateGDOP(txPos, rxPos, trueTarget, c) * 0.1);
end
%% 计算双基地观测
function [ranges, dopplers, angles] = calculateBistaticObservations(target, txPos, rxPos, c, lambda, velocity)
numTx = size(txPos, 1);
numRx = size(rxPos, 1);
ranges = zeros(numTx, numRx);
dopplers = zeros(numTx, numRx);
angles = zeros(numTx, numRx);
for i = 1:numTx
for j = 1:numRx
R1 = norm(target - txPos(i,:));
R2 = norm(target - rxPos(j,:));
ranges(i,j) = R1 + R2;
% 双基地角
vec1 = (target - txPos(i,:)) / R1;
vec2 = (target - rxPos(j,:)) / R2;
angles(i,j) = acos(dot(vec1, vec2));
% 多普勒（假设速度矢量）
vTx = dot(velocity, vec1);
vRx = dot(velocity, vec2);
dopplers(i,j) = -(vTx + vRx) / lambda; % 负号表示接近
end
end
end
%% 椭圆定位求解（非线性最小二乘）
function position = ellipticLocalization(txPos, rxPos, measurements, initialGuess)
% 使用Levenberg-Marquardt算法求解
options = optimset('Display', 'off', 'TolFun', 1e-10);
costFunc = @(p) reshape(arrayfun(@(i) norm(p-txPos(i,:)) + ...
arrayfun(@(j) norm(p-rxPos(j,:)), 1:size(rxPos,1)) - ...
measurements(i,:)', 1:size(txPos,1)), [], 1);
position = lsqnonlin(costFunc, initialGuess, [], [], options);
end
%% 绘制椭球面
function plotEllipsoid(p1, p2, sumDist, gridRes)
% 椭球参数：焦点p1, p2，主轴长度sumDist
center = (p1 + p2) / 2;
focalDist = norm(p1 - p2);
a = sumDist / 2; % 半长轴
c = focalDist / 2;
b = sqrt(a^2 - c^2); % 半短轴
% 旋转矩阵（对准焦点连线）
direction = (p2 - p1) / focalDist;
theta = atan2(direction(2), direction(1));
phi = acos(direction(3));
[u, v] = meshgrid(linspace(0, 2*pi, gridRes), linspace(0, pi, gridRes/2));
x = a * cos(u) .* sin(v);
y = b * sin(u) .* sin(v);
z = b * cos(v);
% 旋转并平移
R = [cos(theta)*cos(phi), -sin(theta), cos(theta)*sin(phi);
sin(theta)*cos(phi), cos(theta), sin(theta)*sin(phi);
-sin(phi), 0, cos(phi)];
coords = R * [x(:)'; y(:)'; z(:)'];
x = reshape(coords(1,:), size(u)) + center(1);
y = reshape(coords(2,:), size(u)) + center(2);
z = reshape(coords(3,:), size(u)) + center(3);
surf(x, y, z, 'FaceAlpha', 0.3, 'EdgeColor', 'none', 'FaceColor', 'c');
hold on;
scatter3(p1(1), p1(2), p1(3), 200, 'r^', 'filled');
scatter3(p2(1), p2(2), p2(3), 200, 'bo', 'filled');
end
%% 绘制多普勒椭球（简化：旋转双曲面近似）
function plotDopplerEllipsoid(p1, p2, fd, lambda, gridRes)
% 等频移面为旋转双曲面（简化展示）
center = (p1 + p2) / 2;
[x, y] = meshgrid(linspace(center(1)-20, center(1)+20, gridRes), ...
linspace(center(2)-20, center(2)+20, gridRes));
z = center(3) * ones(size(x)) + fd * 0.1; % 简化可视化
surf(x, y, z, 'FaceAlpha', 0.3, 'EdgeColor', 'none', 'FaceColor', 'm');
hold on;
scatter3(p1(1), p1(2), p1(3), 200, 'r^', 'filled');
scatter3(p2(1), p2(2), p2(3), 200, 'bo', 'filled');
end
%% 计算GDOP地图
function gdopMap = calculateGDOPMap(txPos, rxPos, xGrid, yGrid, zFixed, c)
gdopMap = zeros(length(yGrid), length(xGrid));
for i = 1:length(yGrid)
for j = 1:length(xGrid)
target = [xGrid(j), yGrid(i), zFixed];
gdopMap(i, j) = calculateGDOP(txPos, rxPos, target, c);
end
end
end
%% 计算单点GDOP
function gdop = calculateGDOP(txPos, rxPos, target, c)
numTx = size(txPos, 1);
numRx = size(rxPos, 1);
% 构建几何矩阵
G = zeros(numTx * numRx, 3);
idx = 1;
for i = 1:numTx
for j = 1:numRx
vec1 = (target - txPos(i,:)) / norm(target - txPos(i,:));
vec2 = (target - rxPos(j,:)) / norm(target - rxPos(j,:));
G(idx, :) = (vec1 + vec2) / c;
idx = idx + 1;
end
end
% GDOP = sqrt(trace(inv(G'G)))
try
gdop = sqrt(trace(inv(G' * G)));
catch
gdop = inf;
end
end
%% 计算双基地角地图
function angleMap = calculateBistaticAngleMap(txPos, rxPos, xGrid, yGrid, zFixed)
angleMap = zeros(length(yGrid), length(xGrid));
for i = 1:length(yGrid)
for j = 1:length(xGrid)
target = [xGrid(j), yGrid(i), zFixed];
% 计算平均双基地角
angles = [];
for ti = 1:size(txPos,1)
for rj = 1:size(rxPos,1)
v1 = target - txPos(ti,:);
v2 = target - rxPos(rj,:);
angles = [angles, acos(dot(v1,v2)/(norm(v1)*norm(v2)))];
end
end
angleMap(i, j) = mean(angles) * 180 / pi;
end
end
end
%% 计算误差椭球
function ellipsoid = calculateErrorEllipsoid(txPos, rxPos, target, c, rangeError)
G = [];
for i = 1:size(txPos,1)
for j = 1:size(rxPos,1)
vec1 = (target - txPos(i,:)) / norm(target - txPos(i,:));
vec2 = (target - rxPos(j,:)) / norm(target - rxPos(j,:));
G = [G; (vec1 + vec2) / c];
end
end
C = inv(G' * G) * rangeError^2; % 协方差矩阵
[V, D] = eig(C);
ellipsoid.axes = 3 * sqrt(diag(D)); % 3-sigma轴长
ellipsoid.rotation = V;
ellipsoid.center = target;
end
%% 绘制误差椭球
function plotErrorEllipsoid(ell, target)
[x, y, z] = ellipsoid(0, 0, 0, ell.axes(1), ell.axes(2), ell.axes(3), 30);
% 应用旋转和平移
coords = ell.rotation * [x(:)'; y(:)'; z(:)'];
x = reshape(coords(1,:), size(x)) + target(1);
y = reshape(coords(2,:), size(y)) + target(2);
z = reshape(coords(3,:), size(z)) + target(3);
surf(x, y, z, 'FaceAlpha', 0.3, 'EdgeColor', 'b', 'FaceColor', 'y');
hold on;
scatter3(target(1), target(2), target(3), 100, 'r*', 'LineWidth', 2);
end脚本2：多基地RCS闪烁与NLOS路径利用
内容说明：实现多基地雷达目标截面积的空间分布建模与闪烁特性仿真，包含双基地RCS计算、前向散射增强效应与NLOS多径利用。通过反射面镜像方法构建虚拟发射源，实现非视距路径的参数估计与目标定位。
使用方式：运行main_Multistatic_RCS_NLOS，配置targetGeometry（目标模型）、bistaticAngles（观测角度范围）与reflectors（环境反射面），分析RCS闪烁统计与NLOS定位精度。
matlab
复制
%% 脚本2：多基地RCS闪烁与NLOS路径利用
% 功能：双基地RCS建模，前向散射增强，NLOS多径利用与镜像定位
% 使用：配置targetGeometry(目标结构)、bistaticAngles(观测角)、reflectors(反射面)
% 输出：RCS方向图、闪烁统计、NLOS路径增益、镜像定位结果
function main_Multistatic_RCS_NLOS()
clear; close all; clc;
%% 系统配置
fc = 10e9; % 10GHz
wavelength = 3e8 / fc;
% 目标配置（简化舰船模型：三个散射中心）
targetCenters = [0, 0, 0; % 船体中心
50, 0, 10; % 舰桥
-30, 20, 5; % 桅杆
20, -15, 8]; % 雷达罩
targetRCS = [100, 50, 30, 40]; % 各部件RCS (m^2)
% 多基地配置
txPos = [0, -5000, 100]; % 发射站 5km外
rxPosGrid = generateRxGrid(0, 0, 10000, 50); % 半圆接收阵
% 反射面配置（NLOS场景）
reflectors = [1000, 2000, 0, 3000, 2000, 0]; % 墙面起点(x1,y1,z1),终点(x2,y2,z2)
%% 计算双基地RCS
bistaticRCS = calculateBistaticRCS(targetCenters, targetRCS, txPos, rxPosGrid, wavelength);
%% NLOS路径分析
[nlosPaths, virtualTx] = calculateNLOSPaths(txPos, reflectors, [0, 0, 0]); % 假设目标在原点
%% 可视化
figure('Name', '多基地RCS与NLOS利用', 'Position', [50 50 1400 900]);
% 双基地RCS方向图（极坐标）
subplot(3, 3, 1);
polarplot(linspace(0, 2*pi, length(bistaticRCS)), 10*log10(bistaticRCS), 'b-', 'LineWidth', 2);
title('双基地RCS方向图 (dBsm)');
% RCS vs 双基地角
subplot(3, 3, 2);
bistaticAngles = calculateBistaticAngles(txPos, rxPosGrid, [0,0,0]);
plot(bistaticAngles * 180/pi, 10*log10(bistaticRCS), 'r.-', 'LineWidth', 1.5);
xline(180, 'g--', '前向散射');
title('RCS vs 双基地角');
xlabel('双基地角 (度)'); ylabel('RCS (dBsm)');
grid on;
% 闪烁统计（模拟姿态变化）
subplot(3, 3, 3);
flashSamples = simulateRCSFluctuation(targetCenters, targetRCS, txPos, rxPosGrid(10,:), wavelength, 100);
histogram(10*log10(flashSamples), 20);
title('RCS闪烁分布');
xlabel('RCS (dBsm)'); ylabel('频次');
% 前向散射增强验证
subplot(3, 3, 4);
fwdScatter = calculateForwardScatter(txPos, rxPosGrid, targetCenters, wavelength);
plot(fwdScatter.angles * 180/pi, 10*log10(fwdScatter.rcs), 'm-', 'LineWidth', 2);
title('前向散射增强 (Po常数)');
xlabel('双基地角 (度)'); ylabel('增强RCS (dBsm)');
grid on;
% 几何配置与反射面
subplot(3, 3, 5);
scatter(txPos(1)/1000, txPos(2)/1000, 300, 'r^', 'filled'); hold on;
scatter(rxPosGrid(:,1)/1000, rxPosGrid(:,2)/1000, 100, 'bo', 'filled');
scatter(targetCenters(:,1)/1000, targetCenters(:,2)/1000, 200, 'gs', 'filled');
% 绘制反射面
for i = 1:size(reflectors,1)
plot([reflectors(i,1), reflectors(i,4)]/1000, ...
[reflectors(i,2), reflectors(i,5)]/1000, 'k-', 'LineWidth', 3);
end
title('多基地几何与反射面');
xlabel('X (km)'); ylabel('Y (km)');
legend('发射站', '接收站', '目标'); axis equal; grid on;
% NLOS路径可视化
subplot(3, 3, 6);
scatter(txPos(1)/1000, txPos(2)/1000, 300, 'r^', 'filled'); hold on;
scatter(virtualTx(:,1)/1000, virtualTx(:,2)/1000, 300, 'm^', 'filled');
plot([txPos(1), virtualTx(1,1)]/1000, [txPos(2), virtualTx(1,2)]/1000, 'm--');
title('虚拟发射源（镜像）');
xlabel('X (km)'); ylabel('Y (km)');
legend('真实发射', '虚拟发射'); axis equal; grid on;
% NLOS路径增益
subplot(3, 3, 7);
bar(nlosPaths.gains);
title('NLOS路径增益');
xlabel('路径编号'); ylabel('相对增益 (dB)');
% LOS vs NLOS定位精度
subplot(3, 3, 8);
[losErr, nlosErr] = compareLOSNLOS(txPos, rxPosGrid(1,:), reflectors, targetCenters, wavelength);
bar([losErr, nlosErr]);
title('定位误差对比 (m)');
set(gca, 'XTickLabel', {'仅有LOS', '利用NLOS'});
% 多径分集增益
subplot(3, 3, 9);
diversityGain = calculateDiversityGain(nlosPaths.gains);
plot(diversityGain, 'b-o', 'LineWidth', 2);
title('多径分集增益 vs 路径数');
xlabel('路径数'); ylabel('累积增益 (dB)');
grid on;
%% 统计输出
fprintf('=== 多基地RCS与NLOS分析 ===\n');
fprintf('平均双基地RCS: %.2f dBsm\n', 10*log10(mean(bistaticRCS)));
fprintf('前向散射峰值: %.2f dBsm @ 180°\n', 10*log10(max(fwdScatter.rcs)));
fprintf('NLOS路径数: %d\n', length(nlosPaths.gains));
fprintf('NLOS定位改善: %.1f%% (相对纯LOS)\n', (1-nlosErr/losErr)*100);
end
%% 计算双基地RCS（简化物理光学）
function rcs = calculateBistaticRCS(centers, rcsValues, txPos, rxPos, lambda)
numRx = size(rxPos, 1);
rcs = zeros(numRx, 1);
for i = 1:numRx
% 计算各散射中心贡献
totalField = 0;
for j = 1:size(centers, 1)
r1 = norm(centers(j,:) - txPos);
r2 = norm(centers(j,:) - rxPos(i,:));
% 相位
phase = 2*pi*(r1 + r2) / lambda;
% 幅度（简化：各向同性）
amp = sqrt(rcsValues(j)) / (r1 * r2);
totalField = totalField + amp * exp(1j * phase);
end
rcs(i) = abs(totalField)^2 * (4*pi);
end
end
%% 生成接收站网格（半圆）
function grid = generateRxGrid(cx, cy, radius, numPoints)
angles = linspace(0, pi, numPoints);
grid = [cx + radius * cos(angles'), cy + radius * sin(angles'), zeros(numPoints, 1)];
end
%% 计算双基地角
function angles = calculateBistaticAngles(txPos, rxPos, target)
numRx = size(rxPos, 1);
angles = zeros(numRx, 1);
v1 = target - txPos;
for i = 1:numRx
v2 = target - rxPos(i,:);
angles(i) = acos(dot(v1, v2) / (norm(v1) * norm(v2)));
end
end
%% 模拟RCS闪烁（姿态扰动）
function samples = simulateRCSFluctuation(centers, rcsVals, txPos, rxPos, lambda, numSamples)
samples = zeros(numSamples, 1);
for k = 1:numSamples
% 随机旋转目标（小角度扰动）
theta = rand() * 0.1; % 随机偏航
R = [cos(theta), -sin(theta), 0; sin(theta), cos(theta), 0; 0, 0, 1];
rotatedCenters = (R * centers')';
samples(k) = calculateBistaticRCS(rotatedCenters, rcsVals, txPos, rxPos, lambda);
end
end
%% 计算前向散射（Po常数近似）
function result = calculateForwardScatter(txPos, rxPos, targetShape, lambda)
% 使用物理光学前向散射近似
numRx = size(rxPos, 1);
angles = calculateBistaticAngles(txPos, rxPos, mean(targetShape));
% 前向散射RCS（双基地角接近180°）
% sigma_fs = 4*pi*A^2/lambda^2, A为目标投影面积
A = 1000; % 投影面积 m^2（简化）
sigma_fs = 4 * pi * A^2 / lambda^2;
% 随角度衰减（简化模型）
rcs = sigma_fs * exp(-((angles - pi).^2) / 0.01); % 高斯衰减
result.angles = angles;
result.rcs = rcs;
end
%% 计算NLOS路径（镜像法）
function [paths, virtualTx] = calculateNLOSPaths(txPos, reflectors, targetPos)
numReflectors = size(reflectors, 1);
virtualTx = zeros(numReflectors, 3);
paths.gains = zeros(numReflectors, 1);
paths.delays = zeros(numReflectors, 1);
for i = 1:numReflectors
% 计算反射面法向
p1 = reflectors(i, 1:3);
p2 = reflectors(i, 4:6);
normal = cross(p2 - p1, [0, 0, 1]); % 假设反射面垂直
normal = normal / norm(normal);
% 镜像点
d = dot(p1 - txPos, normal);
virtualTx(i, :) = txPos + 2 * d * normal;
% 路径增益（反射损耗）
pathLen1 = norm(virtualTx(i,:) - targetPos);
pathLen2 = norm(txPos - targetPos);
paths.gains(i) = 20*log10(pathLen2/pathLen1) - 3; % 反射损耗3dB
paths.delays(i) = (pathLen1 - pathLen2) / 3e8;
end
end
%% 比较LOS与NLOS定位
function [losErr, nlosErr] = compareLOSNLOS(txPos, rxPos, reflectors, target, lambda)
% 纯LOS定位（单站）
trueRange = norm(target(1,:) - txPos) + norm(target(1,:) - rxPos);
% 添加误差
losErr = 10; % 假设10m误差（单站）
% NLOS定位（多径融合）
[paths, virtualTx] = calculateNLOSPaths(txPos, reflectors, target(1,:));
% 多径联合定位降低误差
nlosErr = losErr / sqrt(1 + length(paths.gains)); % 分集增益
end
%% 计算多径分集增益
function gain = calculateDiversityGain(pathGains)
% MRC分集增益累积
linearGains = 10.^(pathGains/10);
cumsumGain = cumsum(linearGains);
gain = 10*log10(cumsumGain / linearGains(1));
end脚本3：直达波与目标回波时间同步
内容说明：实现多基地雷达系统的直达波提取与目标回波时间同步算法，包含直达波检测、时间延迟估计、互相关处理与相位同步误差补偿。支持双向时间传递与基于直达波的时钟校准。
使用方式：运行main_Direct_Path_Sync，配置snrDirect（直达波信噪比）、timeOffset（初始时间偏移）与clockDrift（时钟漂移率），系统执行同步算法并评估残余误差。
matlab
复制
%% 脚本3：直达波与目标回波时间同步
% 功能：直达波提取，双向时间传递，时钟漂移补偿，相位同步
% 使用：配置snrDirect(dB)、timeOffset(ns)、clockDrift(ppm)、signalBandwidth(MHz)
% 输出：同步误差收敛曲线、相位噪声谱、TDOA/FDOA估计精度
function main_Direct_Path_Sync()
clear; close all; clc;
%% 系统配置
fs = 100e6; % 采样率 100MHz
T = 1e-3; % 观测时间 1ms
N = fs * T; % 采样点数
bw = 20e6; % 信号带宽 20MHz
% 时钟参数
trueOffset = 100e-9; % 真实时间偏移 100ns
trueDrift = 1e-6; % 真实时钟漂移 1ppm
phaseNoiseLevel = -90; % dBc/Hz @ 1kHz
% 信道参数
directPathDelay = 100e-6; % 直达波延迟 100us（33km）
targetDelay = 150e-6; % 目标回波延迟 150us
snrDirect = 20; % 直达波SNR
snrTarget = 10; % 目标回波SNR
%% 生成信号
% 发射信号（LFM）
t = (0:N-1)'/fs;
txSignal = chirp(t, 0, T, bw, 'linear');
% 接收信号（含时钟偏移与漂移）
rxClock = (1+trueDrift)*t + trueOffset;
directWave = interp1(t, txSignal, rxClock - directPathDelay, 'linear', 0);
targetEcho = 0.5 * interp1(t, txSignal, rxClock - targetDelay, 'linear', 0);
% 添加相位噪声
phaseNoise = generatePhaseNoise(N, fs, phaseNoiseLevel);
directWave = directWave .* exp(1j * phaseNoise);
targetEcho = targetEcho .* exp(1j * phaseNoise);
% 添加噪声
directWave = awgn(directWave, snrDirect, 'measured');
targetEcho = awgn(targetEcho, snrTarget, 'measured');
rxSignal = directWave + targetEcho;
%% 时间同步处理
% 1. 直达波提取（匹配滤波）
[directFiltered, directTemplate] = extractDirectPath(rxSignal, txSignal, fs, directPathDelay, bw);
% 2. 时间延迟估计（Early-Late门）
[estDelay, trackingError] = earlyLateGate(directFiltered, txSignal, fs, 100);
% 3. 双向时间传递（模拟）
[syncError, driftEst] = twoWayTimeTransfer(txSignal, rxSignal, fs, trueOffset, trueDrift);
% 4. 目标回波对齐
alignedTarget = alignTargetEcho(rxSignal, directFiltered, targetDelay, estDelay);
%% 可视化
figure('Name', '直达波与目标回波时间同步', 'Position', [50 50 1400 900]);
% 时域波形
subplot(3, 3, 1);
plot(t*1e6, real(rxSignal), 'b', 'LineWidth', 0.5); hold on;
plot(t*1e6, real(directWave), 'r', 'LineWidth', 1);
xlim([90, 160]);
title('接收信号时域（含直达波与目标）');
xlabel('时间 (\mus)'); ylabel('幅度');
legend('总接收', '直达波'); grid on;
% 匹配滤波输出（直达波检测）
subplot(3, 3, 2);
mfOutput = abs(xcorr(rxSignal, txSignal, 'coeff'));
lag = (-(N-1):(N-1))'/fs * 1e6;
plot(lag, mfOutput, 'b-', 'LineWidth', 1.5);
xline(directPathDelay*1e6, 'r--', '真实直达波');
xline(targetDelay*1e6, 'g--', '真实目标');
title('匹配滤波输出（直达波检测）');
xlabel('时延 (\mus)'); ylabel('相关幅度');
grid on; xlim([80, 200]);
% Early-Late跟踪误差
subplot(3, 3, 3);
plot(1:length(trackingError), trackingError*1e9, 'b-', 'LineWidth', 1.5);
title('Early-Late跟踪误差');
xlabel('迭代'); ylabel('误差 (ns)');
grid on;
% 双向时间传递收敛
subplot(3, 3, 4);
plot(syncError*1e9, 'b-', 'LineWidth', 2);
title('双向时间传递误差收敛');
xlabel('迭代'); ylabel('时间误差 (ns)');
grid on;
yline(0, 'r--', '零误差');
% 时钟漂移估计
subplot(3, 3, 5);
bar([trueDrift*1e6, driftEst*1e6]);
title('时钟漂移估计');
set(gca, 'XTickLabel', {'真实', '估计'});
ylabel('漂移 (ppm)');
% 相位噪声谱
subplot(3, 3, 6);
[pxx, f] = pwelch(phaseNoise, hamming(256), 128, 1024, fs, 'centered');
semilogy(f/1e3, pxx, 'b-', 'LineWidth', 1.5);
title('相位噪声功率谱');
xlabel('频率偏移 (kHz)'); ylabel('PSD (rad^2/Hz)');
grid on;
% 目标回波对齐验证
subplot(3, 3, 7);
plot(t*1e6, real(alignedTarget), 'b-', 'LineWidth', 1); hold on;
plot(t*1e6, real(targetEcho), 'r--', 'LineWidth', 1);
xlim([140, 160]);
title('目标回波对齐验证');
xlabel('时间 (\mus)'); ylabel('幅度');
legend('对齐后', '理想位置'); grid on;
% TDOA估计精度
subplot(3, 3, 8);
tdoaError = abs((estDelay - directPathDelay) - (targetDelay - directPathDelay));
bar(tdoaError*1e9);
title('TDOA估计误差');
ylabel('误差 (ns)');
% 同步后相参积累增益
subplot(3, 3, 9);
phaseErrors = linspace(0, pi, 100);
gainLoss = 20*log10(abs(sinc(phaseErrors/pi))); % 简化模型
plot(phaseErrors*180/pi, gainLoss, 'b-', 'LineWidth', 2);
title('相位误差 vs 积累损失');
xlabel('相位误差 (度)'); ylabel('增益损失 (dB)');
grid on;
xline(10, 'r--', '目标容限');
%% 统计输出
fprintf('=== 时间同步性能 ===\n');
fprintf('初始时间偏移: %.2f ns\n', trueOffset*1e9);
fprintf('估计时间偏移: %.2f ns (误差: %.2f ns)\n', ...
(estDelay-directPathDelay)*1e9, abs(estDelay-directPathDelay-trueOffset)*1e9);
fprintf('时钟漂移估计误差: %.3f ppm\n', abs(driftEst - trueDrift)*1e6);
fprintf('同步后TDOA精度: %.2f ns\n', tdoaError*1e9);
fprintf('等效距离精度: %.2f m\n', tdoaError*3e8);
end
%% 生成相位噪声（Leeson模型简化）
function phi = generatePhaseNoise(N, fs, noiseLevel)
% 简化：白噪声+闪烁噪声
whiteNoise = randn(N, 1) * 10^(noiseLevel/20);
% 积分得到相位（布朗运动）
phi = cumsum(whiteNoise) * sqrt(1/fs);
end
%% 直达波提取（带通滤波+匹配滤波）
function [filtered, template] = extractDirectPath(rx, tx, fs, expectedDelay, bw)
% 粗略检测
corr = xcorr(rx, tx);
[~, peakIdx] = max(abs(corr));
lag = -(length(rx)-1):(length(rx)-1);
estDelay = lag(peakIdx) / fs;
% 提取窗口
windowStart = round(estDelay * fs) + 1;
windowLen = length(tx);
if windowStart > 0 && windowStart + windowLen - 1 <= length(rx)
filtered = rx(windowStart:windowStart+windowLen-1);
else
filtered = zeros(size(tx));
end
template = tx;
end
%% Early-Late门跟踪
function [estDelay, errorHistory] = earlyLateGate(signal, template, fs, numIter)
delta = 1/fs; % 早/晚门间隔
delay = 0;
errorHistory = zeros(numIter, 1);
for i = 1:numIter
% Early
sEarly = interp1(1:length(signal), signal, (1:length(template))' - delta, 'linear', 0);
eCorr = abs(sum(sEarly .* conj(template)));
% Late
sLate = interp1(1:length(signal), signal, (1:length(template))' + delta, 'linear', 0);
lCorr = abs(sum(sLate .* conj(template)));
% 误差
error = (eCorr - lCorr) / (eCorr + lCorr + eps);
errorHistory(i) = error;
% 更新（PI控制器简化）
delay = delay + 0.1 * error * delta;
if abs(error) < 0.01
break;
end
end
estDelay = delay;
end
%% 双向时间传递（模拟）
function [errorHistory, driftEst] = twoWayTimeTransfer(tx, rx, fs, trueOffset, trueDrift)
% 迭代估计偏移与漂移
estOffset = 0;
estDrift = 0;
numIter = 50;
errorHistory = zeros(numIter, 1);
for i = 1:numIter
% 模拟测量（含噪声）
noise = randn() * 1e-9;
meas = trueOffset + rand() * trueDrift * 1e-3 + noise;
% 卡尔曼滤波简化
alpha = 0.1;
estOffset = (1-alpha) * estOffset + alpha * meas;
estDrift = estDrift + 0.01 * (meas - estOffset);
errorHistory(i) = estOffset - trueOffset;
end
driftEst = estDrift;
end
%% 对齐目标回波
function aligned = alignTargetEcho(rx, directRef, targetDelayTrue, estDirectDelay)
% 根据估计的直达波位置对齐目标
offset = round((targetDelayTrue - estDirectDelay) * 100e6); % 简化
aligned = circshift(rx, offset);
end脚本4：相参积累相位同步误差补偿
内容说明：实现分布式相参积累的相位同步与误差补偿算法，包含相位误差估计、自适应相位对齐与相参积累增益分析。支持多接收站信号的联合处理与合成波束形成。
使用方式：运行main_Coherent_Integration_PhaseSync，配置numNodes（节点数）、phaseErrorStd（相位误差标准差）与snr（信噪比），算法执行相位同步并展示积累增益与误码率改善。
matlab
复制
%% 脚本4：相参积累相位同步误差补偿
% 功能：分布式相参积累，相位误差估计与补偿，合成波束形成
% 使用：配置numNodes(节点数)、phaseErrorStd(度)、snr(dB)、coherentTime(ms)
% 输出：相位收敛曲线、积累增益、合成方向图、检测性能提升
function main_Coherent_Integration_PhaseSync()
clear; close all; clc;
%% 网络配置
numNodes = 8; % 接收节点数
nodePositions = 1000 * [cos(2*pi*(0:numNodes-1)/numNodes);
sin(2*pi*(0:numNodes-1)/numNodes);
zeros(1, numNodes)]'; % 圆阵，半径1km
targetPos = [5000, 3000, 1000]; % 目标位置 5km外
fc = 10e9; % 10GHz
wavelength = 3e8 / fc;
snr = -10; % 单通道SNR (dB)
% 相位误差
truePhaseErrors = deg2rad(30 * randn(numNodes, 1)); % 初始误差±30度
%% 生成接收信号
[rxSignals, steeringVector] = generateArraySignals(numNodes, nodePositions, targetPos, wavelength, snr, truePhaseErrors);
%% 相位同步算法
% 1. 基于参考源的相位估计（假设有校准信号）
estPhaseErrors = estimatePhaseErrors(rxSignals, steeringVector);
% 2. 相位补偿与相参积累
compensatedSignals = compensatePhaseErrors(rxSignals, estPhaseErrors);
% 3. 合成波束形成
beamformedOutput = coherentBeamforming(compensatedSignals, steeringVector);
%% 性能分析
% 相参积累增益
cohGain = calculateCoherentGain(compensatedSignals, numNodes);
% 检测性能（ROC曲线）
[pdNonCoh, pdCoh, pdIdeal] = calculateDetectionPerformance(numNodes, snr, truePhaseErrors, estPhaseErrors);
%% 可视化
figure('Name', '相参积累相位同步', 'Position', [50 50 1400 900]);
% 网络几何
subplot(3, 3, 1);
scatter(nodePositions(:,1)/1000, nodePositions(:,2)/1000, 300, 'bo', 'filled'); hold on;
scatter(targetPos(1)/1000, targetPos(2)/1000, 200, 'r*', 'LineWidth', 3);
plot([nodePositions(:,1), targetPos(1)*ones(numNodes,1)]'/1000, ...
[nodePositions(:,2), targetPos(2)*ones(numNodes,1)]'/1000, 'g:');
title('分布式阵列几何');
xlabel('X (km)'); ylabel('Y (km)');
legend('接收节点', '目标'); axis equal; grid on;
% 真实 vs 估计相位误差
subplot(3, 3, 2);
plot(rad2deg(truePhaseErrors), 'bo-', 'LineWidth', 2); hold on;
plot(rad2deg(estPhaseErrors), 'rs--', 'LineWidth', 2);
title('相位误差估计');
xlabel('节点'); ylabel('相位误差 (度)');
legend('真实', '估计'); grid on;
% 残余相位误差
subplot(3, 3, 3);
residual = rad2deg(truePhaseErrors - estPhaseErrors);
bar(residual);
title('残余相位误差');
xlabel('节点'); ylabel('误差 (度)');
grid on;
% 相参积累增益
subplot(3, 3, 4);
bar([numNodes, numNodes^2 * 10^(-var(truePhaseErrors)/2), cohGain]);
title('积累增益对比');
set(gca, 'XTickLabel', {'非相参(N)', '理想相参(N^2)', '实际相参'});
ylabel('增益');
% 合成方向图
subplot(3, 3, 5);
angles = linspace(-pi, pi, 360);
pattern = zeros(size(angles));
for i = 1:length(angles)
testDir = [cos(angles(i)), sin(angles(i)), 0];
sv = exp(-1j*2*pi/wavelength * nodePositions * testDir');
pattern(i) = abs(sum(exp(-1j*estPhaseErrors) .* sv))^2;
end
polarplot(angles, 10*log10(pattern/max(pattern)), 'b-', 'LineWidth', 2);
title('合成波束方向图 (dB)');
% 相位收敛过程（迭代）
subplot(3, 3, 6);
phaseHistory = simulatePhaseConvergence(numNodes, truePhaseErrors, 50);
plot(1:50, rad2deg(phaseHistory'), 'LineWidth', 1.5);
title('相位同步收敛过程');
xlabel('迭代'); ylabel('相位估计 (度)');
grid on;
% ROC曲线
subplot(3, 3, 7);
snrRange = -20:2:0;
for i = 1:length(snrRange)
[pdc, pdn] = calculateROC(numNodes, snrRange(i), truePhaseErrors, estPhaseErrors);
pdCoh(i) = pdc;
pdNonCoh(i) = pdn;
end
plot(snrRange, pdNonCoh, 'b-o', 'LineWidth', 2); hold on;
plot(snrRange, pdCoh, 'r-s', 'LineWidth', 2);
title('检测性能对比 (ROC)');
xlabel('单通道SNR (dB)'); ylabel('检测概率');
legend('非相参', '相参'); grid on;
% 相位噪声影响
subplot(3, 3, 8);
phaseNoiseRange = 0:5:60;
gainLoss = 20*log10(exp(-(deg2rad(phaseNoiseRange)).^2/2));
plot(phaseNoiseRange, gainLoss, 'b-', 'LineWidth', 2);
title('相位误差 vs 积累损失');
xlabel('相位误差标准差 (度)'); ylabel('增益损失 (dB)');
grid on;
% 信号矢量图（星座图）
subplot(3, 3, 9);
scatter(real(rxSignals(:,100)), imag(rxSignals(:,100)), 'b', 'filled'); hold on;
scatter(real(compensatedSignals(:,100)), imag(compensatedSignals(:,100)), 'r', 'filled');
plot([0, sum(real(compensatedSignals(:,100)))], [0, sum(imag(compensatedSignals(:,100)))], 'g-', 'LineWidth', 3);
title('信号矢量合成');
xlabel('实部'); ylabel('虚部');
legend('原始', '补偿后', '合成');
axis equal; grid on;
%% 统计输出
fprintf('=== 相参积累相位同步 ===\n');
fprintf('节点数: %d\n', numNodes);
fprintf('初始相位误差(RMS): %.2f°\n', rms(rad2deg(truePhaseErrors)));
fprintf('残余相位误差(RMS): %.2f°\n', rms(abs(residual)));
fprintf('相参积累增益: %.2f (理想: %d)\n', cohGain, numNodes^2);
fprintf('增益损失: %.2f dB\n', 20*log10(numNodes^2/cohGain));
end
%% 生成阵列信号
function [signals, sv] = generateArraySignals(numNodes, positions, target, lambda, snr, phaseErrors)
% 导向矢量
direction = (target - mean(positions)) / norm(target - mean(positions));
sv = exp(-1j * 2*pi/lambda * positions * direction');
% 生成信号
N = 1000; % 采样点
signalPower = 10^(snr/10);
s = sqrt(signalPower) * (randn(N, 1) + 1j*randn(N, 1));
% 各通道接收（含相位误差与噪声）
signals = zeros(numNodes, N);
for i = 1:numNodes
noise = (randn(N, 1) + 1j*randn(N, 1)) / sqrt(2);
signals(i, :) = (sv(i) * exp(1j*phaseErrors(i)) * s.' + noise).';
end
end
%% 估计相位误差（基于参考信号或互相关）
function est = estimatePhaseErrors(signals, sv)
numNodes = size(signals, 1);
est = zeros(numNodes, 1);
% 使用第一节点作为参考
ref = signals(1, :);
for i = 1:numNodes
% 互相关估计相位差
corr = sum(signals(i, :) .* conj(ref));
est(i) = angle(corr) - angle(sv(i) * conj(sv(1)));
end
est = est - est(1); % 相对相位
end
%% 补偿相位误差
function comp = compensatePhaseErrors(signals, est)
comp = signals .* exp(-1j * est);
end
%% 合成波束形成
function out = coherentBeamforming(signals, sv)
% 最优合并（MVDR简化）
weights = sv / norm(sv);
out = sum(weights .* signals, 1);
end
%% 计算相参增益
function gain = calculateCoherentGain(signals, N)
% 实际增益
combined = sum(signals, 1);
signalPower = mean(abs(combined).^2);
noisePower = N; % 简化
gain = signalPower / noisePower;
end
%% 模拟相位收敛
function history = simulatePhaseConvergence(numNodes, trueErrors, numIter)
history = zeros(numNodes, numIter);
est = zeros(numNodes, 1);
for i = 1:numIter
% 带噪声的相位测量
noise = deg2rad(5) * randn(numNodes, 1);
meas = trueErrors + noise;
% 滤波估计
alpha = 0.2;
est = (1-alpha) * est + alpha * meas;
history(:, i) = est;
end
end
%% 计算ROC性能
function [pdCoh, pdNonCoh] = calculateROC(N, snr, trueErr, estErr)
% 简化计算（高斯近似）
snrLin = 10^(snr/10);
% 非相参：功率积累
snrNonCoh = N * snrLin;
pdNonCoh = 0.5 * erfc(-sqrt(snrNonCoh/2));
% 相参：电压积累，考虑残余相位误差
residual = trueErr - estErr;
coherency = abs(sum(exp(-1j*residual))) / N;
snrCoh = (N * coherency)^2 * snrLin;
pdCoh = 0.5 * erfc(-sqrt(snrCoh/2));
end脚本5：非相参积累GLRT与CFAR检测
内容说明：实现多基地雷达非相参积累的广义似然比检验（GLRT）与恒虚警率（CFAR）检测，包含平方律检波、单元平均CFAR与有序统计CFAR，分析多通道联合检测性能。
使用方式：运行main_Noncoherent_GLRT_CFAR，配置numChannels（通道数）、snr（信噪比）、numReferenceCells（参考单元数）与guardCells（保护单元数），系统生成检测统计量并计算虚警率与检测概率。
matlab
复制
%% 脚本5：非相参积累GLRT与CFAR检测
% 功能：多基地非相参积累，平方律检波，GLRT实现，CA-CFAR与OS-CFAR
% 使用：配置numChannels(通道数)、snr(dB)、numReferenceCells、guardCells、Pfa
% 输出：检测统计量分布、CFAR门限、ROC曲线、虚警控制
function main_Noncoherent_GLRT_CFAR()
clear; close all; clc;
%% 参数配置
numChannels = 4; % 接收通道数
numSamples = 1000; % 距离单元数
numTargets = 5; % 目标数
snr = 15; % 单通道SNR (dB)
PfaTarget = 1e-6; % 目标虚警率
numRefCells = 32; % 参考单元数（每侧16）
numGuardCells = 4; % 保护单元数（每侧2）
%% 生成数据
% 噪声背景
noise = randn(numChannels, numSamples) + 1j*randn(numChannels, numSamples);
% 添加目标（随机位置）
signal = noise;
targetIdx = randperm(numSamples, numTargets);
for i = 1:numTargets
amp = 10^(snr/20);
phase = 2*pi*rand(numChannels, 1); % 非相参：随机相位
signal(:, targetIdx(i)) = signal(:, targetIdx(i)) + amp * exp(1j*phase);
end
%% 平方律检波（每通道）
squaredData = abs(signal).^2;
%% 非相参积累（各通道求和）
integratedData = sum(squaredData, 1); % 非相参积累
%% CFAR处理（单元平均）
[detectionCA, thresholdCA] = caCFAR(integratedData, numRefCells, numGuardCells, PfaTarget);
%% CFAR处理（有序统计OS-CFAR）
[detectionOS, thresholdOS] = osCFAR(integratedData, numRefCells, numGuardCells, PfaTarget, 0.75);
%% GLRT检测（假设检验框架）
glrtStatistic = calculateGLRT(squaredData, snr);
%% 性能评估
[pdCA, pfaCA] = evaluatePerformance(detectionCA, targetIdx, numSamples);
[pdOS, pfaOS] = evaluatePerformance(detectionOS, targetIdx, numSamples);
%% 可视化
figure('Name', '非相参积累GLRT与CFAR', 'Position', [50 50 1400 900]);
% 积累后数据
subplot(3, 3, 1);
plot(1:numSamples, 10*log10(integratedData), 'b-', 'LineWidth', 0.8); hold on;
plot(targetIdx, 10*log10(integratedData(targetIdx)), 'ro', 'MarkerSize', 8, 'LineWidth', 2);
plot(1:numSamples, 10*log10(thresholdCA), 'r--', 'LineWidth', 1.5);
title('非相参积累与CA-CFAR门限');
xlabel('距离单元'); ylabel('功率 (dB)');
legend('积累输出', '目标', 'CFAR门限'); grid on;
% 各通道数据（瀑布图）
subplot(3, 3, 2);
imagesc(1:numSamples, 1:numChannels, 10*log10(squaredData));
title('各通道平方律检波输出');
xlabel('距离单元'); ylabel('通道');
colorbar; colormap hot;
% CA-CFAR门限
subplot(3, 3, 3);
plot(1:numSamples, 10*log10(thresholdCA), 'b-', 'LineWidth', 1.5); hold on;
plot(1:numSamples, 10*log10(thresholdOS), 'r--', 'LineWidth', 1.5);
title('CFAR门限对比');
xlabel('距离单元'); ylabel('门限 (dB)');
legend('CA-CFAR', 'OS-CFAR'); grid on;
% 检测统计量分布（H0 vs H1）
subplot(3, 3, 4);
h0Data = integratedData(setdiff(1:numSamples, targetIdx));
h1Data = integratedData(targetIdx);
histogram(h0Data, 30, 'Normalization', 'pdf', 'FaceColor', 'b', 'FaceAlpha', 0.5); hold on;
histogram(h1Data, 10, 'Normalization', 'pdf', 'FaceColor', 'r', 'FaceAlpha', 0.5);
xline(mean(thresholdCA), 'g--', '门限');
title('检测统计量分布');
xlabel('积累功率'); ylabel('PDF');
legend('H0 (噪声)', 'H1 (信号+噪声)');
% GLRT统计量
subplot(3, 3, 5);
plot(1:numSamples, glrtStatistic, 'b-', 'LineWidth', 0.8); hold on;
plot(targetIdx, glrtStatistic(targetIdx), 'ro', 'MarkerSize', 8);
title('GLRT检验统计量');
xlabel('距离单元'); ylabel('似然比');
grid on;
% ROC曲线（不同SNR）
subplot(3, 3, 6);
snrRange = 0:2:20;
pdCA = zeros(size(snrRange));
pdOS = zeros(size(snrRange));
for i = 1:length(snrRange)
[pdCA(i), ~] = simulateROC(numChannels, snrRange(i), PfaTarget, 'CA');
[pdOS(i), ~] = simulateROC(numChannels, snrRange(i), PfaTarget, 'OS');
end
plot(snrRange, pdCA, 'b-o', 'LineWidth', 2); hold on;
plot(snrRange, pdOS, 'r-s', 'LineWidth', 2);
title('ROC曲线（不同SNR）');
xlabel('SNR (dB)'); ylabel('检测概率');
legend('CA-CFAR', 'OS-CFAR'); grid on;
% 虚警率控制验证
subplot(3, 3, 7);
pfaRange = logspace(-6, -3, 20);
actualPfaCA = zeros(size(pfaRange));
for i = 1:length(pfaRange)
[~, actualPfaCA(i)] = simulateROC(numChannels, snr, pfaRange(i), 'CA');
end
loglog(pfaRange, actualPfaCA, 'b-o', 'LineWidth', 2); hold on;
loglog(pfaRange, pfaRange, 'k--');
title('虚警率控制验证');
xlabel('目标Pfa'); ylabel('实际Pfa');
legend('实际', '理想'); grid on; axis equal;
% 多通道分集增益
subplot(3, 3, 8);
channelRange = 1:8;
reqSNR = zeros(size(channelRange));
for i = 1:length(channelRange)
reqSNR(i) = findRequiredSNR(channelRange(i), 0.9, 1e-6, 'noncoherent');
end
plot(channelRange, reqSNR, 'b-o', 'LineWidth', 2);
title('所需SNR vs 通道数 (Pd=0.9, Pfa=1e-6)');
xlabel('通道数'); ylabel('所需SNR (dB)');
grid on;
% 检测总结
subplot(3, 3, 9);
text(0.1, 0.8, '检测性能总结:', 'FontSize', 12, 'FontWeight', 'bold');
text(0.1, 0.6, sprintf('CA-CFAR: Pd=%.2f, Pfa=%.2e', pdCA, pfaCA), 'FontSize', 10);
text(0.1, 0.5, sprintf('OS-CFAR: Pd=%.2f, Pfa=%.2e', pdOS, pfaOS), 'FontSize', 10);
text(0.1, 0.3, sprintf('非相参积累增益: %.1fdB (理论%.1fdB)', ...
10*log10(numChannels), 10*log10(numChannels)), 'FontSize', 10);
axis off;
end
%% 单元平均CFAR
function [detection, threshold] = caCFAR(data, numRef, numGuard, Pfa)
N = length(data);
detection = zeros(size(data));
threshold = zeros(size(data));
alpha = numRef * (Pfa^(-1/numRef) - 1); % 门限因子
for i = 1:N
% 参考窗口
leftStart = max(1, i - numGuard - numRef);
leftEnd = max(1, i - numGuard);
rightStart = min(N, i + numGuard + 1);
rightEnd = min(N, i + numGuard + numRef);
refCells = [data(leftStart:leftEnd), data(rightStart:rightEnd)];
% 门限
meanLevel = mean(refCells);
threshold(i) = alpha * meanLevel;
% 检测
if data(i) > threshold(i)
detection(i) = 1;
end
end
end
%% 有序统计CFAR
function [detection, threshold] = osCFAR(data, numRef, numGuard, Pfa, rankPercent)
N = length(data);
detection = zeros(size(data));
threshold = zeros(size(data));
k = round(rankPercent * numRef); % 排序后第k个值
% OS-CFAR门限因子（近似）
alpha = k * (Pfa^(-1/k) - 1) * 1.2; % 修正因子
for i = 1:N
leftStart = max(1, i - numGuard - numRef);
leftEnd = max(1, i - numGuard);
rightStart = min(N, i + numGuard + 1);
rightEnd = min(N, i + numGuard + numRef);
refCells = [data(leftStart:leftEnd), data(rightStart:rightEnd)];
sortedCells = sort(refCells);
if length(sortedCells) >= k
osLevel = sortedCells(k);
else
osLevel = sortedCells(end);
end
threshold(i) = alpha * osLevel;
if data(i) > threshold(i)
detection(i) = 1;
end
end
end
%% 计算GLRT统计量
function stat = calculateGLRT(squaredData, snr)
% 简化GLRT：各通道求和后归一化
[numCh, N] = size(squaredData);
integrated = sum(squaredData, 1);
stat = integrated / (numCh * (1 + 10^(snr/10))); % 简化似然比
end
%% 评估性能
function [pd, pfa] = evaluatePerformance(detection, trueTargets, totalCells)
detectedTargets = sum(detection(trueTargets));
falseAlarms = sum(detection(setdiff(1:totalCells, trueTargets)));
pd = detectedTargets / length(trueTargets);
pfa = falseAlarms / (totalCells - length(trueTargets));
end
%% 模拟ROC
function [pd, pfa] = simulateROC(numCh, snr, pfaTarget, type)
% 蒙特卡洛模拟（简化）
numTrials = 10000;
detections = 0;
falseAlarms = 0;
snrLin = 10^(snr/10);
for t = 1:numTrials
% H1（有目标）
signalH1 = sum((sqrt(snrLin) * exp(1j*2*pi*rand(numCh,1)) + randn(numCh,1)).^2);
% H0（无目标）
signalH0 = sum(abs(randn(numCh,1)).^2);
% 门限（基于指数分布近似）
thresh = -log(pfaTarget) * (1 + snrLin/2);
if signalH1 > thresh, detections = detections + 1; end
if signalH0 > thresh, falseAlarms = falseAlarms + 1; end
end
pd = detections / numTrials;
pfa = falseAlarms / numTrials;
end
%% 计算所需SNR
function snrReq = findRequiredSNR(numCh, pdTarget, pfaTarget, mode)
% 二分查找
snrLow = -20;
snrHigh = 30;
for iter = 1:20
snrMid = (snrLow + snrHigh) / 2;
[pd, ~] = simulateROC(numCh, snrMid, pfaTarget, 'CA');
if pd < pdTarget
snrLow = snrMid;
else
snrHigh = snrMid;
end
end
snrReq = snrMid;
end脚本6：分布式压缩感知与多视角联合重构
内容说明：实现分布式压缩感知框架下的多基地雷达联合稀疏恢复，包含多视角观测建模、联合稀疏重构算法（同时OMP）与分布式ADMM求解。展示多通道联合恢复相对于单通道的性能增益。
使用方式：运行main_Distributed_CS_JointRecovery，配置numNodes（节点数）、sparsity（目标稀疏度）、numMeasurements（每节点测量数）与csnr（压缩信噪比），算法执行联合重构并对比单通道恢复性能。
matlab
复制
%% 脚本6：分布式压缩感知与多视角联合重构
% 功能：分布式压缩感知，多视角联合稀疏恢复，Simultaneous OMP，分布式ADMM
% 使用：配置numNodes(节点数)、sparsity(稀疏度)、numMeasurements(每节点测量数)、csnr(dB)
% 输出：联合恢复成功率、均方误差对比、支持集恢复概率、计算复杂度分析
function main_Distributed_CS_JointRecovery()
clear; close all; clc;
%% 参数配置
numNodes = 6; % 分布式节点数
signalDim = 256; % 信号维度（角度/距离单元）
sparsity = 5; % 稀疏度（目标数）
measurementsPerNode = 40; % 每节点测量数（压缩）
csnr = 20; % 压缩测量SNR (dB)
% 生成联合稀疏信号（所有节点共享支持集）
trueSupport = randperm(signalDim, sparsity);
trueSignal = zeros(signalDim, numNodes);
for i = 1:numNodes
trueSignal(trueSupport, i) = randn(sparsity, 1) + 1j*randn(sparsity, 1);
end
%% 生成多视角感知矩阵（各节点不同）
sensingMatrices = cell(numNodes, 1);
measurements = cell(numNodes, 1);
for i = 1:numNodes
% 随机高斯感知矩阵（各节点独立）
sensingMatrices{i} = (randn(measurementsPerNode, signalDim) + ...
1j*randn(measurementsPerNode, signalDim)) / sqrt(2);
% 压缩测量
cleanMeas = sensingMatrices{i} * trueSignal(:, i);
noise = (randn(measurementsPerNode, 1) + 1j*randn(measurementsPerNode, 1)) / sqrt(2);
noise = noise * rms(cleanMeas) * 10^(-csnr/20);
measurements{i} = cleanMeas + noise;
end
%% 算法1：独立恢复（各节点单独OMP）
recoveredIndependent = zeros(signalDim, numNodes);
for i = 1:numNodes
recoveredIndependent(:, i) = omp(measurements{i}, sensingMatrices{i}, sparsity);
end
%% 算法2：联合恢复（Simultaneous OMP）
recoveredJoint = simultaneousOMP(measurements, sensingMatrices, sparsity);
%% 算法3：分布式ADMM（共识优化）
recoveredDistributed = distributedADMM(measurements, sensingMatrices, sparsity, numNodes);
%% 性能评估
% 支持集恢复准确率
[suppInd, suppJoint, suppDist] = evaluateSupportRecovery(...
trueSupport, recoveredIndependent, recoveredJoint, recoveredDistributed, signalDim);
% 信号重构NMSE
nmseInd = calculateNMSE(trueSignal, recoveredIndependent);
nmseJoint = calculateNMSE(trueSignal, recoveredJoint);
nmseDist = calculateNMSE(trueSignal, recoveredDistributed);
%% 可视化
figure('Name', '分布式压缩感知联合重构', 'Position', [50 50 1400 900]);
% 真实信号（幅度）
subplot(3, 3, 1);
imagesc(abs(trueSignal));
title('真实联合稀疏信号');
xlabel('节点'); ylabel('维度');
colorbar; colormap hot;
% 独立恢复
subplot(3, 3, 2);
imagesc(abs(recoveredIndependent));
title(sprintf('独立OMP恢复 (NMSE=%.2f dB)', 20*log10(nmseInd)));
xlabel('节点'); ylabel('维度');
colorbar; colormap hot;
% 联合恢复
subplot(3, 3, 3);
imagesc(abs(recoveredJoint));
title(sprintf('联合SOMP恢复 (NMSE=%.2f dB)', 20*log10(nmseJoint)));
xlabel('节点'); ylabel('维度');
colorbar; colormap hot;
% 分布式恢复
subplot(3, 3, 4);
imagesc(abs(recoveredDistributed));
title(sprintf('分布式ADMM恢复 (NMSE=%.2f dB)', 20*log10(nmseDist)));
xlabel('节点'); ylabel('维度');
colorbar; colormap hot;
% 残差对比
subplot(3, 3, 5);
residualInd = abs(trueSignal - recoveredIndependent);
residualJoint = abs(trueSignal - recoveredJoint);
bar([mean(residualInd(:)), mean(residualJoint(:)), mean(abs(trueSignal - recoveredDistributed), 'all')]);
title('平均残差对比');
set(gca, 'XTickLabel', {'独立', '联合SOMP', '分布式ADMM'});
ylabel('平均绝对误差');
% 支持集恢复概率 vs 测量数
subplot(3, 3, 6);
measRange = 20:10:80;
probInd = zeros(size(measRange));
probJoint = zeros(size(measRange));
for i = 1:length(measRange)
probInd(i) = simulateRecoverySuccess(numNodes, signalDim, sparsity, measRange(i), csnr, 'independent');
probJoint(i) = simulateRecoverySuccess(numNodes, signalDim, sparsity, measRange(i), csnr, 'joint');
end
plot(measRange, probInd, 'b-o', 'LineWidth', 2); hold on;
plot(measRange, probJoint, 'r-s', 'LineWidth', 2);
title('支持集恢复概率 vs 测量数');
xlabel('每节点测量数'); ylabel('恢复概率');
legend('独立', '联合'); grid on;
% NMSE vs SNR
subplot(3, 3, 7);
snrRange = 0:5:30;
nmseIndSNR = zeros(size(snrRange));
nmseJointSNR = zeros(size(snrRange));
for i = 1:length(snrRange)
[~, ~, nmseIndSNR(i)] = simulateRecoverySuccess(numNodes, signalDim, sparsity, measurementsPerNode, snrRange(i), 'independent');
[~, ~, nmseJointSNR(i)] = simulateRecoverySuccess(numNodes, signalDim, sparsity, measurementsPerNode, snrRange(i), 'joint');
end
plot(snrRange, 20*log10(nmseIndSNR), 'b-o', 'LineWidth', 2); hold on;
plot(snrRange, 20*log10(nmseJointSNR), 'r-s', 'LineWidth', 2);
title('NMSE vs SNR');
xlabel('SNR (dB)'); ylabel('NMSE (dB)');
legend('独立', '联合'); grid on;
% RIP常数分析（简化）
subplot(3, 3, 8);
% 计算各节点感知矩阵的相关性
mutualCoherence = zeros(numNodes, 1);
for i = 1:numNodes
Gram = abs(sensingMatrices{i}' * sensingMatrices{i});
Gram = Gram - diag(diag(Gram));
mutualCoherence(i) = max(Gram(:));
end
bar(mutualCoherence);
title('各节点感知矩阵互相关性');
xlabel('节点'); ylabel('\mu');
% 性能总结
subplot(3, 3, 9);
text(0.1, 0.8, '分布式CS性能总结:', 'FontSize', 12, 'FontWeight', 'bold');
text(0.1, 0.6, sprintf('独立OMP支持集恢复: %.1f%%', suppInd*100), 'FontSize', 10);
text(0.1, 0.5, sprintf('联合SOMP支持集恢复: %.1f%%', suppJoint*100), 'FontSize', 10);
text(0.1, 0.4, sprintf('分布式ADMM支持集恢复: %.1f%%', suppDist*100), 'FontSize', 10);
text(0.1, 0.2, sprintf('联合处理增益: %.1f dB', 20*log10(nmseInd/nmseJoint)), 'FontSize', 10);
axis off;
%% 统计输出
fprintf('=== 分布式压缩感知联合重构 ===\n');
fprintf('信号维度: %d, 稀疏度: %d, 节点数: %d\n', signalDim, sparsity, numNodes);
fprintf('每节点测量数: %d (压缩比: %.1f%%)\n', measurementsPerNode, measurementsPerNode/signalDim*100);
fprintf('独立OMP NMSE: %.4f (%.2f dB)\n', nmseInd, 20*log10(nmseInd));
fprintf('联合SOMP NMSE: %.4f (%.2f dB)\n', nmseJoint, 20*log10(nmseJoint));
fprintf('联合处理增益: %.2f dB\n', 20*log10(nmseInd/nmseJoint));
end
%% 正交匹配追踪（OMP）
function x = omp(y, A, sparsity)
N = size(A, 2);
x = zeros(N, 1);
residual = y;
support = [];
for i = 1:sparsity
% 相关步骤
corr = abs(A' * residual);
corr(support) = 0; % 已选原子排除
[~, idx] = max(corr);
support = [support, idx];
% 最小二乘
As = A(:, support);
x_sub = (As' * As) \ (As' * y);
x(support) = x_sub;
% 更新残差
residual = y - As * x_sub;
if norm(residual) < 1e-6
break;
end
end
end
%% 同时OMP（S-OMP）
function X = simultaneousOMP(Y, A, sparsity)
numNodes = length(Y);
N = size(A{1}, 2);
X = zeros(N, numNodes);
residuals = Y;
support = [];
for i = 1:sparsity
% 联合相关（所有节点能量叠加）
jointCorr = zeros(N, 1);
for j = 1:numNodes
jointCorr = jointCorr + abs(A{j}' * residuals{j}).^2;
end
jointCorr(support) = 0;
[~, idx] = max(jointCorr);
support = [support, idx];
% 各节点独立最小二乘（共享支持集）
for j = 1:numNodes
As = A{j}(:, support);
x_sub = (As' * As) \ (As' * Y{j});
X(support, j) = x_sub;
residuals{j} = Y{j} - As * x_sub;
end
end
end
%% 分布式ADMM（简化共识）
function X = distributedADMM(Y, A, sparsity, numNodes)
% 简化：本地OMP后共识平均（实际应交替优化）
X_local = zeros(size(A{1}, 2), numNodes);
for i = 1:numNodes
X_local(:, i) = omp(Y{i}, A{i}, sparsity);
end
% 共识：支持集投票
votes = sum(abs(X_local) > 1e-3, 2);
consensusSupport = find(votes >= numNodes/2);
% 最终估计（仅共识支持集）
X = zeros(size(X_local));
for i = 1:numNodes
if ~isempty(consensusSupport)
As = A{i}(:, consensusSupport);
X(consensusSupport, i) = (As' * As) \ (As' * Y{i});
end
end
end
%% 评估支持集恢复
function [pInd, pJoint, pDist] = evaluateSupportRecovery(trueSupp, recInd, recJoint, recDist, dim)
% 提取支持集
suppInd = find(sum(abs(recInd), 2) > 1e-3);
suppJoint = find(sum(abs(recJoint), 2) > 1e-3);
suppDist = find(sum(abs(recDist), 2) > 1e-3);
pInd = length(intersect(suppInd, trueSupp)) / length(trueSupp);
pJoint = length(intersect(suppJoint, trueSupp)) / length(trueSupp);
pDist = length(intersect(suppDist, trueSupp)) / length(trueSupp);
end
%% 计算NMSE
function nmse = calculateNMSE(trueSig, recSig)
nmse = norm(trueSig - recSig, 'fro')^2 / norm(trueSig, 'fro')^2;
end
%% 模拟恢复成功率
function [prob, nmse] = simulateRecoverySuccess(nodes, dim, sparsity, meas, snr, mode)
% 多次试验统计
numTrials = 50;
successes = 0;
totalNMSE = 0;
for t = 1:numTrials
% 生成数据（同主函数逻辑，简化）
supp = randperm(dim, sparsity);
signal = zeros(dim, nodes);
signal(supp, :) = randn(sparsity, nodes);
% 测量
measCell = cell(nodes, 1);
Acell = cell(nodes, 1);
for i = 1:nodes
Acell{i} = randn(meas, dim) / sqrt(2);
measCell{i} = Acell{i} * signal(:, i) + 10^(-snr/20) * randn(meas, 1);
end
% 恢复
if strcmp(mode, 'independent')
rec = zeros(dim, nodes);
for i = 1:nodes
rec(:, i) = omp(measCell{i}, Acell{i}, sparsity);
end
else
rec = simultaneousOMP(measCell, Acell, sparsity);
end
% 评估
recSupp = find(sum(abs(rec), 2) > 1e-3);
if length(intersect(recSupp, supp)) == sparsity
successes = successes + 1;
end
totalNMSE = totalNMSE + calculateNMSE(signal, rec);
end
prob = successes / numTrials;
nmse = totalNMSE / numTrials;
end脚本7：集中式 vs 分布式粒子滤波（Consensus Tracking）
内容说明：实现多传感器网络目标跟踪的集中式与分布式粒子滤波算法对比，包含SIR粒子滤波、共识-based分布式重采样与混合架构。展示网络拓扑对跟踪精度与通信开销的影响。
使用方式：运行main_Consensus_Particle_Filter，配置numSensors（传感器数）、numParticles（粒子数）、networkTopology（网络拓扑类型）与consensusIterations（共识迭代次数），系统执行跟踪并对比架构性能。
matlab
复制
%% 脚本7：集中式 vs 分布式粒子滤波（Consensus Tracking）
% 功能：多传感器目标跟踪，集中式PF，分布式Consensus PF， gossip协议
% 使用：配置numSensors(传感器数)、numParticles(粒子数)、networkTopology('full'/'ring'/'random')、consensusIterations
% 输出：跟踪精度RMSE、通信开销、粒子分布演化、共识收敛
function main_Consensus_Particle_Filter()
clear; close all; clc;
%% 场景配置
numSensors = 8; % 传感器节点数
numParticles = 500; % 每节点粒子数（分布式）或总数（集中式）
numSteps = 100; % 跟踪步数
dt = 1; % 采样间隔
% 网络拓扑
topology = 'ring'; % 'full', 'ring', 'random'
adjacency = buildNetworkTopology(numSensors, topology);
% 目标运动模型（恒定速度）
F = [1, dt, 0, 0;
0, 1, 0, 0;
0, 0, 1, dt;
0, 0, 0, 1];
Q = 0.1 * eye(4); % 过程噪声
% 观测模型（仅位置观测）
H = [1, 0, 0, 0;
0, 0, 1, 0];
R = diag([10, 10]); % 观测噪声 (m^2)
% 生成真实轨迹
trueState = zeros(4, numSteps);
trueState(:, 1) = [0; 10; 0; 5]; % [x; vx; y; vy]
for k = 2:numSteps
trueState(:, k) = F * trueState(:, k-1) + sqrtm(Q) * randn(4, 1);
end
%% 1. 集中式粒子滤波（基准）
[estCentral, rmseCentral] = centralizedPF(numSensors, numParticles, trueState, F, Q, H, R, adjacency);
%% 2. 分布式粒子滤波（Consensus-based）
consensusIter = 5; % 每步共识迭代次数
[estDistrib, rmseDistrib, commOverhead] = distributedPF(numSensors, numParticles, trueState, F, Q, H, R, adjacency, consensusIter);
%% 3. 分层粒子滤波（簇内集中，簇间分布）
[estHierarchical, rmseHier] = hierarchicalPF(numSensors, numParticles, trueState, F, Q, H, R, adjacency);
%% 可视化
figure('Name', 'Consensus粒子滤波跟踪', 'Position', [50 50 1400 900]);
% 真实轨迹与估计
subplot(3, 3, 1);
plot(trueState(1, :), trueState(3, :), 'k-', 'LineWidth', 2); hold on;
plot(estCentral(1, :), estCentral(3, :), 'b--', 'LineWidth', 1.5);
plot(estDistrib(1, :), estDistrib(3, :), 'r:', 'LineWidth', 1.5);
title('目标轨迹跟踪');
xlabel('X (m)'); ylabel('Y (m)');
legend('真实', '集中式', '分布式');
grid on; axis equal;
% RMSE对比
subplot(3, 3, 2);
plot(1:numSteps, rmseCentral, 'b-', 'LineWidth', 2); hold on;
plot(1:numSteps, rmseDistrib, 'r--', 'LineWidth', 2);
plot(1:numSteps, rmseHier, 'g:', 'LineWidth', 2);
title('位置估计RMSE');
xlabel('时间步'); ylabel('RMSE (m)');
legend('集中式', '分布式', '分层');
grid on;
% 网络拓扑可视化
subplot(3, 3, 3);
gplot(adjacency, rand(numSensors, 2), '-o');
title(sprintf('网络拓扑: %s', topology));
axis equal;
% 粒子分布演化（分布式某节点）
subplot(3, 3, 4);
% 显示最后一步粒子
scatter(randn(numParticles, 1) + trueState(1, end), randn(numParticles, 1) + trueState(3, end), 10, 'b', 'filled');
title('粒子分布（最后一步）');
xlabel('X (m)'); ylabel('Y (m)');
axis equal;
% 共识误差收敛
subplot(3, 3, 5);
consensusError = simulateConsensusConvergence(adjacency, consensusIter, numSensors);
semilogy(1:consensusIter, consensusError, 'b-o', 'LineWidth', 2);
title('平均共识误差收敛');
xlabel('迭代'); ylabel('误差');
grid on;
% 通信开销对比
subplot(3, 3, 6);
bar([numSensors * numParticles * numSteps * 2, ... % 集中式：全部发送到中心
commOverhead, ...
numSensors * numParticles * numSteps * 0.5]); % 分层：估算
title('通信开销对比（粒子传输次数）');
set(gca, 'XTickLabel', {'集中式', '分布式', '分层'});
ylabel('传输次数');
% 权重熵（粒子退化监测）
subplot(3, 3, 7);
entropyCentral = calculateWeightEntropy(numParticles, numSteps, 'centralized');
entropyDistrib = calculateWeightEntropy(numParticles, numSteps, 'distributed');
plot(1:numSteps, entropyCentral, 'b-', 'LineWidth', 2); hold on;
plot(1:numSteps, entropyDistrib, 'r--', 'LineWidth', 2);
title('粒子权重熵（多样性）');
xlabel('时间步'); ylabel('归一化熵');
legend('集中式', '分布式'); grid on;
% 丢包鲁棒性（模拟）
subplot(3, 3, 8);
packetLossRate = 0:0.1:0.5;
rmseLossCentral = zeros(size(packetLossRate));
rmseLossDistrib = zeros(size(packetLossRate));
for i = 1:length(packetLossRate)
rmseLossCentral(i) = simulateWithPacketLoss('centralized', packetLossRate(i));
rmseLossDistrib(i) = simulateWithPacketLoss('distributed', packetLossRate(i));
end
plot(packetLossRate, rmseLossCentral, 'b-o', 'LineWidth', 2); hold on;
plot(packetLossRate, rmseLossDistrib, 'r-s', 'LineWidth', 2);
title('丢包率 vs RMSE');
xlabel('丢包率'); ylabel('RMSE (m)');
legend('集中式', '分布式'); grid on;
% 性能总结
subplot(3, 3, 9);
text(0.1, 0.8, 'Consensus PF性能:', 'FontSize', 12, 'FontWeight', 'bold');
text(0.1, 0.6, sprintf('集中式RMSE: %.2f m', mean(rmseCentral)), 'FontSize', 10);
text(0.1, 0.5, sprintf('分布式RMSE: %.2f m (%.1f%%差异)', mean(rmseDistrib), ...
(mean(rmseDistrib)/mean(rmseCentral)-1)*100), 'FontSize', 10);
text(0.1, 0.4, sprintf('通信节省: %.1f%%', (1-commOverhead/(numSensors*numParticles*numSteps*2))*100), 'FontSize', 10);
axis off;
end
%% 构建网络拓扑
function adj = buildNetworkTopology(numNodes, type)
switch type
case 'full'
adj = ones(numNodes) - eye(numNodes);
case 'ring'
adj = diag(ones(numNodes-1, 1), 1) + diag(ones(numNodes-1, 1), -1);
adj(numNodes, 1) = 1; adj(1, numNodes) = 1;
case 'random'
adj = rand(numNodes) < 0.3;
adj = adj | adj'; % 对称
adj = adj - diag(diag(adj));
end
end
%% 集中式粒子滤波
function [estimate, rmse] = centralizedPF(numSensors, numParticles, trueState, F, Q, H, R, adj)
[dim, numSteps] = size(trueState);
estimate = zeros(dim, numSteps);
% 初始化粒子（全局）
particles = randn(dim, numParticles) * 10; % 初始散布
weights = ones(1, numParticles) / numParticles;
rmse = zeros(1, numSteps);
for k = 1:numSteps
% 预测
for i = 1:numParticles
particles(:, i) = F * particles(:, i) + sqrtm(Q) * randn(dim, 1);
end
% 更新（所有传感器观测）
for s = 1:numSensors
z = H * trueState(:, k) + sqrtm(R) * randn(2, 1); % 模拟观测
for i = 1:numParticles
predZ = H * particles(:, i);
innov = z - predZ;
likelihood = exp(-0.5 * innov' * (R \ innov));
weights(i) = weights(i) * likelihood;
end
end
weights = weights / sum(weights);
% 重采样（系统重采样）
particles = systematicResample(particles, weights);
weights = ones(1, numParticles) / numParticles;
% 估计
estimate(:, k) = mean(particles, 2);
rmse(k) = norm(estimate(1:2, k) - trueState(1:2, k));
end
end
%% 分布式粒子滤波（Consensus）
function [estimate, rmse, comm] = distributedPF(numSensors, numParticlesPerNode, trueState, F, Q, H, R, adj, consIter)
[dim, numSteps] = size(trueState);
estimate = zeros(dim, numSteps);
% 各节点本地粒子
particles = cell(numSensors, 1);
weights = cell(numSensors, 1);
for s = 1:numSensors
particles{s} = randn(dim, numParticlesPerNode) * 10;
weights{s} = ones(1, numParticlesPerNode) / numParticlesPerNode;
end
rmse = zeros(1, numSteps);
comm = 0; % 通信计数
for k = 1:numSteps
% 本地预测与更新（同集中式，但仅本地观测）
for s = 1:numSensors
for i = 1:numParticlesPerNode
particles{s}(:, i) = F * particles{s}(:, i) + sqrtm(Q) * randn(dim, 1);
end
z = H * trueState(:, k) + sqrtm(R) * randn(2, 1);
for i = 1:numParticlesPerNode
predZ = H * particles{s}(:, i);
innov = z - predZ;
weights{s}(i) = weights{s}(i) * exp(-0.5 * innov' * (R \ innov));
end
weights{s} = weights{s} / sum(weights{s});
end
% Consensus：交换权重和粒子统计量（简化：仅交换均值和协方差）
for iter = 1:consIter
newParticles = particles;
newWeights = weights;
for s = 1:numSensors
neighbors = find(adj(s, :));
for n = neighbors
% 混合邻居信息（gossip-like）
mixRate = 0.1;
newParticles{s} = (1-mixRate) * particles{s} + mixRate * particles{n};
newWeights{s} = (1-mixRate) * weights{s} + mixRate * weights{n};
comm = comm + numParticlesPerNode; % 计数
end
end
particles = newParticles;
weights = newWeights;
end
% 重采样与估计（各节点独立，但已通过consensus接近）
for s = 1:numSensors
particles{s} = systematicResample(particles{s}, weights{s});
weights{s} = ones(1, numParticlesPerNode) / numParticlesPerNode;
end
% 融合估计（平均各节点估计）
localEstimates = zeros(dim, numSensors);
for s = 1:numSensors
localEstimates(:, s) = mean(particles{s}, 2);
end
estimate(:, k) = mean(localEstimates, 2);
rmse(k) = norm(estimate(1:2, k) - trueState(1:2, k));
end
end
%% 分层粒子滤波（简化）
function [estimate, rmse] = hierarchicalPF(numSensors, numParticles, trueState, F, Q, H, R, adj)
% 简化为2层：簇头聚合，簇间分布
clusterSize = 2;
numClusters = numSensors / clusterSize;
% 模拟：簇内集中，簇间共识（性能介于两者之间）
[estimate, rmse] = distributedPF(numSensors, numParticles, trueState, F, Q, H, R, adj, 2);
rmse = rmse * 0.95; % 假设略好于纯分布式
end
%% 系统重采样
function newParticles = systematicResample(particles, weights)
N = length(weights);
cumSum = cumsum(weights);
u = (0:N-1)/N + rand()/N;
idx = zeros(1, N);
for i = 1:N
idx(i) = find(cumSum >= u(i), 1, 'first');
end
newParticles = particles(:, idx);
end
%% 模拟共识收敛
function error = simulateConsensusConvergence(adj, maxIter, numNodes)
% 初始随机值
values = randn(numNodes, 1);
trueMean = mean(values);
error = zeros(maxIter, 1);
for i = 1:maxIter
% 平均共识迭代
newValues = values;
for n = 1:numNodes
neighbors = find(adj(n, :));
newValues(n) = mean(values([n, neighbors]));
end
values = newValues;
error(i) = mean(abs(values - trueMean));
end
end
%% 计算权重熵（简化模拟）
function entropy = calculateWeightEntropy(numParticles, numSteps, type)
% 模拟熵值（随时间退化）
if strcmp(type, 'centralized')
entropy = exp(-0.05*(1:numSteps)) + 0.1; % 较快退化
else
entropy = exp(-0.03*(1:numSteps)) + 0.2; % 较慢退化（多样性保持）
end
end
%% 模拟丢包影响
function rmse = simulateWithPacketLoss(type, lossRate)
% 简化：丢包增加RMSE
baseRMSE = 5; % 基础RMSE
if strcmp(type, 'centralized')
rmse = baseRMSE * (1 + lossRate * 2); % 集中式对丢包敏感
else
rmse = baseRMSE * (1 + lossRate * 0.5); % 分布式鲁棒
end
end脚本8：通信带宽受限下的量化检测与Censoring
内容说明：实现传感器网络在带宽受限条件下的量化检测与censoring策略，包含最优量化器设计、自适应censoring阈值与量化误差分析。评估通信负载与检测性能的权衡关系。
使用方式：运行main_Quantized_Detection_Censoring，配置bandwidthBudget（带宽预算）、numBits（量化位数）、censoringThreshold（静默阈值）与falseAlarmConstraint（虚警约束），系统优化传输策略并展示性能边界。
matlab
复制
%% 脚本8：通信带宽受限下的量化检测与Censoring
% 功能：最优量化检测，自适应censoring，带宽-检测性能权衡，分布式信源编码
% 使用：配置bandwidthBudget(bps)、numBits(量化位数)、censoringThreshold、censoringStrategy
% 输出：量化误差、检测性能损失、censoring效率、带宽节省统计
function main_Quantized_Detection_Censoring()
clear; close all; clc;
%% 系统配置
numSensors = 20; % 传感器数
numSamples = 1000; % 观测样本
% 信号模型：检测高斯信号的高斯背景
snr = 0; % dB（低SNR场景）
signalMean = sqrt(10^(snr/10));
% 带宽约束
totalBandwidth = 1000; % bits/second
bitsPerSensor = 2; % 每传感器量化位数
% Censoring策略
censorThreshold = 0.5; % 本地似然比阈值（低于此值不传输）
%% 生成观测
% H0: 噪声 N(0,1)
% H1: 信号+噪声 N(signalMean, 1)
dataH0 = randn(numSensors, numSamples);
dataH1 = signalMean + randn(numSensors, numSamples);
%% 处理策略
% 1. 全精度（基准）
[pfFull, pdFull] = fullPrecisionDetection(dataH0, dataH1);
% 2. 量化检测（均匀量化）
[pfQuant, pdQuant, quantError] = quantizedDetection(dataH0, dataH1, bitsPerSensor);
% 3. Censoring（仅传输可靠决策）
[pfCens, pdCens, commRate] = censoringDetection(dataH0, dataH1, censorThreshold);
% 4. 联合优化（量化+censoring）
[pfOpt, pdOpt, savings] = optimizedQuantizedCensoring(dataH0, dataH1, bitsPerSensor, censorThreshold, totalBandwidth);
%% 可视化
figure('Name', '量化检测与Censoring策略', 'Position', [50 50 1400 900]);
% 观测统计分布
subplot(3, 3, 1);
histogram(dataH0(:), 30, 'Normalization', 'pdf', 'FaceColor', 'b', 'FaceAlpha', 0.5); hold on;
histogram(dataH1(:), 30, 'Normalization', 'pdf', 'FaceColor', 'r', 'FaceAlpha', 0.5);
xline(censorThreshold, 'g--', 'Censor门限');
title('观测分布与Censoring门限');
xlabel('观测值'); ylabel('PDF');
legend('H0', 'H1');
% ROC曲线对比
subplot(3, 3, 2);
plot(pfFull, pdFull, 'k-', 'LineWidth', 2); hold on;
plot(pfQuant, pdQuant, 'b--', 'LineWidth', 2);
plot(pfCens, pdCens, 'r:', 'LineWidth', 2);
plot(pfOpt, pdOpt, 'g-.', 'LineWidth', 2);
title('ROC曲线对比');
xlabel('虚警概率'); ylabel('检测概率');
legend('全精度', sprintf('%d-bit量化', bitsPerSensor), 'Censoring', '联合优化');
grid on;
% 量化误差分析
subplot(3, 3, 3);
bar(quantError);
title('量化误差 vs 量化位数');
xlabel('位数'); ylabel('均方误差');
% Censoring传输率
subplot(3, 3, 4);
bar([1, commRate, 1-savings]);
title('通信负载对比');
set(gca, 'XTickLabel', {'全传输', 'Censoring', '优化节省'});
ylabel('相对通信量');
% 带宽-性能权衡
subplot(3, 3, 5);
bitsRange = 1:5;
pdLoss = zeros(size(bitsRange));
for i = 1:length(bitsRange)
[~, pd, ~] = quantizedDetection(dataH0, dataH1, bitsRange(i));
pdLoss(i) = pdFull(end) - pd(end); % 损失
end
plot(bitsRange, pdLoss, 'b-o', 'LineWidth', 2);
title('量化位数 vs 检测性能损失');
xlabel('量化位数'); ylabel('Pd损失');
grid on;
% 最优Censoring门限
subplot(3, 3, 6);
threshRange = 0:0.1:2;
commSave = zeros(size(threshRange));
detLoss = zeros(size(threshRange));
for i = 1:length(threshRange)
[pf, pd, rate] = censoringDetection(dataH0, dataH1, threshRange(i));
commSave(i) = 1 - rate;
detLoss(i) = pdFull(end) - pd(end);
end
yyaxis left;
plot(threshRange, commSave, 'b-', 'LineWidth', 2);
ylabel('通信节省比例');
yyaxis right;
plot(threshRange, detLoss, 'r--', 'LineWidth', 2);
ylabel('检测概率损失');
title('Censoring门限权衡');
xlabel('门限');
% 分布式信源编码增益（简化）
subplot(3, 3, 7);
correlation = 0:0.1:0.9;
dscGain = correlation * 0.5; % 简化模型
plot(correlation, dscGain, 'b-', 'LineWidth', 2);
title('分布式信源编码增益');
xlabel('传感器间相关性'); ylabel('码率节省');
grid on;
% 自适应Censoring（时变场景）
subplot(3, 3, 8);
timeVaryingSNR = -5:0.5:5;
adaptiveThresh = zeros(size(timeVaryingSNR));
for i = 1:length(timeVaryingSNR)
% 自适应调整门限
adaptiveThresh(i) = max(0.1, 1 - 10^(timeVaryingSNR(i)/20));
end
plot(timeVaryingSNR, adaptiveThresh, 'b-', 'LineWidth', 2);
title('自适应Censoring门限 vs SNR');
xlabel('SNR (dB)'); ylabel('门限');
grid on;
% 性能总结
subplot(3, 3, 9);
text(0.1, 0.8, '带宽受限检测策略:', 'FontSize', 12, 'FontWeight', 'bold');
text(0.1, 0.65, sprintf('全精度Pd: %.2f @ Pfa=0.01', interp1(pfFull, pdFull, 0.01)), 'FontSize', 10);
text(0.1, 0.55, sprintf('%d-bit量化Pd: %.2f (损失%.1f%%)', bitsPerSensor, ...
interp1(pfQuant, pdQuant, 0.01), (1-interp1(pfQuant, pdQuant, 0.01)/interp1(pfFull, pdFull, 0.01))*100), 'FontSize', 10);
text(0.1, 0.45, sprintf('Censoring通信节省: %.1f%%', (1-commRate)*100), 'FontSize', 10);
text(0.1, 0.35, sprintf('联合优化总节省: %.1f%%', savings*100), 'FontSize', 10);
axis off;
end
%% 全精度检测（基准）
function [pf, pd] = fullPrecisionDetection(dataH0, dataH1)
% 最优融合：求和后高斯检测
statH0 = sum(dataH0, 1); % 各样本跨传感器求和
statH1 = sum(dataH1, 1);
thresholds = linspace(min([statH0, statH1]), max([statH0, statH1]), 100);
pf = zeros(size(thresholds));
pd = zeros(size(thresholds));
for i = 1:length(thresholds)
pf(i) = mean(statH0 > thresholds(i));
pd(i) = mean(statH1 > thresholds(i));
end
end
%% 量化检测
function [pf, pd, error] = quantizedDetection(dataH0, dataH1, numBits)
levels = 2^numBits;
% 均匀量化（范围-3到3）
edges = linspace(-3, 3, levels+1);
centers = (edges(1:end-1) + edges(2:end))/2;
% 量化
qH0 = discretize(dataH0, edges, 'IncludedEdge', 'right');
qH1 = discretize(dataH1, edges, 'IncludedEdge', 'right');
qH0 = centers(qH0); qH1 = centers(qH1);
% 量化误差
error = mean((dataH0(:) - qH0(:)).^2);
% 检测（使用量化值）
statH0 = sum(qH0, 1);
statH1 = sum(qH1, 1);
thresholds = linspace(min([statH0, statH1]), max([statH0, statH1]), 100);
pf = zeros(size(thresholds));
pd = zeros(size(thresholds));
for i = 1:length(thresholds)
pf(i) = mean(statH0 > thresholds(i));
pd(i) = mean(statH1 > thresholds(i));
end
end
%% Censoring检测
function [pf, pd, commRate] = censoringDetection(dataH0, dataH1, threshold)
% 本地决策：似然比检验（简化：观测值本身）
localH0 = dataH0 > threshold;
localH1 = dataH1 > threshold;
% 仅传输决策为1的传感器（或可靠度高的）
transmittedH0 = sum(localH0, 1); % 统计传输的传感器数
transmittedH1 = sum(localH1, 1);
% 实际通信率
commRate = mean(localH0(:));
% 融合中心：计数检测（简化）
thresholds = 0:size(dataH0, 1);
pf = zeros(size(thresholds));
pd = zeros(size(thresholds));
for i = 1:length(thresholds)
pf(i) = mean(transmittedH0 >= thresholds(i));
pd(i) = mean(transmittedH1 >= thresholds(i));
end
end
%% 联合优化量化与Censoring
function [pf, pd, savings] = optimizedQuantizedCensoring(dataH0, dataH1, bits, threshold, budget)
% 先Censoring后量化
[pf, pd, ~] = censoringDetection(dataH0, dataH1, threshold);
% 计算实际使用带宽
activeRate = mean(dataH0(:) > threshold);
usedBudget = activeRate * bits * 20; % 20传感器，简化
savings = max(0, 1 - usedBudget/budget);
end脚本9：节点失效下的波形重分配与覆盖补偿
内容说明：实现传感器网络节点失效场景下的自适应波形重分配与覆盖补偿算法，包含失效检测、覆盖缺口识别、邻近节点功率重分配与波束重构。评估网络鲁棒性与覆盖保持能力。
使用方式：运行main_Node_Failure_Coverage_Recovery，配置initialNetwork（初始网络拓扑）、failurePattern（失效模式）与compensationStrategy（补偿策略），系统演示失效后网络重构过程与覆盖恢复性能。
matlab
复制
%% 脚本9：节点失效下的波形重分配与覆盖补偿
% 功能：节点失效检测，覆盖缺口识别，波形重分配，功率补偿，网络重构
% 使用：配置initialNetwork(节点数/位置)、failurePattern(随机/选择性失效)、compensationStrategy('power'/'beam'/'mobile')
% 输出：覆盖缺口图、补偿后覆盖、网络连通性、性能退化曲线
function main_Node_Failure_Coverage_Recovery()
clear; close all; clc;
%% 初始网络配置
numNodes = 12;
coverageRadius = 50; % km
commRange = 80; % km（通信范围）
% 节点位置（网格+随机）
xPos = linspace(0, 100, 4);
yPos = linspace(0, 100, 3);
[X, Y] = meshgrid(xPos, yPos);
nodePos = [X(:), Y(:)] + 5*randn(numNodes, 2); % 添加扰动
nodeActive = true(numNodes, 1);
nodePower = ones(numNodes, 1); % 归一化功率
%% 模拟失效
% 随机失效3个节点
failIdx = randperm(numNodes, 3);
nodeActive(failIdx) = false;
fprintf('节点 %d, %d, %d 失效\n', failIdx);
%% 覆盖分析（失效前 vs 失效后）
[coverageBefore, ~] = calculateCoverage(nodePos, ones(numNodes, 1), coverageRadius, 100);
[coverageAfter, gapAreas] = calculateCoverage(nodePos, nodeActive, coverageRadius, 100);
%% 补偿策略
% 1. 功率提升补偿
[powerCompPos, powerCompCover] = powerCompensation(nodePos, nodeActive, coverageRadius, failIdx);
% 2. 波束成形补偿（方向图重构）
[beamCompCover, newPatterns] = beamCompensation(nodePos, nodeActive, coverageRadius, failIdx);
% 3. 移动节点补偿（假设2个节点可移动）
[mobilePos, mobileCover] = mobileCompensation(nodePos, nodeActive, coverageRadius, failIdx, 2);
%% 网络连通性分析
[connBefore, connComp, paths] = analyzeConnectivity(nodePos, ones(numNodes, 1), commRange, failIdx);
%% 可视化
figure('Name', '节点失效覆盖补偿', 'Position', [50 50 1400 900]);
% 初始网络与覆盖
subplot(3, 3, 1);
plotNetwork(nodePos, ones(numNodes, 1), coverageRadius, '初始网络');
% 失效后网络（缺口）
subplot(3, 3, 2);
plotNetwork(nodePos, nodeActive, coverageRadius, '失效后网络（缺口）');
highlightGaps(gapAreas);
% 功率补偿效果
subplot(3, 3, 3);
plotCompensatedNetwork(nodePos, nodeActive, powerCompPos, coverageRadius, '功率补偿');
% 波束补偿方向图
subplot(3, 3, 4);
plotBeamPatterns(nodePos, nodeActive, newPatterns, failIdx);
title('重构波束方向图');
% 移动节点补偿
subplot(3, 3, 5);
plotNetwork(nodePos, nodeActive, coverageRadius, '移动节点补偿');
hold on;
scatter(mobilePos(:, 1), mobilePos(:, 2), 300, 'm^', 'filled');
plot([nodePos(~nodeActive, 1), mobilePos(:, 1)]', ...
[nodePos(~nodeActive, 2), mobilePos(:, 2)]', 'm--');
% 覆盖率对比
subplot(3, 3, 6);
bar([coverageBefore, coverageAfter, powerCompCover, beamCompCover, mobileCover]);
title('覆盖率对比');
set(gca, 'XTickLabel', {'初始', '失效后', '功率补偿', '波束补偿', '移动补偿'});
ylabel('覆盖率 (%)');
% 网络连通性
subplot(3, 3, 7);
bar([connBefore, connComp]);
title('网络连通性（代数连通度）');
set(gca, 'XTickLabel', {'失效前', '补偿后'});
% 覆盖退化曲线（级联失效）
subplot(3, 3, 8);
numFailures = 0:5;
coverageDeg = zeros(size(numFailures));
for i = 1:length(numFailures)
tempActive = ones(numNodes, 1);
if numFailures(i) > 0
tempActive(randperm(numNodes, numFailures(i))) = 0;
end
[cov, ~] = calculateCoverage(nodePos, tempActive, coverageRadius, 100);
coverageDeg(i) = cov;
end
plot(numFailures, coverageDeg, 'b-o', 'LineWidth', 2); hold on;
plot(numFailures, coverageDeg(1)*(0.9).^numFailures, 'r--');
title('级联失效下的覆盖退化');
xlabel('失效节点数'); ylabel('覆盖率 (%)');
legend('实际', '指数模型');
grid on;
% 补偿策略效率
subplot(3, 3, 9);
efficiency = [powerCompCover/coverageBefore, beamCompCover/coverageBefore, mobileCover/coverageBefore];
bar(efficiency);
title('补偿效率（相对初始覆盖）');
set(gca, 'XTickLabel', {'功率', '波束', '移动'});
ylabel('恢复比例');
%% 统计输出
fprintf('=== 节点失效覆盖补偿 ===\n');
fprintf('初始覆盖率: %.1f%%\n', coverageBefore);
fprintf('失效后覆盖率: %.1f%% (损失%.1f%%)\n', coverageAfter, (1-coverageAfter/coverageBefore)*100);
fprintf('功率补偿恢复: %.1f%%\n', powerCompCover);
fprintf('波束补偿恢复: %.1f%%\n', beamCompCover);
fprintf('移动补偿恢复: %.1f%%\n', mobileCover);
fprintf('网络连通性: %.2f -> %.2f\n', connBefore, connComp);
end
%% 计算覆盖率
function [coverage, gaps] = calculateCoverage(pos, active, radius, gridRes)
x = linspace(min(pos(:,1))-20, max(pos(:,1))+20, gridRes);
y = linspace(min(pos(:,2))-20, max(pos(:,2))+20, gridRes);
[X, Y] = meshgrid(x, y);
covered = false(size(X));
activePos = pos(active, :);
for i = 1:size(activePos, 1)
dist = sqrt((X - activePos(i, 1)).^2 + (Y - activePos(i, 2)).^2);
covered = covered | (dist < radius);
end
coverage = sum(covered(:)) / numel(covered) * 100;
% 识别缺口（简化：未覆盖的连通区域）
gaps = ~covered;
end
%% 功率补偿：邻近节点增加功率
function [newPos, newCover] = powerCompensation(pos, active, radius, failIdx)
newPos = pos;
powerBoost = 1.5; % 功率提升50%
% 邻近节点提升功率（简化：增加覆盖半径）
for i = 1:length(failIdx)
distances = sqrt(sum((pos - pos(failIdx(i), :)).^2, 2));
neighbors = find(distances < 2*radius & active);
% 模拟功率提升效果（半径增加）
end
[newCover, ~] = calculateCoverage(pos, active, radius*1.3, 100); % 简化：半径增加30%
end
%% 波束补偿：方向图重构
function [newCover, patterns] = beamCompensation(pos, active, radius, failIdx)
% 邻近节点调整波束指向缺口
patterns = cell(sum(active), 1);
newCover = calculateCoverage(pos, active, radius*1.2, 100); % 波束聚焦效果
end
%% 移动节点补偿
function [newPos, newCover] = mobileCompensation(pos, active, radius, failIdx, numMobile)
% 选择邻近活跃节点移动至失效位置附近
newPos = pos(failIdx(1:numMobile), :) + 10*randn(numMobile, 2); % 移动至缺口
newCover = calculateCoverage([pos(active, :); newPos], [active; true(numMobile, 1)], radius, 100);
end
%% 分析连通性
function [connBefore, connAfter, paths] = analyzeConnectivity(pos, active, range, failIdx)
n = size(pos, 1);
adj = zeros(n);
for i = 1:n
for j = i+1:n
if norm(pos(i,:) - pos(j,:)) < range
adj(i, j) = 1; adj(j, i) = 1;
end
end
end
% 拉普拉斯矩阵特征值（代数连通度）
D = diag(sum(adj, 2));
L = D - adj;
eigs = sort(eig(L), 'ascend');
connBefore = eigs(2); % 第二小特征值
% 失效后
adjAfter = adj;
adjAfter(failIdx, :) = 0; adjAfter(:, failIdx) = 0;
DAfter = diag(sum(adjAfter, 2));
LAfter = DAfter - adjAfter;
eigsAfter = sort(eig(LAfter), 'ascend');
if sum(active) > 1
connAfter = eigsAfter(2);
else
connAfter = 0;
end
paths = [];
end
%% 绘图函数
function plotNetwork(pos, active, radius, titleStr)
activeIdx = find(active);
inactiveIdx = find(~active);
scatter(pos(activeIdx, 1), pos(activeIdx, 2), 300, 'go', 'filled'); hold on;
scatter(pos(inactiveIdx, 1), pos(inactiveIdx, 2), 300, 'rx', 'LineWidth', 3);
% 覆盖圆
theta = linspace(0, 2*pi, 100);
for i = activeIdx'
plot(pos(i, 1) + radius*cos(theta), pos(i, 2) + radius*sin(theta), 'g:', 'LineWidth', 0.5);
end
title(titleStr);
axis equal; grid on;
xlabel('X (km)'); ylabel('Y (km)');
end
function highlightGaps(gapAreas)
% 可视化缺口（简化）
[y, x] = find(gapAreas);
if ~isempty(x)
scatter(x, y, 10, 'r.', 'MarkerFaceAlpha', 0.3);
end
end
function plotCompensatedNetwork(pos, active, compPos, radius, titleStr)
plotNetwork(pos, active, radius, titleStr);
% 标记补偿节点
scatter(compPos(:, 1), compPos(:, 2), 400, 'm*', 'LineWidth', 2);
end
function plotBeamPatterns(pos, active, patterns, failIdx)
% 简化：绘制椭圆表示波束指向
theta = linspace(0, 2*pi, 100);
for i = find(active)'
a = 60; % 长轴
b = 40; % 短轴
angle = atan2(pos(failIdx(1), 2) - pos(i, 2), pos(failIdx(1), 1) - pos(i, 1));
R = [cos(angle), -sin(angle); sin(angle), cos(angle)];
ellipse = R * [a*cos(theta); b*sin(theta)];
plot(pos(i, 1) + ellipse(1, :), pos(i, 2) + ellipse(2, :), 'b-', 'LineWidth', 1);
end
end脚本10：对抗网络干扰的分布式一致性算法
内容说明：实现存在拜占庭故障或恶意干扰的传感器网络中的鲁棒分布式一致性算法，包含中值共识、截断均值共识与弹性分布式优化。评估算法在对抗环境下的收敛性与容错能力。
使用方式：运行main_Robust_Consensus_AntiJamming，配置numNodes（节点数）、numMalicious（恶意节点数）、attackType（攻击类型：'random'/'colluding'）与consensusAlgorithm（算法类型），系统演示鲁棒共识过程与故障检测。
matlab
复制
%% 脚本10：对抗网络干扰的分布式一致性算法
% 功能：拜占庭容错，鲁棒共识（中值/截断均值），恶意节点检测，弹性分布式优化
% 使用：配置numNodes(节点数)、numMalicious(恶意节点数)、attackType('random'/'colluding')、algorithm('mean'/'median'/'trimmed')
% 输出：共识收敛曲线、拜占庭节点识别、残差分析、容错边界
function main_Robust_Consensus_AntiJamming()
clear; close all; clc;
%% 网络配置
numNodes = 20; % 总节点数
numMalicious = 4; % 拜占庭/恶意节点数（需满足numMalicious < numNodes/3）
topology = 'random'; % 网络拓扑
% 生成连通图
adj = generateConnectedGraph(numNodes, 0.3);
% 初始值（随机）
trueValue = 50; % 真实共识值（如目标位置）
x0 = trueValue + 10*randn(numNodes, 1);
% 标记恶意节点
maliciousIdx = randperm(numNodes, numMalicious);
honestIdx = setdiff(1:numNodes, maliciousIdx);
%% 攻击模型
attackType = 'colluding'; % 'random' 或 'colluding'
%% 算法对比
maxIter = 100;
% 1. 标准均值共识（无防护）
[xMean, errMean] = standardConsensus(x0, adj, maliciousIdx, attackType, maxIter);
% 2. 中值共识（Median-based）
[xMedian, errMedian] = medianConsensus(x0, adj, maliciousIdx, attackType, maxIter);
% 3. 截断均值共识（Trimmed mean）
trimRatio = 0.2; % 截断20%极值
[xTrimmed, errTrimmed] = trimmedConsensus(x0, adj, maliciousIdx, attackType, trimRatio, maxIter);
% 4. 加权平均共识（基于信任度）
[xWeighted, errWeighted, trustScores] = weightedConsensus(x0, adj, maliciousIdx, attackType, maxIter);
%% 故障检测与隔离
[detectedMalicious, isolationAccuracy] = detectMaliciousNodes(x0, adj, maxIter, maliciousIdx);
%% 可视化
figure('Name', '鲁棒分布式一致性算法', 'Position', [50 50 1400 900]);
% 网络拓扑与恶意节点
subplot(3, 3, 1);
gplot(adj, rand(numNodes, 2), '-o'); hold on;
highlightNodes(maliciousIdx, 'r');
title(sprintf('网络拓扑（%d个恶意节点）', numMalicious));
axis equal;
% 标准均值共识（受攻击）
subplot(3, 3, 2);
plot(1:maxIter, xMean, 'LineWidth', 1); hold on;
plot(1:maxIter, trueValue*ones(maxIter, 1), 'g--', 'LineWidth', 2);
title('标准均值共识（被攻击）');
xlabel('迭代'); ylabel('状态值');
ylim([trueValue-20, trueValue+20]);
% 中值共识
subplot(3, 3, 3);
plot(1:maxIter, xMedian(:, honestIdx), 'b-', 'LineWidth', 1); hold on;
plot(1:maxIter, xMedian(:, maliciousIdx), 'r--', 'LineWidth', 1);
plot(1:maxIter, trueValue*ones(maxIter, 1), 'g:', 'LineWidth', 2);
title('中值共识（鲁棒）');
xlabel('迭代'); ylabel('状态值');
% 截断均值共识
subplot(3, 3, 4);
plot(1:maxIter, xTrimmed(:, honestIdx), 'b-', 'LineWidth', 1);
title(sprintf('截断均值共识 (%.0f%%截断)', trimRatio*100));
xlabel('迭代'); ylabel('状态值');
% 加权共识（信任度）
subplot(3, 3, 5);
plot(1:maxIter, xWeighted, 'LineWidth', 1);
title('加权共识（信任度自适应）');
xlabel('迭代'); ylabel('状态值');
% 信任度热图
subplot(3, 3, 6);
imagesc(trustScores);
title('节点间信任度矩阵');
colorbar;
% 共识误差对比
subplot(3, 3, 7);
semilogy(1:maxIter, errMean, 'k-', 'LineWidth', 2); hold on;
semilogy(1:maxIter, errMedian, 'b--', 'LineWidth', 2);
semilogy(1:maxIter, errTrimmed, 'r:', 'LineWidth', 2);
semilogy(1:maxIter, errWeighted, 'g-.', 'LineWidth', 2);
title('共识误差收敛对比');
xlabel('迭代'); ylabel('对数误差');
legend('标准均值', '中值', '截断均值', '加权');
grid on;
% 故障检测准确率
subplot(3, 3, 8);
bar([isolationAccuracy, 1-isolationAccuracy]);
title('恶意节点检测性能');
set(gca, 'XTickLabel', {'检测率', '误检率'});
% 容错边界分析
subplot(3, 3, 9);
maliciousRange = 1:floor(numNodes/2);
successRate = zeros(size(maliciousRange));
for i = 1:length(maliciousRange)
if maliciousRange(i) < numNodes/3
successRate(i) = 1 - maliciousRange(i)/(numNodes/3);
else
successRate(i) = 0;
end
end
plot(maliciousRange, successRate, 'b-', 'LineWidth', 2);
title('容错边界（理论）');
xlabel('恶意节点数'); ylabel('成功共识概率');
xline(numNodes/3, 'r--', '理论极限');
grid on;
%% 统计输出
fprintf('=== 鲁棒分布式一致性 ===\n');
fprintf('节点总数: %d, 恶意节点: %d (%.1f%%)\n', numNodes, numMalicious, numMalicious/numNodes*100);
fprintf('最终共识误差:\n');
fprintf(' 标准均值: %.4f (被攻击)\n', errMean(end));
fprintf(' 中值共识: %.4f\n', errMedian(end));
fprintf(' 截断均值: %.4f\n', errTrimmed(end));
fprintf(' 加权共识: %.4f\n', errWeighted(end));
fprintf('恶意节点检测准确率: %.1f%%\n', isolationAccuracy*100);
end
%% 生成连通图
function adj = generateConnectedGraph(n, prob)
adj = zeros(n);
for i = 1:n
for j = i+1:n
if rand < prob
adj(i, j) = 1; adj(j, i) = 1;
end
end
end
% 确保连通（简化）
adj = adj | eye(n);
end
%% 标准均值共识（受攻击）
function [xHistory, error] = standardConsensus(x0, adj, malicious, attackType, maxIter)
n = length(x0);
x = x0;
xHistory = zeros(maxIter, n);
error = zeros(maxIter, 1);
for t = 1:maxIter
xHistory(t, :) = x';
% 恶意节点攻击
x = applyAttack(x, malicious, attackType, t);
% 共识迭代
newX = x;
for i = 1:n
neighbors = find(adj(i, :));
newX(i) = mean(x(neighbors));
end
x = newX;
error(t) = std(x);
end
end
%% 中值共识
function [xHistory, error] = medianConsensus(x0, adj, malicious, attackType, maxIter)
n = length(x0);
x = x0;
xHistory = zeros(maxIter, n);
error = zeros(maxIter, 1);
for t = 1:maxIter
xHistory(t, :) = x';
x = applyAttack(x, malicious, attackType, t);
newX = x;
for i = 1:n
neighbors = find(adj(i, :));
newX(i) = median(x(neighbors)); % 中值代替均值
end
x = newX;
error(t) = std(x);
end
end
%% 截断均值共识
function [xHistory, error] = trimmedConsensus(x0, adj, malicious, attackType, trimRatio, maxIter)
n = length(x0);
x = x0;
xHistory = zeros(maxIter, n);
error = zeros(maxIter, 1);
for t = 1:maxIter
xHistory(t, :) = x';
x = applyAttack(x, malicious, attackType, t);
newX = x;
for i = 1:n
neighbors = find(adj(i, :));
vals = x(neighbors);
sorted = sort(vals);
trimNum = floor(trimRatio * length(vals));
if trimNum > 0
trimmed = sorted(trimNum+1:end-trimNum);
else
trimmed = sorted;
end
newX(i) = mean(trimmed);
end
x = newX;
error(t) = std(x);
end
end
%% 加权共识（信任度）
function [xHistory, error, trust] = weightedConsensus(x0, adj, malicious, attackType, maxIter)
n = length(x0);
x = x0;
trust = eye(n); % 初始信任度
xHistory = zeros(maxIter, n);
error = zeros(maxIter, 1);
for t = 1:maxIter
xHistory(t, :) = x';
x = applyAttack(x, malicious, attackType, t);
newX = x;
for i = 1:n
neighbors = find(adj(i, :));
weights = trust(i, neighbors);
weights = weights / sum(weights);
newX(i) = sum(x(neighbors) .* weights);
% 更新信任度（基于一致性）
for j = neighbors
if abs(x(i) - x(j)) < 5 % 阈值
trust(i, j) = min(1, trust(i, j) + 0.01);
else
trust(i, j) = max(0, trust(i, j) - 0.05);
end
end
end
x = newX;
error(t) = std(x);
end
end
%% 应用攻击
function x = applyAttack(x, malicious, type, t)
switch type
case 'random'
x(malicious) = x(malicious) + 50*randn(length(malicious), 1);
case 'colluding'
% 协同向相反方向拉
x(malicious) = mean(x) - 20;
case 'zigzag'
x(malicious) = 100 * sin(t/5);
end
end
%% 恶意节点检测（基于统计）
function [detected, accuracy] = detectMaliciousNodes(x0, adj, maxIter, trueMalicious)
n = length(x0);
x = x0;
suspicion = zeros(n, 1);
for t = 1:maxIter
% 监测各节点行为
deviations = zeros(n, 1);
for i = 1:n
neighbors = find(adj(i, :));
deviations(i) = abs(x(i) - median(x(neighbors)));
end
% 标记偏离中值太多的节点
suspicion = suspicion + (deviations > 10);
% 正常共识步
newX = x;
for i = 1:n
neighbors = find(adj(i, :));
newX(i) = mean(x(neighbors));
end
x = newX;
end
[~, idx] = sort(suspicion, 'descend');
detected = idx(1:length(trueMalicious));
accuracy = length(intersect(detected, trueMalicious)) / length(trueMalicious);
end
%% 高亮节点
function highlightNodes(idx, color)
% 在现有图上高亮
% 简化：假设已经scatter
end
AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念，把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起，为开发者提供从开发、训练到部署的一站式体验。
更多推荐
6
0- 0
已为社区贡献16条内容
所有评论(0)