Python算法
基于 Python 的单脉冲雷达导引头回波生成技术
基于 Python 的单脉冲雷达导引头回波高保真建模方法。通过构建目标 - 雷达几何模型、天线方向图合成及通道特性仿真，实现了从发射波形到接收基带信号的完整工具链。重点解决了多通道相干信号生成、系统误差注入（如幅度相位失配、RCS 起伏、角闪烁）、脉冲压缩集成等关键技术。提供了模块化代码实现和验证框架，支持算法研发、性能评估及故障诊断，为单脉冲测角算法研究提供标准化数据接口。
基于 Python 的单脉冲雷达导引头回波高保真建模方法。通过构建目标 - 雷达几何模型、天线方向图合成及通道特性仿真，实现了从发射波形到接收基带信号的完整工具链。重点解决了多通道相干信号生成、系统误差注入（如幅度相位失配、RCS 起伏、角闪烁）、脉冲压缩集成等关键技术。提供了模块化代码实现和验证框架，支持算法研发、性能评估及故障诊断，为单脉冲测角算法研究提供标准化数据接口。
本文系统阐述了单脉冲雷达导引头回波信号的高保真建模方法，通过 Python 实现了从目标 - 雷达几何建模、天线方向图合成、通道特性仿真到完整基带信号生成的完整工具链。文章重点解决了多通道相干信号生成、典型系统误差注入、脉冲压缩集成等关键技术难题，提供了可验证的仿真框架和模块化代码实现，为后续单脉冲测角算法研究提供了标准化的数据接口和测试基准。
单脉冲雷达技术因其在单次脉冲内即可获取角度信息的卓越特性，已成为现代精确制导、靶场测量、天文观测等领域的核心技术。与传统的圆锥扫描、顺序波束转换技术相比，单脉冲技术具备更高的角度测量精度、更强的抗干扰能力和更快的目标捕获速度。
然而，单脉冲处理算法的性能高度依赖于输入数据的质量。在实际工程应用中，和差通道间的幅度/相位不一致性、接收机非线性、目标角闪烁、多径效应等因素会严重恶化测角性能。算法研究人员若仅使用理想数学模型进行验证，将无法充分评估算法在实际系统中的鲁棒性。
教科书中的单脉冲公式往往基于若干理想假设：
现实系统则存在：
因此，构建一个能够可控地注入各类典型误差的高保真回波生成器，对于算法研发、系统性能预估和故障诊断具有不可替代的价值。
Python 在科学计算和快速原型开发中的独特优势使其成为雷达系统仿真的理想选择：
技术优势：
在本项目中的应用定位：
建立清晰的坐标系是精确建模的基础。需要定义以下关键坐标系：
# 坐标转换基础类定义
class CoordinateSystem:
"""坐标系管理与转换基类"""
def __init__(self):
# 定义各坐标系间的转换矩阵
self.T_inertial_to_body = np.eye(3) # 惯性系到弹体系
self.T_body_to_antenna = np.eye(3) # 弹体系到天线系
def inertial_to_antenna(self, pos_inertial):
"""将惯性系中的位置转换到天线坐标系"""
pos_body = np.dot(self.T_inertial_to_body, pos_inertial)
pos_antenna = np.dot(self.T_body_to_antenna, pos_body)
return pos_antenna
def compute_line_of_sight(self, target_pos, radar_pos):
"""计算视线向量及其在天线系中的角度"""
# 相对位置矢量
rel_pos = target_pos - radar_pos
R = np.linalg.norm(rel_pos) # 斜距
# 转换到天线坐标系
rel_pos_antenna = self.inertial_to_antenna(rel_pos)
# 计算方位角和俯仰角
azimuth = np.arctan2(rel_pos_antenna[1], rel_pos_antenna[0])
elevation = np.arcsin(rel_pos_antenna[2] / R)
return azimuth, elevation, R
在实际系统中，视线角的计算需要考虑：
和波束用于目标检测和测距，要求具有高增益和低副瓣特性。常用的数学模型包括：
高斯近似模型：
适用于快速仿真和理论分析
def sum_pattern_gaussian(theta, phi, beamwidth):
"""高斯型方向图函数
参数：
theta: 方位角偏移（弧度）
phi: 俯仰角偏移（弧度）
beamwidth: 3dB 波束宽度（弧度）
返回：
归一化的方向图增益
"""
# 将二维角度偏移转换为标量偏移
offset = np.sqrt(theta**2 + phi**2)
# 高斯函数参数计算
# 使在 beamwidth/2 处增益下降 3dB
sigma = beamwidth / (2 * np.sqrt(2 * np.log(2)))
pattern = np.exp(-offset**2 / (2 * sigma**2))
return pattern
实际阵列方向图模型：
基于阵列天线理论，更接近真实系统
def sum_pattern_array(theta, phi, freq, element_positions):
"""基于阵列理论的和波束方向图
参数：
theta, phi: 角度（弧度）
freq: 工作频率（Hz）
element_positions: N×3 阵列，每个阵元的位置
返回：
阵列方向图增益
"""
c = 3e8 # 光速
wavelength = c / freq
# 计算波数矢量
k = 2 * np.pi / wavelength
u = np.array([
np.sin(theta) * np.cos(phi),
np.sin(theta) * np.sin(phi),
np.cos(theta)
])
# 计算阵元间的相位差
phase_delays = np.dot(element_positions, u) * k
# 阵元加权（可设计为低副瓣加权）
weights = np.hamming(element_positions.shape[0]) # 海明加权
# 方向图计算
pattern = np.sum(weights * np.exp(1j * phase_delays))
return pattern / np.sum(weights) # 归一化
差波束用于角度测量，要求在波束中心形成零陷，两侧形成正负对称的响应。
一维差波束（适用于单平面测角）：
def difference_pattern_gaussian(theta, phi, beamwidth):
"""高斯型差波束方向图
采用奇对称设计，在波束中心为零
"""
offset = np.sqrt(theta**2 + phi**2)
sigma = beamwidth / (2 * np.sqrt(2 * np.log(2)))
# 奇函数设计：与角度成正比
pattern = (offset / (sigma * np.sqrt(2))) * np.exp(-offset**2 / (2 * sigma**2))
# 添加符号信息，使两侧响应反相
if theta != 0:
pattern = pattern * np.sign(theta)
return pattern
二维差波束（方位与俯仰独立）：
def dual_difference_pattern(theta, phi, beamwidth_az, beamwidth_el):
"""二维差波束，分别对应方位和俯仰通道
返回：
d_az: 方位差方向图
d_el: 俯仰差方向图
"""
sigma_az = beamwidth_az / (2 * np.sqrt(2 * np.log(2)))
sigma_el = beamwidth_el / (2 * np.sqrt(2 * np.log(2)))
# 方位差：对 theta 敏感，对 phi 不敏感
d_az = (theta / (sigma_az * np.sqrt(2))) * \
np.exp(-theta**2/(2*sigma_az**2) - phi**2/(2*sigma_el**2))
# 俯仰差：对 phi 敏感，对 theta 不敏感
d_el = (phi / (sigma_el * np.sqrt(2))) * \
np.exp(-theta**2/(2*sigma_az**2) - phi**2/(2*sigma_el**2))
return d_az, d_el
考虑一个简单的点目标场景，目标位于方位角θ_t、俯仰角φ_t 处：
class SignalPropagator:
"""信号传播与接收建模"""
def __init__(self, radar_params):
self.fc = radar_params['carrier_freq'] # 载频
self.c = 3e8
self.lambda_ = self.c / self.fc
def compute_channel_signals(self, theta_t, phi_t, R, target_rcs, sum_pattern, diff_pattern_az, diff_pattern_el):
"""计算和差通道的接收信号
参数：
theta_t, phi_t: 目标视线角
R: 斜距
target_rcs: 目标雷达截面积
各方向图函数
返回：
s_sum: 和通道复基带信号
s_diff_az: 方位差通道信号
s_diff_el: 俯仰差通道信号
"""
# 1. 基本传播衰减
# 自由空间路径损耗
path_loss = (self.lambda_**2) / ((4*np.pi)**3 * R**4)
# 2. 天线方向图调制
G_sum = sum_pattern(theta_t, phi_t) # 和波束增益
G_diff_az = diff_pattern_az(theta_t, phi_t) # 方位差增益
G_diff_el = diff_pattern_el(theta_t, phi_t) # 俯仰差增益
# 3. 目标反射特性（考虑 RCS 起伏）
sigma = self._apply_rcs_fluctuation(target_rcs)
# 4. 信号幅度计算（雷达方程）
A_sum = np.sqrt(path_loss * G_sum**2 * sigma)
A_diff_az = np.sqrt(path_loss * G_diff_az**2 * sigma)
A_diff_el = np.sqrt(path_loss * G_diff_el**2 * sigma)
# 5. 相位项（距离延迟和多普勒）
phase_delay = 2 * np.pi * 2 * R / self.lambda_ # 往返相位
doppler_phase = 0 # 多普勒项，需根据径向速度计算
total_phase = phase_delay + doppler_phase
# 6. 构建复基带信号
s_sum = A_sum * np.exp(1j * total_phase)
s_diff_az = A_diff_az * np.exp(1j * total_phase)
s_diff_el = A_diff_el * np.exp(1j * total_phase)
return s_sum, s_diff_az, s_diff_el
def _apply_rcs_fluctuation(self, mean_rcs):
"""应用 RCS 起伏模型（Swerling 模型）"""
# 此处可实现 Swerling I-IV 模型
return mean_rcs # 暂返回恒定值
对于运动目标，必须考虑多普勒频移的影响：
def apply_doppler_effect(signal, radial_velocity, fs, pulse_width):
"""应用多普勒频移到信号
参数：
signal: 输入信号
radial_velocity: 径向速度（m/s）
fs: 采样率
pulse_width: 脉冲宽度（s）
返回：
包含多普勒频移的信号
"""
# 计算多普勒频移
fd = 2 * radial_velocity / self.lambda_
# 生成时间向量
t = np.arange(len(signal)) / fs
# 应用相位调制
doppler_phase = 2 * np.pi * fd * t
signal_doppler = signal * np.exp(1j * doppler_phase)
return signal_doppler
线性调频信号是实现距离高分辨率的基石：
class LFMWaveform:
"""线性调频脉冲波形生成器"""
def __init__(self, bandwidth, pulse_width, fs, fc):
self.B = bandwidth # 带宽（Hz）
self.T = pulse_width # 脉冲宽度（s）
self.fs = fs # 采样率（Hz）
self.fc = fc # 载频（Hz）
self.K = self.B / self.T # 调频率（Hz/s）
def generate_pulse(self):
"""生成 LFM 脉冲"""
# 采样点数
N = int(self.T * self.fs)
t = np.arange(N) / self.fs - self.T/2
# 生成复基带 LFM 信号
baseband = np.exp(1j * np.pi * self.K * t**2)
# 上变频到载频（可选，取决于仿真需求）
carrier = np.exp(1j * 2 * np.pi * self.fc * t)
rf_signal = baseband * carrier
return baseband, rf_signal, t
def matched_filter(self):
"""生成匹配滤波器"""
# 匹配滤波器是发射信号的共轭反转
baseband, _, _ = self.generate_pulse()
matched_filter = np.conj(baseband[::-1])
return matched_filter
def ambiguity_function(self, delay_range, doppler_range):
"""计算模糊函数"""
baseband, _, t = self.generate_pulse()
N = len(baseband)
delays = np.linspace(-delay_range, delay_range, 100)
dopplers = np.linspace(-doppler_range, doppler_range, 100)
ambiguity = np.zeros((len(delays), len(dopplers)), dtype=complex)
for i, tau in enumerate(delays):
for j, fd in enumerate(dopplers):
# 延迟版本
delayed = np.roll(baseband, int(tau * self.fs))
# 多普勒调制
doppler_phase = 2 * np.pi * fd * t
doppler_modulated = delayed * np.exp(1j * doppler_phase)
# 互相关
ambiguity[i, j] = np.sum(baseband * np.conj(doppler_modulated))
return np.abs(ambiguity), delays, dopplers
脉冲压缩将长脉冲的能量压缩到窄脉冲，实现距离分辨率：
def pulse_compression(echo_signal, matched_filter):
"""执行脉冲压缩
参数：
echo_signal: 接收的回波信号
matched_filter: 匹配滤波器系数
返回：
compressed: 压缩后的信号
"""
# 频域卷积（更快）
N_echo = len(echo_signal)
N_filter = len(matched_filter)
N_fft = 2**np.ceil(np.log2(N_echo + N_filter - 1)).astype(int)
# 频域相乘
echo_fft = np.fft.fft(echo_signal, N_fft)
filter_fft = np.fft.fft(matched_filter, N_fft)
compressed_fft = echo_fft * filter_fft
compressed = np.fft.ifft(compressed_fft)[:N_echo]
return compressed
class MonopulseEchoSimulator:
"""单脉冲回波仿真器主类"""
def __init__(self, config_file=None):
# 加载配置
self.config = self._load_config(config_file)
# 初始化各子系统
self.geometry = GeometryEngine(self.config['geometry'])
self.antenna = AntennaModel(self.config['antenna'])
self.waveform = LFMWaveform(**self.config['waveform'])
self.target = TargetModel(self.config['target'])
self.channel = ChannelModel(self.config['channel'])
self.effects = EffectsManager(self.config['effects'])
# 结果缓存
self.results = {}
def _load_config(self, config_file):
"""加载仿真配置文件"""
if config_file is None:
# 默认配置
config = {
'simulation': {
'duration': 10.0, # 仿真时长（秒）
'fs': 100e6, # 采样率（Hz）
'prf': 1000, # 脉冲重复频率（Hz）
},
'radar': {
'fc': 10e9, # 载频（Hz）
'beamwidth_az': 3.0, # 方位波束宽度（度）
'beamwidth_el': 3.0, # 俯仰波束宽度（度）
},
# ... 其他配置项
}
else:
with open(config_file, 'r') as f:
config = json.load(f)
return config
def run_simulation(self):
"""执行完整仿真流程"""
print("开始单脉冲回波仿真...")
# 步骤 1: 目标轨迹生成
print("生成目标轨迹...")
target_trajectory = self.target.generate_trajectory(
self.config['simulation']['duration']
)
# 步骤 2: 对每个脉冲时刻进行处理
prf = self.config['simulation']['prf']
num_pulses = int(self.config['simulation']['duration'] * prf)
# 初始化数据存储
sum_channel_data = []
diff_az_data = []
diff_el_data = []
for pulse_idx in range(num_pulses):
# 计算当前时间
t = pulse_idx / prf
# 步骤 3: 计算几何关系
az, el, R = self.geometry.compute_angles(t, target_trajectory)
# 步骤 4: 生成发射波形
tx_pulse, _, _ = self.waveform.generate_pulse()
# 步骤 5: 计算和差通道响应
s_sum, s_diff_az, s_diff_el = self.antenna.compute_response(az, el)
# 步骤 6: 应用传播效应
echo_sum = self._apply_propagation(tx_pulse, s_sum, R)
echo_diff_az = self._apply_propagation(tx_pulse, s_diff_az, R)
echo_diff_el = self._apply_propagation(tx_pulse, s_diff_el, R)
# 步骤 7: 应用系统误差和环境影响
echo_sum = self.effects.apply_all(echo_sum, channel_type='sum')
echo_diff_az = self.effects.apply_all(echo_diff_az, channel_type='diff_az')
echo_diff_el = self.effects.apply_all(echo_diff_el, channel_type='diff_el')
# 存储数据
sum_channel_data.append(echo_sum)
diff_az_data.append(echo_diff_az)
diff_el_data.append(echo_diff_el)
# 进度显示
if pulse_idx % 100 == 0:
print(f"进度：{pulse_idx}/{num_pulses}")
# 步骤 8: 整理输出数据
self.results = {
'sum_channel': np.array(sum_channel_data),
'diff_az': np.array(diff_az_data),
'diff_el': np.array(diff_el_data),
'time_axis': np.arange(num_pulses) / prf,
'parameters': self.config
}
print("仿真完成！")
return self.results
def _apply_propagation(self, tx_signal, antenna_gain, R):
"""应用传播效应"""
# 距离衰减
attenuation = 1.0 / (R**2) # 简化模型
# 延迟效应
delay_samples = int(2 * R / 3e8 * self.config['simulation']['fs'])
# 应用衰减和延迟
echo = tx_signal * antenna_gain * attenuation
echo = np.roll(echo, delay_samples)
# 补零处理延迟溢出
if delay_samples > 0:
echo[:delay_samples] = 0
return echo
本架构采用模块化设计，具备以下优点：
单脉冲处理依赖于和差通道信号的相干性。相干性破坏将导致测角性能急剧下降。保持相干性的关键：
class CoherentSignalGenerator:
"""相干信号生成器"""
def generate_coherent_channels(self, base_signal, azimuth, elevation):
"""从基础信号生成相干的多通道信号"""
# 1. 生成公共的随机相位（模拟目标反射的随机初相）
common_phase = np.random.uniform(0, 2*np.pi)
# 2. 计算各通道的幅度调制因子
sum_factor = self.antenna.sum_pattern(azimuth, elevation)
diff_az_factor = self.antenna.diff_az_pattern(azimuth, elevation)
diff_el_factor = self.antenna.diff_el_pattern(azimuth, elevation)
# 3. 应用相同的相位噪声源（保持相位关系）
phase_noise = self._generate_phase_noise(len(base_signal))
# 4. 生成各通道信号
sum_signal = base_signal * sum_factor * np.exp(1j*common_phase) * phase_noise
diff_az_signal = base_signal * diff_az_factor * np.exp(1j*common_phase) * phase_noise
diff_el_signal = base_signal * diff_el_factor * np.exp(1j*common_phase) * phase_noise
return {
'sum': sum_signal,
'diff_az': diff_az_signal,
'diff_el': diff_el_signal,
'common_phase': common_phase
}
def _generate_phase_noise(self, length, f0=100e3, phase_noise_dbc=-80):
"""生成相位噪声，模拟振荡器不稳定性"""
t = np.arange(length) / self.fs
# 生成 1/f 相位噪声
f = np.fft.fftfreq(length, 1/self.fs)
f[0] = 1e-6 # 避免除以零
# 相位噪声功率谱密度
psd = 10**(phase_noise_dbc/10) / (np.abs(f)/f0 + 1e-6)
# 生成随机相位
random_phase = np.random.randn(length) + 1j*np.random.randn(length)
random_phase_fft = np.fft.fft(random_phase)
# 应用功率谱形状
phase_noise_fft = random_phase_fft * np.sqrt(psd)
phase_noise = np.fft.ifft(phase_noise_fft)
# 归一化并转换为相位因子
phase_factor = np.exp(1j * np.angle(phase_noise))
return phase_factor
def evaluate_channel_coherence(channel_data):
"""评估通道间的相干性"""
sum_sig = channel_data['sum']
diff_az_sig = channel_data['diff_az']
diff_el_sig = channel_data['diff_el']
# 计算互相关系数
corr_sum_diff_az = np.corrcoef(sum_sig, diff_az_sig)[0, 1]
corr_sum_diff_el = np.corrcoef(sum_sig, diff_el_sig)[0, 1]
# 计算相位一致性
phase_diff_az = np.angle(sum_sig * np.conj(diff_az_sig))
phase_std_az = np.std(phase_diff_az)
phase_diff_el = np.angle(sum_sig * np.conj(diff_el_sig))
phase_std_el = np.std(phase_diff_el)
coherence_metrics = {
'correlation_sum_diff_az': abs(corr_sum_diff_az),
'correlation_sum_diff_el': abs(corr_sum_diff_el),
'phase_std_az_deg': np.degrees(phase_std_az),
'phase_std_el_deg': np.degrees(phase_std_el),
'is_coherent': (abs(corr_sum_diff_az) > 0.9 and abs(corr_sum_diff_el) > 0.9 and phase_std_az < 0.1 and phase_std_el < 0.1)
}
return coherence_metrics
class ChannelErrorModel:
"""通道误差建模"""
def __init__(self, config):
# 幅度不平衡（dB 转换为线性比）
self.gain_imbalance_db = config.get('gain_imbalance_db', 0.5)
self.gain_imbalance_linear = 10**(self.gain_imbalance_db/20)
# 相位不平衡（度转换为弧度）
self.phase_imbalance_deg = config.get('phase_imbalance_deg', 3.0)
self.phase_imbalance_rad = np.radians(self.phase_imbalance_deg)
# 时间延迟不平衡（秒）
self.delay_imbalance_sec = config.get('delay_imbalance_sec', 10e-12)
# I/Q 不平衡
self.iq_imbalance_db = config.get('iq_imbalance_db', 0.3)
self.iq_phase_imbalance_deg = config.get('iq_phase_imbalance_deg', 2.0)
def apply_to_channel(self, signal, channel_type, fs):
"""应用误差到特定通道"""
# 复制信号避免修改原数据
corrupted_signal = signal.copy()
# 1. 幅度误差
if channel_type == 'diff_az':
gain_error = self.gain_imbalance_linear
elif channel_type == 'diff_el':
gain_error = 1.0 / self.gain_imbalance_linear
else: # sum 通道
gain_error = 1.0
corrupted_signal *= gain_error
# 2. 相位误差
if channel_type == 'diff_az':
phase_error = np.exp(1j * self.phase_imbalance_rad)
elif channel_type == 'diff_el':
phase_error = np.exp(-1j * self.phase_imbalance_rad)
else:
phase_error = 1.0
corrupted_signal *= phase_error
# 3. 时间延迟误差
if self.delay_imbalance_sec != 0:
delay_samples = int(self.delay_imbalance_sec * fs)
if delay_samples != 0:
corrupted_signal = np.roll(corrupted_signal, delay_samples)
if delay_samples > 0:
corrupted_signal[:delay_samples] = 0
else:
corrupted_signal[delay_samples:] = 0
# 4. I/Q 不平衡（对于基带信号）
if np.iscomplexobj(corrupted_signal):
corrupted_signal = self._apply_iq_imbalance(corrupted_signal)
return corrupted_signal
def _apply_iq_imbalance(self, complex_signal):
"""应用 I/Q 不平衡"""
I = np.real(complex_signal)
Q = np.imag(complex_signal)
# 增益不平衡
i_gain = 10**(self.iq_imbalance_db/20)
q_gain = 1.0 / i_gain
I_imbalanced = I * i_gain
Q_imbalanced = Q * q_gain
# 相位不平衡
phi = np.radians(self.iq_phase_imbalance_deg)
I_corrected = I_imbalanced
Q_corrected = Q_imbalanced * np.cos(phi) + I_imbalanced * np.sin(phi)
return I_corrected + 1j * Q_corrected
def generate_error_report(self):
"""生成误差配置报告"""
report = f"""
通道误差配置：
----------------------------
幅度不平衡：{self.gain_imbalance_db:.2f} dB (线性比：{self.gain_imbalance_linear:.4f})
相位不平衡：{self.phase_imbalance_deg:.2f} 度 ({self.phase_imbalance_rad:.4f} 弧度)
延迟不平衡：{self.delay_imbalance_sec*1e12:.2f} ps
I/Q 幅度不平衡：{self.iq_imbalance_db:.2f} dB
I/Q 相位不平衡：{self.iq_phase_imbalance_deg:.2f} 度
"""
return report
def analyze_error_impact(error_model, angle_range=np.radians(np.linspace(-5, 5, 100))):
"""分析误差对测角 S 曲线的影响"""
ideal_s_curve = []
corrupted_s_curve = []
for theta in angle_range:
# 理想情况
sum_sig_ideal = np.exp(1j * 0) # 简化的和信号
diff_sig_ideal = theta * np.exp(1j * 0) # 简化的差信号
s_ideal = np.real(diff_sig_ideal * np.conj(sum_sig_ideal))
ideal_s_curve.append(s_ideal)
# 加误差情况
diff_sig_corrupted = error_model.apply_to_channel(
diff_sig_ideal, 'diff_az', fs=1e6
)
s_corrupted = np.real(diff_sig_corrupted * np.conj(sum_sig_ideal))
corrupted_s_curve.append(s_corrupted)
# 计算误差指标
ideal_slope = np.polyfit(angle_range[:20], ideal_s_curve[:20], 1)[0]
corrupted_slope = np.polyfit(angle_range[:20], corrupted_s_curve[:20], 1)[0]
slope_error_percent = abs(corrupted_slope - ideal_slope) / ideal_slope * 100
zero_crossing_ideal = np.interp(0, ideal_s_curve, angle_range)
zero_crossing_corrupted = np.interp(0, corrupted_s_curve, angle_range)
bias_error = zero_crossing_corrupted - zero_crossing_ideal
return {
'ideal_s_curve': ideal_s_curve,
'corrupted_s_curve': corrupted_s_curve,
'angle_range': angle_range,
'slope_error_percent': slope_error_percent,
'bias_error_rad': bias_error,
'bias_error_deg': np.degrees(bias_error)
}
class SwerlingModel:
"""Swerling 起伏模型实现"""
def __init__(self, model_type='Swerling1', mean_rcs=1.0):
self.model_type = model_type
self.mean_rcs = mean_rcs
# 模型参数
self.models = {
'Swerling1': {'pdf': 'exponential', 'correlation': 'pulse_to_pulse'},
'Swerling2': {'pdf': 'exponential', 'correlation': 'independent'},
'Swerling3': {'pdf': 'chi_square_4', 'correlation': 'pulse_to_pulse'},
'Swerling4': {'pdf': 'chi_square_4', 'correlation': 'independent'},
'Swerling0': {'pdf': 'constant', 'correlation': 'constant'} # 非起伏
}
def generate_rcs_sequence(self, num_pulses, prf):
"""生成 RCS 起伏序列"""
model_info = self.models.get(self.model_type, self.models['Swerling1'])
if model_info['pdf'] == 'exponential':
# Swerling I/II：指数分布
rcs_values = np.random.exponential(self.mean_rcs, num_pulses)
# 对于 Swerling I，脉冲间相关
if model_info['correlation'] == 'pulse_to_pulse':
# 应用一阶马尔可夫过程
alpha = 0.9 # 相关系数
rcs_values_correlated = np.zeros(num_pulses)
rcs_values_correlated[0] = rcs_values[0]
for i in range(1, num_pulses):
rcs_values_correlated[i] = (
alpha * rcs_values_correlated[i-1] + (1-alpha) * rcs_values[i]
)
rcs_values = rcs_values_correlated
elif model_info['pdf'] == 'chi_square_4':
# Swerling III/IV：4 自由度的卡方分布
# 通过两个独立高斯变量的平方和实现
gauss1 = np.random.randn(num_pulses)
gauss2 = np.random.randn(num_pulses)
rcs_values = (gauss1**2 + gauss2**2) * self.mean_rcs / 2
if model_info['correlation'] == 'pulse_to_pulse':
# 类似处理相关
alpha = 0.9
rcs_values_correlated = np.zeros(num_pulses)
rcs_values_correlated[0] = rcs_values[0]
for i in range(1, num_pulses):
rcs_values_correlated[i] = (
alpha * rcs_values_correlated[i-1] + (1-alpha) * rcs_values[i]
)
rcs_values = rcs_values_correlated
else:
# 恒定 RCS
rcs_values = np.ones(num_pulses) * self.mean_rcs
return rcs_values
def plot_rcs_statistics(self, num_samples=10000):
"""绘制 RCS 统计特性"""
rcs_samples = self.generate_rcs_sequence(num_samples, 1000)
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
# 1. 时间序列
axes[0, 0].plot(rcs_samples[:200])
axes[0, 0].set_xlabel('脉冲序号')
axes[0, 0].set_ylabel('RCS')
axes[0, 0].set_title(f'{self.model_type} RCS 时间序列')
axes[0, 0].grid(True)
# 2. 直方图
axes[0, 1].hist(rcs_samples, bins=50, density=True, alpha=0.7)
# 绘制理论 PDF
x = np.linspace(0, np.max(rcs_samples), 100)
if self.models[self.model_type]['pdf'] == 'exponential':
theoretical_pdf = np.exp(-x/self.mean_rcs) / self.mean_rcs
axes[0, 1].plot(x, theoretical_pdf, 'r-', linewidth=2, label='理论指数分布')
axes[0, 1].set_xlabel('RCS')
axes[0, 1].set_ylabel('概率密度')
axes[0, 1].set_title('RCS 概率分布')
axes[0, 1].legend()
axes[0, 1].grid(True)
# 3. 自相关函数
autocorr = np.correlate(rcs_samples, rcs_samples, mode='full')
autocorr = autocorr[len(autocorr)//2:]
autocorr = autocorr[:100] / autocorr[0]
axes[1, 0].plot(autocorr)
axes[1, 0].set_xlabel('滞后')
axes[1, 0].set_ylabel('自相关系数')
axes[1, 0].set_title('RCS 自相关函数')
axes[1, 0].grid(True)
# 4. 累积分布函数
sorted_rcs = np.sort(rcs_samples)
cdf = np.arange(1, len(sorted_rcs)+1) / len(sorted_rcs)
axes[1, 1].plot(sorted_rcs, cdf)
axes[1, 1].set_xlabel('RCS')
axes[1, 1].set_ylabel('累积概率')
axes[1, 1].set_title('RCS 累积分布函数')
axes[1, 1].grid(True)
plt.tight_layout()
return fig
class GlintModel:
"""目标角闪烁效应建模"""
def __init__(self, target_length=10.0, wavelength=0.03):
"""
参数：
target_length: 目标特征长度（米）
wavelength: 波长（米）
"""
self.target_length = target_length
self.wavelength = wavelength
# 闪烁频率与目标旋转相关
self.rotation_rate = 1.0 # 弧度/秒
def generate_glint_error(self, num_samples, fs, target_aspect):
"""生成角闪烁误差序列"""
# 闪烁谱的带宽与目标尺寸相关
glint_bandwidth = 2 * self.target_length / self.wavelength * self.rotation_rate
# 生成带限高斯噪声
t = np.arange(num_samples) / fs
# 1. 生成白噪声
white_noise = np.random.randn(num_samples)
# 2. 设计带通滤波器（模拟闪烁谱）
from scipy.signal import butter, filtfilt
nyquist = fs / 2
lowcut = 0.1 * glint_bandwidth / nyquist
highcut = 2.0 * glint_bandwidth / nyquist
if highcut >= 1.0:
highcut = 0.99
b, a = butter(4, [lowcut, highcut], btype='band')
glint_noise = filtfilt(b, a, white_noise)
# 3. 归一化到典型的闪烁幅度
# 闪烁误差通常为波束宽度的 10-30%
glint_std = 0.15 # 波束宽度的 15%
glint_noise = glint_noise / np.std(glint_noise) * glint_std
# 4. 调制与目标方位相关
# 侧面视角闪烁更严重
aspect_factor = np.abs(np.sin(target_aspect))
glint_noise = glint_noise * (0.5 + 0.5 * aspect_factor)
return glint_noise
def apply_glint_to_angles(self, true_angles, fs):
"""将闪烁误差应用到真实角度上"""
num_samples = len(true_angles)
if isinstance(true_angles[0], (list, tuple, np.ndarray)):
# 二维角度
az_glint = self.generate_glint_error(num_samples, fs, true_angles[0])
el_glint = self.generate_glint_error(num_samples, fs, true_angles[1])
observed_az = true_angles[0] + az_glint
observed_el = true_angles[1] + el_glint
return observed_az, observed_el
else:
# 一维角度
glint_error = self.generate_glint_error(num_samples, fs, true_angles[0])
observed_angles = true_angles + glint_error
return observed_angles
class PulseCompressionChain:
"""完整的脉冲压缩处理链"""
def __init__(self, bandwidth, pulse_width, fs, window_type='hamming'):
self.bandwidth = bandwidth
self.pulse_width = pulse_width
self.fs = fs
self.window_type = window_type
# 生成波形和匹配滤波器
self.waveform = LFMWaveform(bandwidth, pulse_width, fs, 0)
self.tx_pulse, _, self.t = self.waveform.generate_pulse()
self.matched_filter = self._design_matched_filter()
# 性能参数
self.range_resolution = 3e8 / (2 * bandwidth)
self.unambiguous_range = 3e8 / (2 * fs) * len(self.tx_pulse)
def _design_matched_filter(self):
"""设计匹配滤波器（可加窗降低旁瓣）"""
mf = np.conj(self.tx_pulse[::-1])
# 加窗处理
if self.window_type == 'hamming':
window = np.hamming(len(mf))
elif self.window_type == 'chebyshev':
from scipy.signal import chebwin
window = chebwin(len(mf), at=50)
elif self.window_type == 'taylor':
# Taylor 窗实现
nbar = 4
sll = -35 # 旁瓣电平 dB
window = self._taylor_window(len(mf), nbar, sll)
else:
window = np.ones(len(mf))
mf_windowed = mf * window
# 重新归一化，避免信噪比损失
mf_windowed = mf_windowed / np.sqrt(np.sum(np.abs(window)**2))
return mf_windowed
def _taylor_window(self, N, nbar=4, sll=-30):
"""Taylor 加权窗函数"""
# Taylor 窗实现代码
# 此处为简化版本，实际需完整实现
from scipy.signal import taylor
return taylor(N, nbar, sll)
def process_echo(self, echo_signal, apply_cfar=True):
"""处理单个回波脉冲"""
# 1. 脉冲压缩
compressed = pulse_compression(echo_signal, self.matched_filter)
# 2. 包络检测
envelope = np.abs(compressed)
# 3. 对数压缩（用于显示）
log_envelope = 20 * np.log10(envelope + 1e-12)
# 4. CFAR 检测（可选）
if apply_cfar:
detections = self._cfar_detection(log_envelope)
else:
detections = None
# 5. 计算性能指标
metrics = self._compute_performance_metrics(compressed)
results = {
'compressed_signal': compressed,
'envelope': envelope,
'log_envelope': log_envelope,
'detections': detections,
'metrics': metrics,
'range_axis': np.arange(len(compressed)) * 3e8 / (2 * self.fs)
}
return results
def _cfar_detection(self, log_data, guard_cells=2, ref_cells=10, pfa=1e-6):
"""CFAR 目标检测"""
num_samples = len(log_data)
detections = np.zeros(num_samples, dtype=bool)
# 线性功率
linear_data = 10**(log_data / 10)
for i in range(num_samples):
# 保护单元
start_left = max(0, i - guard_cells - ref_cells)
end_left = max(0, i - guard_cells)
start_right = min(num_samples, i + guard_cells + 1)
end_right = min(num_samples, i + guard_cells + ref_cells + 1)
# 参考单元
ref_cells_left = linear_data[start_left:end_left]
ref_cells_right = linear_data[start_right:end_right]
ref_cells_all = np.concatenate([ref_cells_left, ref_cells_right])
if len(ref_cells_all) > 0:
# CA-CFAR：取平均
noise_estimate = np.mean(ref_cells_all)
# 计算阈值
from scipy import stats
threshold_factor = stats.chi2.ppf(1-pfa, 2*len(ref_cells_all)) / (2*len(ref_cells_all))
threshold = noise_estimate * threshold_factor
# 检测判决
detections[i] = linear_data[i] > threshold
return detections
def _compute_performance_metrics(self, compressed_signal):
"""计算脉冲压缩性能指标"""
envelope = np.abs(compressed_signal)
# 找到主瓣峰值
peak_idx = np.argmax(envelope)
peak_value = envelope[peak_idx]
# 计算主瓣宽度（3dB）
half_power = peak_value / np.sqrt(2)
left_idx = np.where(envelope[:peak_idx] >= half_power)[0]
right_idx = np.where(envelope[peak_idx:] >= half_power)[0]
if len(left_idx) > 0 and len(right_idx) > 0:
left_3db = left_idx[-1]
right_3db = peak_idx + right_idx[-1]
mainlobe_width_samples = right_3db - left_3db
mainlobe_width_meters = mainlobe_width_samples * 3e8 / (2 * self.fs)
else:
mainlobe_width_meters = None
# 计算峰值旁瓣比（PSLR）
sidelobe_region = np.concatenate([
envelope[:int(0.3*len(envelope))],
envelope[int(0.7*len(envelope)):]
])
if len(sidelobe_region) > 0:
max_sidelobe = np.max(sidelobe_region)
pslr_db = 20 * np.log10(peak_value / max_sidelobe)
else:
pslr_db = None
# 计算积分旁瓣比（ISLR）
mainlobe_energy = np.sum(envelope[max(0, peak_idx-10):min(len(envelope), peak_idx+10)]**2)
total_energy = np.sum(envelope**2)
islr_db = 10 * np.log10(mainlobe_energy / (total_energy - mainlobe_energy))
metrics = {
'peak_power_db': 20 * np.log10(peak_value),
'mainlobe_width_m': mainlobe_width_meters,
'pslr_db': pslr_db,
'islr_db': islr_db,
'compression_ratio': self.pulse_width * self.fs / mainlobe_width_samples if mainlobe_width_samples else None
}
return metrics
class ElectromagneticEnvironment:
"""复杂电磁环境建模"""
def __init__(self, config):
self.config = config
# 噪声设置
self.noise_temperature = config.get('noise_temperature', 290) # K
self.receiver_bandwidth = config.get('receiver_bandwidth', 10e6) # Hz
# 干扰设置
self.jammers = config.get('jammers', [])
self.clutter_model = config.get('clutter_model', 'constant_gamma')
def add_thermal_noise(self, signal, snr_db=None):
"""添加热噪声"""
if snr_db is not None:
# 根据指定 SNR 添加噪声
signal_power = np.mean(np.abs(signal)**2)
noise_power = signal_power / (10**(snr_db/10))
else:
# 根据系统噪声温度计算
k = 1.38e-23 # 玻尔兹曼常数
noise_power = k * self.noise_temperature * self.receiver_bandwidth
# 生成复高斯噪声
noise_std = np.sqrt(noise_power / 2) # 复噪声，实部虚部各一半功率
noise = noise_std * (np.random.randn(len(signal)) + 1j*np.random.randn(len(signal)))
noisy_signal = signal + noise
# 计算实际 SNR
actual_signal_power = np.mean(np.abs(signal)**2)
actual_noise_power = np.mean(np.abs(noise)**2)
actual_snr_db = 10 * np.log10(actual_signal_power / actual_noise_power)
return noisy_signal, actual_snr_db
def add_jamming(self, signal, time_axis):
"""添加干扰信号"""
jammed_signal = signal.copy()
jammer_info = []
for jammer in self.jammers:
jam_type = jammer.get('type', 'noise')
jam_power_db = jammer.get('power_db', 30)
jam_freq = jammer.get('frequency', 10e9)
jam_bw = jammer.get('bandwidth', 10e6)
jam_power_linear = 10**(jam_power_db/10)
if jam_type == 'noise':
# 噪声干扰
jammer_signal = np.sqrt(jam_power_linear) * (
np.random.randn(len(signal)) + 1j*np.random.randn(len(signal))
)
elif jam_type == 'cwi':
# 连续波干扰
jam_freq_offset = jammer.get('freq_offset', 0)
jammer_signal = np.sqrt(jam_power_linear) * np.exp(
1j * 2 * np.pi * jam_freq_offset * time_axis
)
elif jam_type == 'swept':
# 扫频干扰
sweep_rate = jammer.get('sweep_rate', 1e6) # Hz/s
jammer_signal = np.sqrt(jam_power_linear) * np.exp(
1j * 2 * np.pi * (0.5 * sweep_rate * time_axis**2)
)
# 应用干扰
jammed_signal += jammer_signal
jammer_info.append({
'type': jam_type,
'power_db': jam_power_db,
'jamming_to_signal_ratio': 10*np.log10(
np.mean(np.abs(jammer_signal)**2) / np.mean(np.abs(signal)**2)
)
})
return jammed_signal, jammer_info
def add_clutter(self, signal, range_axis, grazing_angle):
"""添加地/海杂波"""
clutter_signal = np.zeros_like(signal, dtype=complex)
if self.clutter_model == 'constant_gamma':
# 恒定γ模型
gamma_db = self.config.get('clutter_gamma_db', -20) # 典型值
gamma_linear = 10**(gamma_db/10)
# 杂波 RCS = γ * A_c，A_c 为照射面积
# 简化：假设与距离平方成正比
for i, R in enumerate(range_axis):
if R > 0:
# 简化的杂波功率模型
clutter_power = gamma_linear * R**2 * np.sin(grazing_angle)
clutter_amplitude = np.sqrt(clutter_power)
# 生成随机相位的杂波
clutter_phase = np.random.uniform(0, 2*np.pi)
clutter_signal[i] = clutter_amplitude * np.exp(1j * clutter_phase)
elif self.clutter_model == 'weibull':
# Weibull 分布杂波
scale = self.config.get('clutter_scale', 1.0)
shape = self.config.get('clutter_shape', 2.0)
# shape=2 时为瑞利分布
for i, R in enumerate(range_axis):
if R > 0:
# Weibull 分布幅度
clutter_amplitude = scale * (-np.log(np.random.rand()))**(1/shape)
clutter_phase = np.random.uniform(0, 2*np.pi)
clutter_signal[i] = clutter_amplitude * np.exp(1j * clutter_phase)
# 添加多普勒扩展（对于运动平台）
if self.config.get('clutter_doppler_spread', 0) > 0:
doppler_spread = self.config.get('clutter_doppler_spread')
clutter_signal = self._apply_doppler_spread(clutter_signal, doppler_spread)
clutter_signal *= np.sqrt(len(signal)) # 保持能量
clutter_power_db = 10 * np.log10(np.mean(np.abs(clutter_signal)**2))
return signal + clutter_signal, clutter_power_db
def _apply_doppler_spread(self, signal, doppler_spread):
"""应用多普勒扩展"""
from scipy.signal import filtfilt, butter
# 设计一个多普勒滤波器模拟杂波谱
b, a = butter(4, doppler_spread, 'low')
I = filtfilt(b, a, np.random.randn(len(signal)))
Q = filtfilt(b, a, np.random.randn(len(signal)))
doppler_modulated = signal * (I + 1j*Q)
return doppler_modulated
def get_environment_summary(self):
"""获取环境配置摘要"""
summary = "电磁环境配置：\n"
summary += f"噪声温度：{self.noise_temperature} K\n"
summary += f"接收机带宽：{self.receiver_bandwidth/1e6:.1f} MHz\n"
summary += f"杂波模型：{self.clutter_model}\n"
if self.jammers:
summary += "干扰设置：\n"
for i, jammer in enumerate(self.jammers):
summary += f" 干扰源{i+1}: 类型={jammer.get('type')}, "
summary += f"功率={jammer.get('power_db')}dB, "
summary += f"频率={jammer.get('frequency', 0)/1e9:.2f}GHz\n"
return summary
def run_complete_simulation_example():
"""运行完整的仿真示例"""
# 1. 仿真参数配置
config = {
'simulation': {
'duration': 0.1, # 仿真时长 100ms
'fs': 50e6, # 采样率 50MHz
'prf': 1000, # 脉冲重复频率 1kHz
'num_pulses': 100, # 脉冲数
},
'radar': {
'fc': 10e9, # 载频 10GHz
'bandwidth': 10e6, # 带宽 10MHz
'pulse_width': 10e-6, # 脉宽 10μs
'beamwidth_az': np.radians(3.0), # 3 度波束宽度
'beamwidth_el': np.radians(3.0),
},
'target': {
'initial_range': 20000, # 初始距离 20km
'radial_velocity': -300, # 径向速度 -300m/s
'rcs_mean': 5.0, # 平均 RCS 5m²
'rcs_type': 'Swerling1', # 起伏类型
'glint_enabled': True, # 启用角闪烁
},
'channel_errors': {
'gain_imbalance_db': 0.8, # 0.8dB 幅度不平衡
'phase_imbalance_deg': 3.5, # 3.5 度相位不平衡
'delay_imbalance_sec': 15e-12, # 15ps 延迟不平衡
},
'environment': {
'noise_temperature': 500, # 噪声温度 500K
'jammers': [
{ 'type': 'noise', 'power_db': 40, 'bandwidth': 20e6, }
],
'clutter_model': 'constant_gamma',
'clutter_gamma_db': -25,
}
}
# 2. 创建仿真器
simulator = MonopulseEchoSimulator(config)
# 3. 运行仿真
print("开始高保真单脉冲回波仿真...")
print("="*50)
print("仿真参数：")
print(f" 载频：{config['radar']['fc']/1e9:.2f} GHz")
print(f" 带宽：{config['radar']['bandwidth']/1e6:.1f} MHz")
print(f" 脉宽：{config['radar']['pulse_width']*1e6:.1f} μs")
print(f" 脉冲数：{config['simulation']['num_pulses']}")
print(f" 目标距离：{config['target']['initial_range']/1000:.1f} km")
print(f" 目标速度：{config['target']['radial_velocity']} m/s")
print("="*50)
results = simulator.run_simulation()
# 4. 显示仿真结果摘要
print("\n仿真结果摘要：")
print(f" 和通道数据形状：{results['sum_channel'].shape}")
print(f" 方位差通道数据形状：{results['diff_az'].shape}")
print(f" 俯仰差通道数据形状：{results['diff_el'].shape}")
print(f" 时间轴长度：{len(results['time_axis'])}")
return results, config
def visualize_simulation_results(results, config):
"""可视化仿真结果"""
# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
# 创建图形
fig = plt.figure(figsize=(20, 16))
# 1. 时域波形
ax1 = plt.subplot(3, 3, 1)
pulse_idx = 0
time_samples = np.arange(1000) / config['simulation']['fs'] * 1e6 # μs
ax1.plot(time_samples, np.real(results['sum_channel'][pulse_idx, :1000]), label='实部', linewidth=1)
ax1.plot(time_samples, np.imag(results['sum_channel'][pulse_idx, :1000]), label='虚部', linewidth=1, alpha=0.7)
ax1.set_xlabel('时间 (μs)')
ax1.set_ylabel('幅度')
ax1.set_title('单脉冲回波时域波形 (第 1 个脉冲)')
ax1.legend()
ax1.grid(True, alpha=0.3)
# 2. 频域分析
ax2 = plt.subplot(3, 3, 2)
fft_sum = np.fft.fft(results['sum_channel'][pulse_idx])
fft_diff_az = np.fft.fft(results['diff_az'][pulse_idx])
freq_axis = np.fft.fftfreq(len(fft_sum), 1/config['simulation']['fs']) / 1e6
ax2.plot(freq_axis[:len(freq_axis)//2], 20*np.log10(np.abs(fft_sum[:len(fft_sum)//2])), label='和通道', linewidth=1.5)
ax2.plot(freq_axis[:len(freq_axis)//2], 20*np.log10(np.abs(fft_diff_az[:len(fft_diff_az)//2])), label='差通道', linewidth=1.5, alpha=0.7)
ax2.set_xlabel('频率 (MHz)')
ax2.set_ylabel('功率谱密度 (dB)')
ax2.set_title('和差通道频谱对比')
ax2.legend()
ax2.grid(True, alpha=0.3)
ax2.set_xlim([-config['radar']['bandwidth']/2e6, config['radar']['bandwidth']/2e6])
# 3. 脉冲压缩结果
ax3 = plt.subplot(3, 3, 3)
from scipy.signal import hilbert
# 对第一个脉冲进行脉冲压缩
waveform = LFMWaveform(
config['radar']['bandwidth'],
config['radar']['pulse_width'],
config['simulation']['fs'],
0
)
matched_filter = waveform.matched_filter()
compressed = pulse_compression(results['sum_channel'][pulse_idx], matched_filter)
envelope = np.abs(compressed)
range_axis = np.arange(len(envelope)) * 3e8 / (2 * config['simulation']['fs'])
ax3.plot(range_axis, 20*np.log10(envelope), linewidth=1.5)
ax3.set_xlabel('距离 (m)')
ax3.set_ylabel('幅度 (dB)')
ax3.set_title('脉冲压缩结果 - 距离像')
ax3.grid(True, alpha=0.3)
ax3.set_xlim([0, 10000])
# 4. 距离 - 多普勒图
ax4 = plt.subplot(3, 3, 4)
# 对多个脉冲进行脉冲压缩
num_pulses = min(64, results['sum_channel'].shape[0])
range_profiles = []
for i in range(num_pulses):
compressed = pulse_compression(results['sum_channel'][i], matched_filter)
range_profiles.append(np.abs(compressed[:500]))
range_profiles = np.array(range_profiles)
# 距离维 FFT
range_fft = np.fft.fft(range_profiles, axis=0)
doppler_axis = np.fft.fftshift(np.fft.fftfreq(num_pulses, 1/config['simulation']['prf']))
im = ax4.imshow(20*np.log10(np.abs(range_fft)), aspect='auto', extent=[0, 500 * 3e8/(2*config['simulation']['fs']), doppler_axis[0], doppler_axis[-1]], cmap='jet', origin='lower')
ax4.set_xlabel('距离 (m)')
ax4.set_ylabel('多普勒频率 (Hz)')
ax4.set_title('距离 - 多普勒图')
plt.colorbar(im, ax=ax4, label='功率 (dB)')
# 5. 和差通道相关性分析
ax5 = plt.subplot(3, 3, 5)
# 取中间距离单元
range_bin = np.argmax(np.mean(np.abs(results['sum_channel']), axis=0))
sum_samples = results['sum_channel'][:, range_bin]
diff_samples = results['diff_az'][:, range_bin]
ax5.scatter(np.real(sum_samples), np.real(diff_samples), s=10, alpha=0.6, label='实部')
ax5.scatter(np.imag(sum_samples), np.imag(diff_samples), s=10, alpha=0.6, label='虚部', color='red')
# 拟合直线
from scipy.stats import linregress
slope_real, intercept_real, r_value_real, _, _ = linregress(
np.real(sum_samples), np.real(diff_samples)
)
slope_imag, intercept_imag, r_value_imag, _, _ = linregress(
np.imag(sum_samples), np.imag(diff_samples)
)
x_fit = np.array([np.min(np.real(sum_samples)), np.max(np.real(sum_samples))])
ax5.plot(x_fit, slope_real*x_fit + intercept_real, 'b-', label=f'实部：R={r_value_real:.3f}')
ax5.plot(x_fit, slope_imag*x_fit + intercept_imag, 'r-', label=f'虚部：R={r_value_imag:.3f}')
ax5.set_xlabel('和通道 (I/Q)')
ax5.set_ylabel('差通道 (I/Q)')
ax5.set_title('和差通道相关性分析')
ax5.legend()
ax5.grid(True, alpha=0.3)
ax5.axis('equal')
# 6. 相位关系分析
ax6 = plt.subplot(3, 3, 6)
phase_diff = np.angle(sum_samples * np.conj(diff_samples))
phase_diff_deg = np.degrees(phase_diff)
ax6.hist(phase_diff_deg, bins=30, alpha=0.7, edgecolor='black')
ax6.axvline(np.mean(phase_diff_deg), color='red', linestyle='--', linewidth=2, label=f'均值：{np.mean(phase_diff_deg):.2f}°')
ax6.axvline(np.std(phase_diff_deg), color='green', linestyle=':', linewidth=2, label=f'标准差：{np.std(phase_diff_deg):.2f}°')
ax6.set_xlabel('和差通道相位差 (度)')
ax6.set_ylabel('频数')
ax6.set_title('和差通道相位差分布')
ax6.legend()
ax6.grid(True, alpha=0.3)
# 7. 通道不平衡分析
ax7 = plt.subplot(3, 3, 7)
# 计算通道幅度比
amplitude_ratio = np.abs(diff_samples) / (np.abs(sum_samples) + 1e-12)
amplitude_ratio_db = 20 * np.log10(amplitude_ratio)
ax7.hist(amplitude_ratio_db, bins=30, alpha=0.7, edgecolor='black')
ax7.axvline(np.mean(amplitude_ratio_db), color='red', linestyle='--', linewidth=2, label=f'均值：{np.mean(amplitude_ratio_db):.2f} dB')
ax7.set_xlabel('差和通道幅度比 (dB)')
ax7.set_ylabel('频数')
ax7.set_title('通道幅度不平衡分析')
ax7.legend()
ax7.grid(True, alpha=0.3)
# 8. 时频分析
ax8 = plt.subplot(3, 3, 8)
from scipy.signal import spectrogram
f, t, Sxx = spectrogram(results['sum_channel'][:, range_bin], fs=config['simulation']['prf'], nperseg=32, noverlap=16)
im = ax8.pcolormesh(t, f, 10*np.log10(Sxx), shading='gouraud', cmap='viridis')
ax8.set_xlabel('时间 (s)')
ax8.set_ylabel('频率 (Hz)')
ax8.set_title('和通道时频分析')
plt.colorbar(im, ax=ax8, label='功率 (dB)')
# 9. 误差统计
ax9 = plt.subplot(3, 3, 9)
# 计算各种误差指标
error_metrics = {
'幅度不平衡 (dB)': np.std(amplitude_ratio_db),
'相位不平衡 (度)': np.std(phase_diff_deg),
'I/Q 不平衡 (dB)': 0, # 需要从数据计算
'信噪比 (dB)': 15, # 需要从数据计算
'信杂比 (dB)': 20, # 需要从数据计算
}
metrics_names = list(error_metrics.keys())
metrics_values = list(error_metrics.values())
bars = ax9.barh(metrics_names, metrics_values, color='steelblue', alpha=0.7)
ax9.set_xlabel('误差值')
ax9.set_title('系统误差统计')
ax9.grid(True, alpha=0.3, axis='x')
# 在柱状图上添加数值
for bar, value in zip(bars, metrics_values):
ax9.text(value, bar.get_y() + bar.get_height()/2, f' {value:.2f}', va='center', ha='left')
plt.tight_layout()
plt.suptitle('单脉冲回波仿真结果综合分析', fontsize=16, y=1.02)
plt.show()
return fig
def evaluate_system_performance(results, config):
"""系统性能综合评估"""
print("="*60)
print("单脉冲回波生成系统性能评估")
print("="*60)
# 1. 数据质量检查
print("\n1. 数据质量检查:")
print("-"*40)
# 检查数据维度
num_pulses, samples_per_pulse = results['sum_channel'].shape
print(f" 数据维度：{num_pulses}个脉冲 × {samples_per_pulse}个采样点")
# 检查 NaN 或 Inf 值
has_nan = np.any(np.isnan(results['sum_channel']))
has_inf = np.any(np.isinf(results['sum_channel']))
print(f" 数据有效性：NaN={has_nan}, Inf={has_inf}")
# 检查动态范围
dynamic_range_db = 20 * np.log10(np.max(np.abs(results['sum_channel'])) / np.std(results['sum_channel'][np.abs(results['sum_channel']) > 0]))
print(f" 动态范围：{dynamic_range_db:.2f} dB")
# 2. 通道一致性评估
print("\n2. 通道一致性评估:")
print("-"*40)
# 选取中间距离单元进行评估
mid_range_bin = samples_per_pulse // 2
sum_channel_mid = results['sum_channel'][:, mid_range_bin]
diff_az_mid = results['diff_az'][:, mid_range_bin]
# 幅度一致性
amplitude_ratio = np.abs(diff_az_mid) / (np.abs(sum_channel_mid) + 1e-12)
amp_mean_db = 20 * np.log10(np.mean(amplitude_ratio))
amp_std_db = 20 * np.log10(np.std(amplitude_ratio))
print(f" 幅度一致性:")
print(f" - 均值：{amp_mean_db:.2f} dB")
print(f" - 标准差：{amp_std_db:.2f} dB")
# 相位一致性
phase_diff = np.angle(sum_channel_mid * np.conj(diff_az_mid))
phase_mean_deg = np.degrees(np.mean(phase_diff))
phase_std_deg = np.degrees(np.std(phase_diff))
print(f" 相位一致性:")
print(f" - 均值：{phase_mean_deg:.2f} °")
print(f" - 标准差：{phase_std_deg:.2f} °")
# 3. 信噪比估计
print("\n3. 信噪比估计:")
print("-"*40)
# 使用脉冲压缩后的峰值进行估计
from scipy.signal import find_peaks
# 对每个脉冲进行脉冲压缩
waveform = LFMWaveform(
config['radar']['bandwidth'],
config['radar']['pulse_width'],
config['simulation']['fs'],
0
)
matched_filter = waveform.matched_filter()
snr_list = []
for i in range(min(10, num_pulses)):
# 检查前 10 个脉冲
compressed = pulse_compression(results['sum_channel'][i], matched_filter)
envelope = np.abs(compressed)
# 找到主瓣峰值
peaks, properties = find_peaks(envelope, height=np.max(envelope)*0.5, distance=10)
if len(peaks) > 0:
main_peak = peaks[0]
main_peak_power = envelope[main_peak]
# 计算噪声（主瓣外区域）
noise_region = np.concatenate([envelope[:main_peak-20], envelope[main_peak+20:]])
noise_power = np.mean(noise_region**2)
snr = 10 * np.log10(main_peak_power**2 / noise_power)
snr_list.append(snr)
if snr_list:
avg_snr = np.mean(snr_list)
std_snr = np.std(snr_list)
print(f" 脉冲压缩后 SNR: {avg_snr:.2f} ± {std_snr:.2f} dB")
# 4. 脉冲压缩性能
print("\n4. 脉冲压缩性能:")
print("-"*40)
# 对第一个脉冲进行详细分析
compressed = pulse_compression(results['sum_channel'][0], matched_filter)
envelope = np.abs(compressed)
# 找到主瓣
peak_idx = np.argmax(envelope)
peak_value = envelope[peak_idx]
# 计算 3dB 主瓣宽度
half_power = peak_value / np.sqrt(2)
left_idx = np.where(envelope[:peak_idx] >= half_power)[0]
right_idx = np.where(envelope[peak_idx:] >= half_power)[0]
if len(left_idx) > 0 and len(right_idx) > 0:
left_3db = left_idx[-1]
right_3db = peak_idx + right_idx[-1]
mainlobe_width_samples = right_3db - left_3db
mainlobe_width_meters = mainlobe_width_samples * 3e8 / (2 * config['simulation']['fs'])
print(f" 距离分辨率：{mainlobe_width_meters:.2f} m")
print(f" 理论分辨率：{3e8/(2*config['radar']['bandwidth']):.2f} m")
# 计算 PSLR
sidelobe_region = np.concatenate([envelope[:peak_idx-20], envelope[peak_idx+20:]])
if len(sidelobe_region) > 0:
max_sidelobe = np.max(sidelobe_region)
pslr_db = 20 * np.log10(peak_value / max_sidelobe)
print(f" 峰值旁瓣比：{pslr_db:.2f} dB")
# 5. 和差通道关系验证
print("\n5. 和差通道关系验证:")
print("-"*40)
# 验证 S 曲线特性
test_angles = np.radians(np.linspace(-5, 5, 21))
s_curve_values = []
for angle in test_angles:
# 简化验证：使用理想方向图
sum_response = np.exp(-(angle**2) / (2*(np.radians(3)/2.35)**2)) # 高斯波束
diff_response = angle * np.exp(-(angle**2) / (2*(np.radians(3)/2.35)**2)) # 奇对称
s_value = diff_response * sum_response # 忽略共轭，因为都是实数
s_curve_values.append(s_value)
# 计算线性度
from scipy.stats import linregress
linear_region = (np.abs(test_angles) < np.radians(2))
slope, intercept, r_value, p_value, std_err = linregress(
test_angles[linear_region], np.array(s_curve_values)[linear_region]
)
print(f" S 曲线线性区域斜率：{slope:.4f}")
print(f" 线性相关系数 R: {r_value:.4f}")
print(f" 线性度误差：{std_err:.6f}")
# 6. 误差源分析
print("\n6. 误差源贡献分析:")
print("-"*40)
error_sources = {
'通道幅度不平衡': config['channel_errors']['gain_imbalance_db'],
'通道相位不平衡': config['channel_errors']['phase_imbalance_deg'],
'通道延迟不平衡': config['channel_errors']['delay_imbalance_sec'] * 1e12,
'I/Q 不平衡': 0.5, # 假设值
'热噪声': config['environment']['noise_temperature'],
}
for source, value in error_sources.items():
unit = 'dB' if 'dB' in source else 'ps' if '延迟' in source else 'K' if '热噪声' in source else '度'
print(f" {source}: {value} {unit}")
# 7. 蒙特卡洛统计分析
print("\n7. 蒙特卡洛统计分析:")
print("-"*40)
# 进行多次仿真，统计性能变化
num_monte_carlo = 10
snr_stats = []
phase_std_stats = []
for mc_idx in range(num_monte_carlo):
# 这里简化，实际应该重新运行仿真
# 使用现有数据加入随机噪声模拟
noise_power = 0.1
noisy_sum = results['sum_channel'] + np.sqrt(noise_power/2) * (
np.random.randn(*results['sum_channel'].shape) + 1j*np.random.randn(*results['sum_channel'].shape)
)
# 计算相位一致性
sum_sample = noisy_sum[:, mid_range_bin]
diff_sample = results['diff_az'][:, mid_range_bin]
phase_diff = np.angle(sum_sample * np.conj(diff_sample))
phase_std = np.std(phase_diff)
phase_std_stats.append(np.degrees(phase_std))
phase_std_mean = np.mean(phase_std_stats)
phase_std_std = np.std(phase_std_stats)
print(f" 相位一致性标准差统计:")
print(f" - 均值：{phase_std_mean:.2f} °")
print(f" - 标准差：{phase_std_std:.2f} °")
print(f" - 变异系数：{phase_std_std/phase_std_mean*100:.1f}%")
print("\n" + "="*60)
print("性能评估完成")
print("="*60)
return {
'dynamic_range_db': dynamic_range_db,
'amplitude_consistency_mean_db': amp_mean_db,
'phase_consistency_std_deg': phase_std_deg,
'snr_db': avg_snr if 'avg_snr' in locals() else None,
'range_resolution_m': mainlobe_width_meters if 'mainlobe_width_meters' in locals() else None,
's_curve_linearity_r': r_value,
}
class MultiTargetSimulator:
"""多目标场景仿真"""
def __init__(self, base_simulator):
self.simulator = base_simulator
self.targets = []
def add_target(self, range_m, radial_velocity, rcs_db, azimuth_deg, elevation_deg):
"""添加目标"""
target = {
'range': range_m,
'velocity': radial_velocity,
'rcs_db': rcs_db,
'azimuth': np.radians(azimuth_deg),
'elevation': np.radians(elevation_deg),
'signal': None
}
self.targets.append(target)
def simulate_multi_target(self):
"""仿真多目标场景"""
num_targets = len(self.targets)
print(f"开始多目标仿真，目标数量：{num_targets}")
# 获取仿真参数
config = self.simulator.config
num_pulses = config['simulation']['num_pulses']
samples_per_pulse = int(config['simulation']['fs'] * config['radar']['pulse_width'])
# 初始化接收信号
sum_channel = np.zeros((num_pulses, samples_per_pulse), dtype=complex)
diff_az_channel = np.zeros_like(sum_channel)
diff_el_channel = np.zeros_like(sum_channel)
# 对每个目标单独仿真
for i, target in enumerate(self.targets):
print(f" 生成目标{i+1}回波...")
# 临时修改仿真器目标参数
self.simulator.target.range = target['range']
self.simulator.target.velocity = target['velocity']
self.simulator.target.rcs_db = target['rcs_db']
# 计算该目标的回波
target_results = self.simulator._simulate_single_target(
target['azimuth'], target['elevation']
)
# 累加到总信号中
sum_channel += target_results['sum']
diff_az_channel += target_results['diff_az']
diff_el_channel += target_results['diff_el']
target['signal'] = target_results['sum']
# 添加噪声
env = ElectromagneticEnvironment(config['environment'])
time_axis = np.arange(samples_per_pulse) / config['simulation']['fs']
for pulse_idx in range(num_pulses):
sum_channel[pulse_idx], _ = env.add_thermal_noise(
sum_channel[pulse_idx], snr_db=20
)
diff_az_channel[pulse_idx], _ = env.add_thermal_noise(
diff_az_channel[pulse_idx], snr_db=20
)
diff_el_channel[pulse_idx], _ = env.add_thermal_noise(
diff_el_channel[pulse_idx], snr_db=20
)
results = {
'sum_channel': sum_channel,
'diff_az': diff_az_channel,
'diff_el': diff_el_channel,
'time_axis': np.arange(num_pulses) / config['simulation']['prf'],
'targets': self.targets,
'config': config
}
return results
def _simulate_single_target(self, azimuth, elevation):
"""仿真单个目标（简化的内部方法）"""
# 这里调用基础仿真器的内部方法
# 注意：实际实现中需要根据基础仿真器的结构进行调整
pass
class WidebandSimulator:
"""宽带信号仿真器，支持高分辨率成像"""
def __init__(self, fc, bandwidth, pulse_width, fs):
self.fc = fc # 中心频率
self.B = bandwidth # 带宽
self.T = pulse_width # 脉冲宽度
self.fs = fs # 采样率
# 子带划分
self.num_subbands = 8
self.subband_width = self.B / self.num_subbands
def generate_stepped_freq_signal(self):
"""生成步进频信号"""
# 步进频信号：多个连续脉冲，每个脉冲频率递增
num_pulses = 64
freq_steps = np.linspace(self.fc - self.B/2, self.fc + self.B/2, num_pulses)
signals = []
for i, freq in enumerate(freq_steps):
t = np.arange(int(self.T * self.fs)) / self.fs
pulse = np.exp(1j * 2 * np.pi * freq * t)
signals.append(pulse)
return np.array(signals), freq_steps
def generate_ofdm_signal(self):
"""生成 OFDM 雷达信号"""
num_subcarriers = 64
num_symbols = 10 # OFDM 参数
subcarrier_spacing = self.B / num_subcarriers
symbol_duration = 1 / subcarrier_spacing
# 生成随机 QPSK 符号
symbols = (
np.random.randint(0, 2, (num_symbols, num_subcarriers)) * 2 - 1
) + 1j * (
np.random.randint(0, 2, (num_symbols, num_subcarriers)) * 2 - 1
)
# IFFT 生成时域信号
ofdm_symbols = np.fft.ifft(symbols, axis=1) * np.sqrt(num_subcarriers)
# 添加循环前缀
cp_length = num_subcarriers // 4
ofdm_symbols_cp = np.hstack([ofdm_symbols[:, -cp_length:], ofdm_symbols])
# 串并转换
ofdm_signal = ofdm_symbols_cp.flatten()
return ofdm_signal, symbols
def range_profile_synthesis(self, stepped_freq_signals):
"""距离像合成（步进频信号处理）"""
num_pulses, samples_per_pulse = stepped_freq_signals.shape
# 对每个脉冲进行 FFT（距离像）
range_profiles = np.zeros((num_pulses, samples_per_pulse), dtype=complex)
for i in range(num_pulses):
range_profiles[i] = np.fft.fft(stepped_freq_signals[i])
# 合成高分辨率距离像
# 这里可以使用逆 FFT 或其他合成算法
synthesized_profile = np.sum(range_profiles, axis=0)
return synthesized_profile, range_profiles
def sar_imaging_simulation(self, platform_velocity, synthetic_aperture_length):
"""合成孔径雷达（SAR）成像仿真"""
# 简化版 SAR 仿真
num_pulses = 256
range_bins = 512
# 生成点目标场景
targets = [
{'range': 5000, 'azimuth': 0, 'rcs': 1.0},
{'range': 5050, 'azimuth': 0.1, 'rcs': 0.8},
{'range': 5100, 'azimuth': -0.1, 'rcs': 0.6},
]
# 初始化数据矩阵
sar_data = np.zeros((num_pulses, range_bins), dtype=complex)
# 平台位置序列
platform_positions = np.linspace(-synthetic_aperture_length/2, synthetic_aperture_length/2, num_pulses)
# 对每个目标生成回波
for target in targets:
for i, platform_pos in enumerate(platform_positions):
# 计算斜距历程
R = np.sqrt(target['range']**2 + (platform_pos - target['azimuth'])**2)
# 计算延迟
delay_samples = int(2 * R / 3e8 * self.fs)
if delay_samples < range_bins:
# 生成 LFM 回波（简化）
t = np.arange(int(self.T * self.fs)) / self.fs
chirp = np.exp(1j * np.pi * self.B / self.T * t**2)
# 应用延迟和相位
sar_data[i, delay_samples:delay_samples+len(chirp)] += \
chirp * target['rcs'] * np.exp(-1j * 4 * np.pi * R / 0.03)
# 添加噪声
noise = 0.1 * (np.random.randn(*sar_data.shape) + 1j*np.random.randn(*sar_data.shape))
sar_data += noise
return sar_data, platform_positions, targets
def isar_imaging_simulation(self, target_rotation_rate, imaging_time):
"""逆合成孔径雷达（ISAR）成像仿真"""
# 旋转目标 ISAR 成像
num_pulses = 128
range_bins = 256
doppler_bins = 128
# 目标散射点模型
scatterers = [
{'x': 0, 'y': 0, 'z': 0, 'rcs': 1.0},
{'x': 2, 'y': 0, 'z': 0, 'rcs': 0.8},
{'x': 0, 'y': 2, 'z': 0, 'rcs': 0.6},
{'x': -2, 'y': 0, 'z': 0, 'rcs': 0.7},
{'x': 0, 'y': -2, 'z': 0, 'rcs': 0.5},
]
# 时间序列
time_axis = np.linspace(0, imaging_time, num_pulses)
# 初始化距离 - 慢时间矩阵
range_time_matrix = np.zeros((num_pulses, range_bins), dtype=complex)
for i, t in enumerate(time_axis):
# 目标旋转角度
rotation_angle = target_rotation_rate * t
for scatterer in scatterers:
# 旋转后的坐标
x_rot = scatterer['x'] * np.cos(rotation_angle) - scatterer['y'] * np.sin(rotation_angle)
y_rot = scatterer['x'] * np.sin(rotation_angle) + scatterer['y'] * np.cos(rotation_angle)
# 计算距离（简化：假设视线沿 x 轴）
R = np.sqrt(x_rot**2 + y_rot**2)
# 计算多普勒频移
# 速度 = 旋转角速度 × 半径（在视线方向的投影）
v = target_rotation_rate * y_rot
doppler_freq = 2 * v / 0.03 # 波长假设为 0.03m
# 生成距离单元响应
range_bin = int(R / (3e8/(2*self.fs)))
if 0 <= range_bin < range_bins:
# 添加散射点贡献
range_time_matrix[i, range_bin] += \
scatterer['rcs'] * np.exp(1j * 2 * np.pi * doppler_freq * t)
# 距离维脉冲压缩（如果使用 LFM）
# 这里简化处理
# 多普勒处理（生成 ISAR 像）
isar_image = np.zeros((range_bins, doppler_bins), dtype=complex)
for r in range(range_bins):
# 对每个距离单元进行多普勒分析
doppler_profile = np.fft.fft(range_time_matrix[:, r], doppler_bins)
isar_image[r, :] = doppler_profile
return isar_image, range_time_matrix, scatterers
class AdaptiveBeamforming:
"""自适应数字波束形成（ADBF）"""
def __init__(self, num_elements, element_spacing, wavelength):
self.N = num_elements # 阵元数
self.d = element_spacing # 阵元间距
self.lambda_ = wavelength # 波长
# 阵列流形矩阵
self.steering_vectors = {}
def compute_steering_vector(self, theta):
"""计算导向矢量"""
# theta: 波达方向（弧度）
n = np.arange(self.N)
phase_shift = 2 * np.pi * self.d / self.lambda_ * np.sin(theta)
a = np.exp(1j * n * phase_shift)
return a.reshape(-1, 1) # 列向量
def conventional_beamforming(self, received_data, look_direction):
"""常规波束形成"""
# received_data: N×M 矩阵，N 阵元，M 快拍
# look_direction: 波束指向（弧度）
w = self.compute_steering_vector(look_direction)
w = w / np.sqrt(self.N) # 归一化
# 波束形成输出
beam_output = w.conj().T @ received_data
return beam_output, w
def mvdr_beamformer(self, received_data, look_direction, interference_angles=None):
"""最小方差无失真响应（MVDR）波束形成器"""
# 计算样本协方差矩阵
R_xx = received_data @ received_data.conj().T / received_data.shape[1]
# 添加对角加载改善数值稳定性
R_xx += 1e-6 * np.eye(self.N)
# 计算导向矢量
a = self.compute_steering_vector(look_direction)
# MVDR 权值
R_inv = np.linalg.inv(R_xx)
w_mvdr = (R_inv @ a) / (a.conj().T @ R_inv @ a)
# 波束形成输出
beam_output = w_mvdr.conj().T @ received_data
return beam_output, w_mvdr, R_xx
def lcmv_beamformer(self, received_data, constraint_matrix, response_vector):
"""线性约束最小方差（LCMV）波束形成器"""
# constraint_matrix: 约束矩阵（N×C）
# response_vector: 期望响应（C×1）
# 样本协方差矩阵
R_xx = received_data @ received_data.conj().T / received_data.shape[1]
R_inv = np.linalg.inv(R_xx + 1e-6*np.eye(self.N))
# LCMV 权值
C = constraint_matrix
f = response_vector
w_lcmv = R_inv @ C @ np.linalg.inv(C.conj().T @ R_inv @ C) @ f
# 波束形成输出
beam_output = w_lcmv.conj().T @ received_data
return beam_output, w_lcmv
def music_doa(self, received_data, num_sources, search_grid=np.linspace(-np.pi/2, np.pi/2, 361)):
"""MUSIC 算法进行波达方向估计"""
# num_sources: 信源数
# 计算样本协方差矩阵
R_xx = received_data @ received_data.conj().T / received_data.shape[1]
# 特征值分解
eigenvalues, eigenvectors = np.linalg.eig(R_xx)
# 按特征值排序
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
# 信号子空间和噪声子空间
signal_subspace = eigenvectors[:, :num_sources]
noise_subspace = eigenvectors[:, num_sources:]
# MUSIC 谱
music_spectrum = np.zeros(len(search_grid))
for i, theta in enumerate(search_grid):
a = self.compute_steering_vector(theta)
music_spectrum[i] = 1 / (a.conj().T @ noise_subspace @ noise_subspace.conj().T @ a + 1e-12)
music_spectrum = np.abs(music_spectrum).flatten()
# 寻找峰值
from scipy.signal import find_peaks
peaks, properties = find_peaks(music_spectrum, height=np.max(music_spectrum)*0.1, distance=10)
doa_estimates = search_grid[peaks]
return music_spectrum, doa_estimates, search_grid, eigenvalues
def plot_beam_pattern(self, weights, title="波束方向图"):
"""绘制波束方向图"""
theta_grid = np.linspace(-np.pi/2, np.pi/2, 361)
beam_pattern = np.zeros(len(theta_grid), dtype=complex)
for i, theta in enumerate(theta_grid):
a = self.compute_steering_vector(theta)
beam_pattern[i] = weights.conj().T @ a
beam_pattern_db = 20 * np.log10(np.abs(beam_pattern) + 1e-12)
plt.figure(figsize=(10, 6))
plt.plot(np.degrees(theta_grid), beam_pattern_db)
plt.xlabel('角度 (度)')
plt.ylabel('增益 (dB)')
plt.title(title)
plt.grid(True)
plt.ylim([-50, 5])
return beam_pattern_db, theta_grid
class PolarimetricRadarSimulator:
"""极化雷达仿真器"""
def __init__(self):
# 极化散射矩阵（Sinclair 矩阵）
self.S_matrix = np.eye(2, dtype=complex) # 初始为单位矩阵（各向同性散射）
# 发射和接收极化基
self.transmit_polarization = np.array([1, 0]) # 水平极化
self.receive_polarization = np.array([1, 0]) # 水平极化
def set_target_scattering_matrix(self, shh, shv, svh, svv):
"""设置目标散射矩阵"""
# shh: 水平发射 - 水平接收
# shv: 水平发射 - 垂直接收
# svh: 垂直发射 - 水平接收
# svv: 垂直发射 - 垂直接收
self.S_matrix = np.array([[shh, shv], [svh, svv]], dtype=complex)
def set_polarization_basis(self, transmit, receive):
"""设置发射和接收极化"""
self.transmit_polarization = transmit / np.linalg.norm(transmit)
self.receive_polarization = receive / np.linalg.norm(receive)
def compute_scattering_coefficient(self):
"""计算散射系数"""
# 接收电压 = receive^H * S * transmit
scattering_coeff = self.receive_polarization.conj().T @ \
self.S_matrix @ \
self.transmit_polarization
return scattering_coeff
def simulate_full_polarimetric_data(self):
"""仿真全极化测量数据"""
# 四种极化组合：HH, HV, VH, VV
polarizations = {
'HH': (np.array([1, 0]), np.array([1, 0])),
'HV': (np.array([1, 0]), np.array([0, 1])),
'VH': (np.array([0, 1]), np.array([1, 0])),
'VV': (np.array([0, 1]), np.array([0, 1]))
}
polarimetric_data = {}
for pol_name, (tx_pol, rx_pol) in polarizations.items():
self.set_polarization_basis(tx_pol, rx_pol)
scattering_coeff = self.compute_scattering_coefficient()
polarimetric_data[pol_name] = scattering_coeff
return polarimetric_data
def compute_polarimetric_descriptors(self, polarimetric_data):
"""计算极化描述子"""
# 获取散射矩阵元素
Shh = polarimetric_data['HH']
Shv = polarimetric_data['HV']
Svh = polarimetric_data['VH']
Svv = polarimetric_data['VV']
# 构建散射矩阵
S = np.array([[Shh, Shv], [Svh, Svv]])
# 计算极化熵
# 首先计算相干矩阵 T
k = np.array([Shh, Shv + Svh, Svv]) / np.sqrt(2) # Pauli 基矢量
T = np.outer(k, k.conj())
# 特征值分解
eigenvalues, _ = np.linalg.eig(T)
eigenvalues = np.real(eigenvalues)
eigenvalues = eigenvalues / np.sum(eigenvalues) # 归一化
# 极化熵 H
H = -np.sum(eigenvalues * np.log(eigenvalues + 1e-12))
# 各向异性度
if len(eigenvalues) > 1:
A = (eigenvalues[1] - eigenvalues[2]) / (eigenvalues[1] + eigenvalues[2])
else:
A = 0
# 平均散射角
alpha = np.sum(eigenvalues * np.array([0, 45, 90])) # 简化计算
descriptors = {
'polarization_entropy': H,
'anisotropy': A,
'mean_scattering_angle': alpha,
'eigenvalues': eigenvalues,
'scattering_matrix': S
}
return descriptors
def generate_polarimetric_signature(self, target_type='dihedral'):
"""生成极化特征图"""
# 不同目标的典型散射矩阵
scattering_models = {
'sphere': np.array([[1, 0], [0, 1]]), # 球体
'dihedral': np.array([[1, 0], [0, -1]]), # 二面角
'dipole': np.array([[1, 0], [0, 0]]), # 偶极子
'helix': np.array([[1, 1j], [1j, 1]]) # 螺旋体
}
if target_type in scattering_models:
self.S_matrix = scattering_models[target_type]
# 生成全极化数据
polarimetric_data = self.simulate_full_polarimetric_data()
return polarimetric_data
class PerformanceOptimizer:
"""性能优化工具"""
@staticmethod
def vectorize_computations():
"""矢量化计算示例"""
# 不好的做法：循环
def slow_computation(N):
result = np.zeros(N)
for i in range(N):
result[i] = np.sin(i) * np.cos(i)
return result
# 好的做法：矢量化
def fast_computation(N):
i = np.arange(N)
return np.sin(i) * np.cos(i)
return fast_computation
@staticmethod
def use_numba_acceleration():
"""使用 Numba 加速"""
from numba import jit, njit
@njit(parallel=True)
def accelerated_pulse_compression(echo_signal, matched_filter):
"""加速的脉冲压缩"""
N = len(echo_signal)
M = len(matched_filter)
result = np.zeros(N, dtype=np.complex128)
for i in range(N):
for j in range(M):
if i - j >= 0:
result[i] += echo_signal[i - j] * matched_filter[j]
return result
return accelerated_pulse_compression
@staticmethod
def memory_optimization():
"""内存优化技巧"""
tips = """
内存优化技巧：
1. 使用 np.float32 而不是 np.float64（如果精度允许）
2. 及时释放不再需要的大数组：del large_array
3. 使用内存映射文件处理超大数组
4. 分批处理数据，避免一次加载所有数据
5. 使用稀疏矩阵存储稀疏数据
"""
return tips
@staticmethod
def gpu_acceleration():
"""GPU 加速示例（使用 CuPy）"""
try:
import cupy as cp
def gpu_pulse_compression(echo_signal, matched_filter):
# 传输到 GPU
echo_gpu = cp.asarray(echo_signal)
filter_gpu = cp.asarray(matched_filter)
# GPU 上的 FFT 卷积
N = len(echo_signal)
M = len(matched_filter)
fft_size = 2**int(np.ceil(np.log2(N + M - 1)))
echo_fft = cp.fft.fft(echo_gpu, fft_size)
filter_fft = cp.fft.fft(filter_gpu, fft_size)
result_gpu = cp.fft.ifft(echo_fft * filter_fft)[:N]
# 传输回 CPU
result = cp.asnumpy(result_gpu)
return result
return gpu_pulse_compression
except ImportError:
print("CuPy 未安装，无法使用 GPU 加速")
return None
class CodeQualityManager:
"""代码质量管理"""
@staticmethod
def unit_test_example():
"""单元测试示例"""
import unittest
class TestMonopulseSimulator(unittest.TestCase):
def setUp(self):
"""测试准备"""
self.simulator = MonopulseEchoSimulator()
def test_geometry_calculation(self):
"""测试几何计算"""
# 测试视线角计算
target_pos = np.array([10000, 0, 0])
radar_pos = np.array([0, 0, 0])
azimuth, elevation, R = self.simulator.geometry.compute_angles(
target_pos, radar_pos
)
self.assertAlmostEqual(azimuth, 0, places=5)
self.assertAlmostEqual(elevation, 0, places=5)
self.assertAlmostEqual(R, 10000, places=1)
def test_antenna_pattern(self):
"""测试天线方向图"""
# 测试波束中心增益为 1
gain = self.simulator.antenna.sum_pattern(0, 0)
self.assertAlmostEqual(gain, 1.0, places=2)
# 测试差波束中心为零
diff_gain = self.simulator.antenna.diff_pattern(0, 0)
self.assertAlmostEqual(diff_gain, 0.0, places=2)
def test_pulse_compression(self):
"""测试脉冲压缩"""
# 生成测试信号
test_signal = np.exp(1j * np.linspace(0, 10, 1000))
matched_filter = np.conj(test_signal[::-1])
compressed = pulse_compression(test_signal, matched_filter)
# 检查主瓣峰值
peak_idx = np.argmax(np.abs(compressed))
self.assertEqual(peak_idx, len(test_signal)-1)
return TestMonopulseSimulator
@staticmethod
def integration_test():
"""集成测试示例"""
def test_full_simulation_pipeline():
"""测试完整仿真流程"""
print("运行集成测试...")
# 创建仿真器
simulator = MonopulseEchoSimulator()
# 运行仿真
results = simulator.run_simulation()
# 检查结果
assert 'sum_channel' in results
assert 'diff_az' in results
assert 'diff_el' in results
assert 'time_axis' in results
# 检查数据维度一致性
sum_shape = results['sum_channel'].shape
diff_az_shape = results['diff_az'].shape
diff_el_shape = results['diff_el'].shape
assert sum_shape == diff_az_shape == diff_el_shape
assert len(results['time_axis']) == sum_shape[0]
print("集成测试通过！")
return True
return test_full_simulation_pipeline
@staticmethod
def performance_benchmark():
"""性能基准测试"""
import time
def benchmark_pulse_compression():
"""脉冲压缩性能测试"""
sizes = [1000, 5000, 10000, 50000]
results = {}
for size in sizes:
print(f"测试大小：{size}")
# 生成测试数据
echo_signal = np.random.randn(size) + 1j*np.random.randn(size)
matched_filter = np.random.randn(size//10) + 1j*np.random.randn(size//10)
# 测试 CPU 版本
start = time.time()
_ = pulse_compression(echo_signal, matched_filter)
cpu_time = time.time() - start
results[size] = {
'cpu_time': cpu_time,
'throughput': size / cpu_time # 每秒处理采样点数
}
print(f" CPU 时间：{cpu_time:.4f}秒")
print(f" 吞吐量：{results[size]['throughput']:.0f} 采样点/秒")
return results
return benchmark_pulse_compression
本文系统性地阐述了单脉冲雷达导引头回波生成的全过程，涵盖了从基础理论到工程实现的各个层面。通过本仿真系统的构建，我们实现了以下核心能力：
本仿真系统在以下场景中具有重要应用价值：
对于希望在实际项目中应用本技术的工程师，建议：
monopulse_simulator/ │ ├── README.md # 项目说明文档 ├── requirements.txt # 依赖包列表 ├── LICENSE # 开源许可证 ├── .gitignore # Git 忽略文件 ├── setup.py # 安装脚本 │ ├── config/ # 配置文件目录 │ ├── default_config.json # 默认配置 │ ├── high_fidelity_config.json # 高保真配置 │ ├── quick_test_config.json # 快速测试配置 │ ├── multi_target_config.json # 多目标场景配置 │ └── calibration_config.json # 校准测试配置 │ ├── src/ # 源代码目录 │ ├── init.py │ ├── core/ # 核心仿真模块 │ │ ├── init.py │ │ ├── simulator.py # 主仿真器类 │ │ ├── geometry.py # 几何计算模块 │ │ ├── antenna.py # 天线方向图模块 │ │ ├── waveform.py # 波形生成模块 │ │ ├── target.py # 目标模型模块 │ │ ├── propagator.py # 信号传播模块 │ │ └── coordinate_system.py # 坐标系管理 │ │ │ ├── effects/ # 效应建模模块 │ │ ├── init.py │ │ ├── channel_errors.py # 通道误差模型 │ │ ├── target_effects.py # 目标效应模型 │ │ ├── environment.py # 环境效应模型 │ │ ├── clutter.py # 杂波生成器 │ │ ├── jammer.py # 干扰生成器 │ │ └── glint.py # 角闪烁模型 │ │ │ ├── processing/ # 信号处理模块 │ │ ├── init.py │ │ ├── pulse_compression.py # 脉冲压缩处理 │ │ ├── beamforming.py # 波束形成算法 │ │ ├── doa_estimation.py # 波达方向估计 │ │ ├── cfar.py # CFAR 检测算法 │ │ └── tracking.py # 目标跟踪算法 │ │ │ ├── utils/ # 工具函数模块 │ │ ├── init.py │ │ ├── signal_processing.py # 信号处理工具 │ │ ├── math_utils.py # 数学工具函数 │ │ ├── file_io.py # 文件读写工具 │ │ ├── performance.py # 性能评估工具 │ │ └── validation.py # 验证工具函数 │ │ │ ├── visualization/ # 可视化模块 │ │ ├── init.py │ │ ├── plotter.py # 基本绘图工具 │ │ ├── radar_display.py # 雷达数据显示 │ │ ├── animation.py # 动画生成工具 │ │ └── dashboard.py # 交互式仪表板 │ │ │ └── interfaces/ # 外部接口模块 │ ├── init.py │ ├── ml_interface.py # 机器学习接口 │ ├── hardware_interface.py # 硬件接口 │ ├── data_export.py # 数据导出接口 │ └── api.py # Web API 接口 │ ├── examples/ # 示例脚本目录 │ ├── basic_simulation.py # 基础仿真示例 │ ├── multi_target_example.py # 多目标示例 │ ├── error_analysis.py # 误差分析示例 │ ├── performance_benchmark.py # 性能基准测试 │ ├── calibration_demo.py # 校准演示 │ ├── real_time_demo.py # 实时仿真演示 │ └── tutorial_notebooks/ # Jupyter 教程 │ ├── 01_basic_concepts.ipynb │ ├── 02_signal_generation.ipynb │ ├── 03_error_modeling.ipynb │ └── 04_advanced_features.ipynb │ ├── tests/ # 测试目录 │ ├── init.py │ ├── unit_tests/ # 单元测试 │ │ ├── test_geometry.py │ │ ├── test_antenna.py │ │ ├── test_waveform.py │ │ ├── test_channel_errors.py │ │ └── test_pulse_compression.py │ │ │ ├── integration_tests/ # 集成测试 │ │ ├── test_simulator_integration.py │ │ ├── test_data_pipeline.py │ │ └── test_performance.py │ │ │ ├── validation_tests/ # 验证测试 │ │ ├── test_s_curve.py │ │ ├── test_channel_coherence.py │ │ └── test_monte_carlo.py │ │ │ └── test_data/ # 测试数据 │ ├── reference_signals.npy │ ├── test_configs/ │ └── expected_results/ │ ├── docs/ # 文档目录 │ ├── index.md # 文档首页 │ ├── installation.md # 安装指南 │ ├── quick_start.md # 快速开始 │ ├── user_guide/ # 用户指南 │ │ ├── configuration.md │ │ ├── simulation_workflow.md │ │ ├── advanced_features.md │ │ └── troubleshooting.md │ │ │ ├── api_reference/ # API 参考 │ │ ├── core.md │ │ ├── effects.md │ │ ├── processing.md │ │ └── utils.md │ │ │ ├── theory/ # 理论文档 │ │ ├── monopulse_principles.md │ │ ├── signal_modeling.md │ │ ├── error_analysis.md │ │ └── performance_metrics.md │ │ │ └── tutorials/ # 教程文档 │ ├── basic_simulation.md │ ├── custom_models.md │ └── real_world_applications.md │ ├── data/ # 数据目录 │ ├── raw/ # 原始数据 │ ├── processed/ # 处理后的数据 │ ├── calibration/ # 校准数据 │ └── results/ # 仿真结果 │ ├── notebooks/ # 分析笔记本目录 │ ├── exploratory_analysis.ipynb # 探索性分析 │ ├── error_sensitivity.ipynb # 误差敏感性分析 │ ├── performance_comparison.ipynb # 性能比较 │ └── real_data_validation.ipynb # 实测数据验证 │ ├── scripts/ # 实用脚本目录 │ ├── batch_simulation.py # 批量仿真脚本 │ ├── parameter_sweep.py # 参数扫描脚本 │ ├── data_generator.py # 数据生成脚本 │ ├── report_generator.py # 报告生成脚本 │ └── deploy_service.py # 部署服务脚本 │ ├── benchmarks/ # 性能基准目录 │ ├── speed_benchmarks.py # 速度基准测试 │ ├── memory_benchmarks.py # 内存使用测试 │ ├── accuracy_benchmarks.py # 精度基准测试 │ └── scalability_tests.py # 可扩展性测试 │ └── deployment/ # 部署相关文件 ├── docker/ │ ├── Dockerfile │ └── docker-compose.yml ├── kubernetes/ │ ├── deployment.yaml │ └── service.yaml ├── cloud/ │ ├── aws_deploy.sh │ └── azure_config.json └── api_server/ # API 服务器 ├── app.py ├── requirements.txt └── config.yaml
{
"simulation": {
"mode": "high_fidelity",
"duration": 1.0,
"fs": 100e6,
"prf": 2000,
"num_pulses": 2000,
"random_seed": 42
},
"radar": {
"name": "X-band monopulse seeker",
"fc": 10e9,
"bandwidth": 20e6,
"pulse_width": 20e-6,
"beamwidth_az": 3.0,
"beamwidth_el": 3.0,
"antenna_type": "gaussian"
},
"target": {
"scenario": "single_moving_target",
"initial_range": 15000,
"radial_velocity": -250,
"cross_range_velocity": 50,
"rcs_mean_db": 10,
"rcs_type": "Swerling3",
"glint_enabled": true
}
}
# 主仿真器类结构
class MonopulseEchoSimulator:
def __init__(self, config=None):
self.config = config or self._load_default_config()
self._initialize_modules()
def _initialize_modules(self):
self.geometry = GeometryEngine(self.config['geometry'])
self.antenna = AntennaModel(self.config['antenna'])
self.waveform = WaveformGenerator(self.config['waveform'])
self.target = TargetModel(self.config['target'])
self.channel = ChannelModel(self.config['channel'])
self.effects = EffectsManager(self.config['effects'])
self.processor = SignalProcessor(self.config['processing'])
def simulate(self, scenario=None):
"""执行仿真"""
# 1. 准备场景
scenario = scenario or self._create_default_scenario()
# 2. 生成目标轨迹
trajectory = self.target.generate_trajectory(scenario)
# 3. 对每个时间步进行处理
results = self._process_time_steps(trajectory)
# 4. 后处理
results = self.processor.post_process(results)
return results
def calibrate(self):
"""系统校准"""
calibration_data = self._run_calibration_sequence()
calibration_params = self._estimate_calibration_params(calibration_data)
self._apply_calibration(calibration_params)
return calibration_params
# 标准数据格式
class RadarDataFormat:
"""定义标准雷达数据格式"""
@staticmethod
def create_data_structure():
return {
'metadata': {
'timestamp': '2024-01-01T00:00:00',
'config': {}, # 仿真配置
'version': '1.0.0'
},
'signals': {
'sum_channel': None, # 和通道数据 [pulses × samples]
'diff_az': None, # 方位差通道数据
'diff_el': None, # 俯仰差通道数据
'time_axis': None, # 时间轴
'range_axis': None, # 距离轴
'doppler_axis': None # 多普勒轴
},
'targets': [
{
'id': 1,
'type': 'aircraft',
'trajectory': None, # 目标轨迹
'parameters': {} # 目标参数
}
],
'processing': {
'pulse_compressed': False,
'filtered': False,
'detected': False,
'tracked': False
},
'metrics': {
'snr_db': None,
'range_resolution_m': None,
'doppler_resolution_hz': None,
'angle_accuracy_deg': None
}
}
微信公众号「极客日志」，在微信中扫描左侧二维码关注。展示文案：极客日志 zeeklog
使用加密算法（如AES、TripleDES、Rabbit或RC4）加密和解密文本明文。 在线工具，加密/解密文本在线工具，online
解析常见 curl 参数并生成 fetch、axios、PHP curl 或 Python requests 示例代码。 在线工具，curl 转代码在线工具，online
将字符串编码和解码为其 Base64 格式表示形式即可。 在线工具，Base64 字符串编码/解码在线工具，online
将字符串、文件或图像转换为其 Base64 表示形式。 在线工具，Base64 文件转换器在线工具，online
将 Markdown（GFM）转为 HTML 片段，浏览器内 marked 解析；与 HTML转Markdown 互为补充。 在线工具，Markdown转HTML在线工具，online
将 HTML 片段转为 GitHub Flavored Markdown，支持标题、列表、链接、代码块与表格等；浏览器内处理，可链接预填。 在线工具，HTML转Markdown在线工具，online