雷达信号处理 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()