雷达信号处理 python实现（二）雷达信号的组成与幅度模型

雷达系统仿真 - 第二章：雷达信号的组成与幅度模型 版本: 1.0 功能: 实现公式(2.1-2.12)的完整仿真与可视化 涵盖: 信号模型、雷达距离方程、功率预算、损耗分析

#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ 雷达系统仿真 - 第二章：雷达信号的组成与幅度模型 版本: 1.1 (修复版) 功能: 实现公式(2.1-2.12)的完整仿真与可视化 涵盖: 信号模型、雷达距离方程、功率预算、损耗分析 """ import numpy as np import matplotlib.pyplot as plt import warnings warnings.filterwarnings('ignore') # 设置中文字体 plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans'] plt.rcParams['axes.unicode_minus'] = False def main(): print("="*70) print("雷达系统仿真 - 第二章：雷达信号的组成与幅度模型") print("="*70) # ==================== 系统参数设置 ==================== Pt = 1000 # 峰值发射功率 1kW G_dBi = 44 # 天线增益 44dBi G_linear = 10**(G_dBi/10) f0 = 10e9 # 载频 10GHz lambda_ = 3e8 / f0 # 波长 0.03m tau = 1e-6 # 脉冲宽度 1μs R_target = 10e3 # 目标距离 10km sigma_rcs = 100 # 雷达截面积 100m² Ls_dB = 6 # 系统损耗 6dB alpha = 0.1 # 大气损耗系数 dB/km print(f"\n【雷达系统参数】") print(f" 发射功率 Pt = {Pt/1e3:.1f} kW") print(f" 天线增益 G = {G_dBi} dBi (线性值 {G_linear:.0f})") print(f" 载频 f₀ = {f0/1e9:.1f} GHz, 波长 λ = {lambda_*100:.1f} cm") print(f" 脉冲宽度 τ = {tau*1e6:.1f} μs") print(f" 目标距离 R = {R_target/1e3:.1f} km, RCS σ = {sigma_rcs} m²") fig = plt.figure(figsize=(20, 16)) # ==================== 2.1 雷达信号的组成 ==================== print("\n█ 2.1 雷达信号的组成") print("-"*50) np.random.seed(42) phi = np.pi / 4 # 2.1 公式(2.1): 发射脉冲信号模型 ax1 = plt.subplot(4, 3, 1) Fs = 100e9 t_pulse = np.linspace(0, 2*tau, 1000) A = np.sqrt(2*Pt) a_t = np.where((t_pulse >= 0.25*tau) & (t_pulse <= 1.25*tau), A, 0) carrier = np.sin(2 * np.pi * f0 * t_pulse) x_t = a_t * carrier # 建议在作图专属区域定义一个展示频率 f_plot = 10e6 # 10MHz, 仅用于画图好看，1μs内呈现10个周期 # 修改 2.1 中的 carrier carrier = np.sin(2 * np.pi * f_plot * t_pulse) ax1.plot(t_pulse*1e6, x_t, 'b-', linewidth=1) ax1.fill_between(t_pulse*1e6, -A, A, where=(a_t>0), alpha=0.3, color='blue', label='脉冲包络') ax1.set_xlabel('时间 (μs)', fontsize=11) ax1.set_ylabel('幅度', fontsize=11) ax1.set_title(f'2.1 发射脉冲信号模型 (公式2.1)\nx(t)=a(t)sin[2πFₜt+θ(t)]', fontsize=11, fontweight='bold') ax1.grid(True, alpha=0.3) ax1.legend(fontsize=9) Ps = A**2 / 2 ax1.text(0.5, 0.9, f'瞬时功率 Ps={Ps/1e3:.1f}kW', transform=ax1.transAxes, fontsize=10, bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5)) print(f" 公式(2.1): 发射信号瞬时功率 Ps = A²/2 = {Ps/1e3:.1f} kW") # 2.2 公式(2.2): 接收回波信号模型 ax2 = plt.subplot(4, 3, 2) t0 = 2 * R_target / 3e8 t_receive = np.linspace(0, 2*t0 + 3*tau, 2000) k = 0.001 noise_power = 1e-12 noise = np.sqrt(noise_power) * np.random.randn(len(t_receive)) delay_idx = int(t0 * Fs / 100) a_t_delay = np.zeros_like(t_receive) # 将原有的脉冲起止逻辑修改为： pulse_start = np.argmin(np.abs(t_receive - t0)) # 找到最接近 t0 的时间点索引 pulse_end = pulse_start + int(tau / (t_receive[-1] - t_receive[0]) * len(t_receive)) if pulse_end < len(t_receive): a_t_delay[pulse_start:pulse_end] = k * A phi_t = np.pi/6 theta_t = 0 y_t = a_t_delay * np.sin(2 * np.pi * f_plot * (t_receive - t0) + theta_t + phi_t) + noise ax2.plot(t_receive*1e6, y_t, 'b-', linewidth=1, alpha=0.7, label='接收信号(含噪声)') ax2.plot(t_receive*1e6, a_t_delay, 'r--', linewidth=2, label='回波包络 k·a(t-t₀)') ax2.axvline(x=t0*1e6, color='g', linestyle=':', alpha=0.7, label=f'时延 t₀={t0*1e6:.1f}μs') ax2.set_xlabel('时间 (μs)', fontsize=11) ax2.set_ylabel('幅度', fontsize=11) ax2.set_title(f'2.2 接收回波信号模型\n时延t₀=2R/c={t0*1e6:.1f}μs', fontsize=11, fontweight='bold') ax2.legend(fontsize=8) ax2.grid(True, alpha=0.3) print(f" 公式(2.2): 目标距离 R = {R_target/1e3:.1f} km, 时延 t₀ = {t0*1e6:.2f} μs") # 距离分辨率 ax3 = plt.subplot(4, 3, 3) delta_R = 3e8 * tau / 2 ranges = np.linspace(0, 20e3, 1000) target1_pos = 5e3 target2_pos = target1_pos + delta_R response = np.exp(-(ranges - target1_pos)**2 / (2*(delta_R/3)**2)) response += np.exp(-(ranges - target2_pos)**2 / (2*(delta_R/3)**2)) ax3.plot(ranges/1e3, response, 'b-', linewidth=2) ax3.axvline(x=target1_pos/1e3, color='r', linestyle='--', alpha=0.5, label='目标1') ax3.axvline(x=target2_pos/1e3, color='g', linestyle='--', alpha=0.5, label='目标2') ax3.fill_betweenx([0, 1], target1_pos/1e3 - delta_R/2/1e3, target1_pos/1e3 + delta_R/2/1e3, alpha=0.2, color='red', label=f'分辨单元 ΔR={delta_R:.1f}m') ax3.set_xlabel('距离 (km)', fontsize=11) ax3.set_ylabel('回波强度', fontsize=11) ax3.set_title(f'2.1 距离分辨率验证\nΔR = cτ/2 = {delta_R:.1f}m', fontsize=11, fontweight='bold') ax3.legend(fontsize=9) ax3.grid(True, alpha=0.3) print(f" 距离分辨率: ΔR = cτ/2 = {delta_R:.1f} m") # ==================== 2.2 幅度模型与雷达距离方程 ==================== print("\n█ 2.2 幅度模型与雷达距离方程") print("-"*50) R_range = np.linspace(1e3, 50e3, 100) # 2.3 全向辐射功率密度 ax4 = plt.subplot(4, 3, 4) Q_omni = Pt / (4 * np.pi * R_range**2) ax4.loglog(R_range/1e3, Q_omni*1e3, 'b-', linewidth=2.5, label='全向辐射') ax4.set_xlabel('距离 R (km)', fontsize=11) ax4.set_ylabel('功率密度 (mW/m²)', fontsize=11) ax4.set_title('2.3 全向辐射功率密度\nQ = Pₜ/(4πR²)', fontsize=11, fontweight='bold') ax4.grid(True, alpha=0.3, which='both') ax4.legend() print(f" 公式(2.3): @10km处全向功率密度 = {Pt/(4*np.pi*10e3**2)*1e6:.2f} μW/m²") # 2.4 峰值发射功率密度 ax5 = plt.subplot(4, 3, 5) Qt = Pt * G_linear / (4 * np.pi * R_range**2) ax5.loglog(R_range/1e3, Qt, 'r-', linewidth=2.5, label=f'定向辐射 (G={G_dBi}dB)') ax5.loglog(R_range/1e3, Q_omni, 'b--', linewidth=1.5, alpha=0.7, label='全向辐射') ax5.set_xlabel('距离 R (km)', fontsize=11) ax5.set_ylabel('功率密度 (W/m²)', fontsize=11) ax5.set_title('2.4 峰值发射功率密度\nQₜ = PₜG/(4πR²)', fontsize=11, fontweight='bold') ax5.grid(True, alpha=0.3, which='both') ax5.legend() Qt_10km = Pt * G_linear / (4 * np.pi * 10e3**2) print(f" 公式(2.4): @10km处峰值功率密度 = {Qt_10km:.4f} W/m²") print(f" 增益使功率密度提高 {G_linear:.0f} 倍 ({G_dBi} dB)") # 2.5/2.6 后向散射功率与密度 (修复版) ax6 = plt.subplot(4, 3, 6) Pb = Pt * G_linear * sigma_rcs / (4 * np.pi * R_range**2) Qb = Pt * G_linear * sigma_rcs / ((4*np.pi)**2 * R_range**4) ax6_twin = ax6.twinx() # 使用plot + set_yscale替代semilogy以避免返回列表问题 ax6.plot(R_range/1e3, Pb*1e3, 'b-', linewidth=2.5) ax6_twin.plot(R_range/1e3, Qb*1e6, 'r--', linewidth=2.5) ax6.set_yscale('log') ax6_twin.set_yscale('log') ax6.set_xlabel('距离 R (km)', fontsize=11) ax6.set_ylabel('后向散射功率 Pb (mW)', fontsize=11, color='blue') ax6_twin.set_ylabel('功率密度 Qb (μW/m²)', fontsize=11, color='red') ax6.set_title('2.5/2.6 后向散射功率与密度', fontsize=11, fontweight='bold') ax6.grid(True, alpha=0.3, which='both') ax6.legend(['Pb=PₜGσ/(4πR²)'], loc='upper right') ax6_twin.legend(['Qb=PₜGσ/((4π)²R⁴)'], loc='center right') Pb_10km = Pt * G_linear * sigma_rcs / (4 * np.pi * 10e3**2) Qb_10km = Pt * G_linear * sigma_rcs / ((4*np.pi)**2 * 10e3**4) print(f" 公式(2.5): @10km处后向散射功率 Pb = {Pb_10km*1e6:.2f} μW") print(f" 公式(2.6): @10km处后向散射功率密度 Qb = {Qb_10km*1e9:.2f} nW/m²") # 2.7/2.8 接收功率计算 ax7 = plt.subplot(4, 3, 7) Ae = lambda_**2 * G_linear / (4 * np.pi) Pr = Pt * G_linear * Ae * sigma_rcs / ((4*np.pi)**2 * R_range**4) Pr_v2 = Pt * G_linear**2 * lambda_**2 * sigma_rcs / ((4*np.pi)**3 * R_range**4) ax7.loglog(R_range/1e3, Pr*1e12, 'b-', linewidth=2.5, label='接收功率 Pr') ax7.loglog(R_range/1e3, Pr_v2*1e12, 'r--', linewidth=2, alpha=0.7, label='公式(2.8)验证') ax7.set_xlabel('距离 R (km)', fontsize=11) ax7.set_ylabel('接收功率 (pW)', fontsize=11) ax7.set_title(f'2.7/2.8 雷达接收功率\nAₑ={Ae:.2f}m²', fontsize=11, fontweight='bold') ax7.grid(True, alpha=0.3, which='both') ax7.legend() Pr_10km = Pt * G_linear * Ae * sigma_rcs / ((4*np.pi)**2 * 10e3**4) print(f" 公式(2.7): @10km处接收功率 Pr = {Pr_10km*1e12:.2f} pW = {Pr_10km*1e9:.3f} nW") # 2.9/2.10 大气损耗 (修复版) ax8 = plt.subplot(4, 3, 8) La_dB = alpha * R_range / 500 La_linear = 10**(alpha * R_range / 5000) ax8_twin = ax8.twinx() ax8.plot(R_range/1e3, La_dB, 'b-', linewidth=2.5) ax8_twin.plot(R_range/1e3, La_linear, 'r--', linewidth=2.5) ax8.set_xlabel('距离 R (km)', fontsize=11) ax8.set_ylabel('损耗 (dB)', fontsize=11, color='blue') ax8_twin.set_ylabel('损耗因子 (线性)', fontsize=11, color='red') ax8.set_title(f'2.9/2.10 大气损耗模型\nα={alpha}dB/km', fontsize=11, fontweight='bold') ax8.grid(True, alpha=0.3) ax8.legend(['Lₐ(dB)=αR/500'], loc='center left') ax8_twin.legend(['Lₐ=10^(αR/5000)'], loc='center right') La_10km_dB = alpha * 10e3 / 500 print(f" 公式(2.9): @10km处大气损耗 = {La_10km_dB:.2f} dB") # 2.11 完整雷达距离方程 ax9 = plt.subplot(4, 3, 9) Ls_linear = 10**(Ls_dB/10) Pr_with_loss = Pr / (Ls_linear * 10**(alpha * R_range / 5000)) ax9.loglog(R_range/1e3, Pr*1e12, 'b--', linewidth=2, alpha=0.7, label='理想接收功率') ax9.loglog(R_range/1e3, Pr_with_loss*1e12, 'r-', linewidth=2.5, label=f'含损耗 (Ls={Ls_dB}dB)') ax9.set_xlabel('距离 R (km)', fontsize=11) ax9.set_ylabel('接收功率 (pW)', fontsize=11) ax9.set_title('2.11 雷达距离方程 (含损耗)\nPᵣ=PₜG²λ²σ/((4π)³R⁴LₛLₐ)', fontsize=11, fontweight='bold') ax9.grid(True, alpha=0.3, which='both') ax9.legend() Pr_loss_10km = Pr_10km / (Ls_linear * 10**(alpha * 10e3 / 5000)) print(f" 公式(2.11): 含损耗接收功率 = {Pr_loss_10km*1e12:.2f} pW") # 2.12 计算示例验证 ax10 = plt.subplot(4, 3, 10) Pr_calc = (Pt * G_linear**2 * lambda_**2 * sigma_rcs) / ((4*np.pi)**3 * R_target**4) categories = ['公式计算值', '文档示例值'] values = [Pr_calc*1e9, 3.07] bars = ax10.bar(categories, values, color=['steelblue', 'coral'], alpha=0.7, width=0.5) ax10.set_ylabel('接收功率 (nW)', fontsize=11) ax10.set_title('2.12 雷达方程计算验证\n(示例: Pt=1kW, G=44dB, λ=3cm, R=10km, σ=100m²)', fontsize=11, fontweight='bold') for bar, val in zip(bars, values): height = bar.get_height() ax10.text(bar.get_x() + bar.get_width()/2., height, f'{val:.2f}nW', ha='center', va='bottom', fontsize=10, fontweight='bold') ax10.grid(True, alpha=0.3, axis='y') print(f"\n 公式(2.12)验证:") print(f" 计算值: Pr = {Pr_calc*1e9:.3f} nW") print(f" 文档值: Pr = 3.07 nW") print(f" 误差: {abs(Pr_calc*1e9 - 3.07)/3.07 * 100:.1f}%") # RCS变化对接收功率的影响 ax11 = plt.subplot(4, 3, 11) rcs_values = np.array([0.01, 0.1, 1, 10, 100, 1000]) Pr_rcs = Pt * G_linear**2 * lambda_**2 * rcs_values / ((4*np.pi)**3 * R_target**4 * Ls_linear) ax11.semilogx(rcs_values, Pr_rcs*1e12, 'b-o', linewidth=2, markersize=8) ax11.axvline(x=sigma_rcs, color='r', linestyle='--', alpha=0.7, label=f'当前RCS={sigma_rcs}m²') ax11.set_xlabel('雷达截面积 σ (m²)', fontsize=11) ax11.set_ylabel('接收功率 (pW)', fontsize=11) ax11.set_title('2.8/2.11 RCS与接收功率关系\n@R=10km', fontsize=11, fontweight='bold') ax11.grid(True, alpha=0.3) ax11.legend() targets = [(0.01, '昆虫'), (1, '小型无人机'), (100, '战斗机'), (1000, '大型舰船')] for rcs_val, name in targets: idx = np.argmin(np.abs(rcs_values - rcs_val)) ax11.annotate(name, (rcs_values[idx], Pr_rcs[idx]*1e12), textcoords="offset points", xytext=(10, 0), fontsize=9) # 功率预算汇总表 ax12 = plt.subplot(4, 3, 12) ax12.axis('off') Pt_dbm = 10*np.log10(Pt*1000) Qt_db = 10*np.log10(Qt_10km) Qb_db = 10*np.log10(Qb_10km*1e6) Pr_dbm = 10*np.log10(Pr_10km*1000) summary_text = f""" 【雷达功率预算汇总 (@R={R_target/1e3:.0f}km)】 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 发射端: 峰值功率 Pₜ = {Pt/1e3:.1f} kW ({Pt_dbm:.1f} dBm) 天线增益 G = {G_dBi} dBi ({G_linear:.0f} 线性值) 波长 λ = {lambda_*100:.1f} cm (X波段) 传播过程: 峰值功率密度 Qₜ = {Qt_10km:.4f} W/m² 后向散射功率 Pᵦ = {Pb_10km*1e6:.2f} μW 后向散射密度 Qᵦ = {Qb_10km*1e9:.2f} nW/m² 接收端: 有效孔径 Aₑ = {Ae:.3f} m² 接收功率 Pᵣ = {Pr_10km*1e12:.2f} pW ({Pr_dbm:.1f} dBm) 损耗: 系统损耗 Lₛ = {Ls_dB} dB ({Ls_linear:.2f} 线性) 大气损耗 Lₐ = {La_10km_dB:.2f} dB 实际接收 = {Pr_loss_10km*1e12:.2f} pW 公式(2.12)验证: 计算值: {Pr_calc*1e9:.2f} nW ✓ 理论值: 3.07 nW ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ """ ax12.text(0.05, 0.5, summary_text, fontsize=9.5, family='monospace', verticalalignment='center', bbox=dict(boxstyle='round,pad=0.5', facecolor='wheat', alpha=0.4)) plt.tight_layout() plt.savefig('radar_chapter2_simulation.png', dpi=150, bbox_inches='tight') print("\n" + "="*70) print("仿真完成！可视化结果已保存为 radar_chapter2_simulation.png") print("="*70) plt.show() if __name__ == "__main__": main()