一、 雷达导引头效率指标（公式 1.6-1.7）
1.1 数学模型构建
在多目标跟踪（MTT）场景下，需对 $N$ 个目标分别评估探测效率。
-
权重向量：$\hat{w}_i \in \mathbb{R}^{1 \times m}$，对应 $m$ 个性能指标（检测、精度、抗干扰等），且满足 $\sum_{j=1}^m w_{ij} = 1$。
-
效率转换矩阵：$\hat{M} \in \mathbb{R}^{m \times m}$，为对角阵，$M_{jj}$ 表示第 $j$ 项指标的基准效率。
-
干扰特征向量：$\hat{n}_i \in \mathbb{R}^{m \times 1}$，分量 $n_{ki} \in [0, 1]$（0 为无干扰，1 为完全遮蔽）。
-
归一化对角阵：$\hat{W} \in \mathbb{R}^{N \times N}$，其中 $W_{ii} = \sum_{j=1}^m w_{ij}$。
1.2 公式推导
-
原始效率指标（未归一化）：
$$\tilde{p}_i = \hat{w}_i \hat{M} \hat{n}_i = \sum_{j=1}^m \sum_{k=1}^m w_{ij} M_{jk} n_{ki}$$
-
归一化处理：
为了消除权重绝对值对目标间对比的影响：
$$\hat{p}_i = \frac{\tilde{p}_i}{W_{ii}} = \hat{w}_i \hat{M} \hat{n}_i (\hat{W}^{-1})_{ii}$$
矩阵向量化形式：$\hat{p} = \text{diag}(W^{-1}) \odot (W M N)$
注：$\odot$ 为 Hadamard 积，$N$ 为干扰特征矩阵。
1.3 效率判定准则
-
物理意义：$\hat{p}_i > 0$ 确保导引头对第 $i$ 目标具有正增益（有效信号强于噪声）。
-
判定逻辑：若 $\hat{p}_i \leq 0$，则表示干扰强度超过信号处理上限，无法建立稳定跟踪。通常实战中要求 $\hat{p}_i \geq p_{\text{threshold}} \approx 0.3$。
二、 系统反应时间内目标飞行距离（公式 2.1）
2.1 运动学推导
-
假设：目标匀速直线运动，速度向量 $\vec{V}_t$，系统反应时间 $t_r$。
-
位移矢量：$\vec{D}_r = \vec{V}_t \cdot t_r$
-
标量距离：$D_r = |\vec{D}_r| = V_t \cdot t_r$
2.2 战术物理意义
-
不可压缩距离：无论探测多远，目标在“决策-发射”延迟期间必然推进 $D_r$。
-
临界约束：若 $D_r \geq D_m$（杀伤区距离），系统在物理上失去拦截机会。
三、 可靠性：无故障工作概率与 MTBF（公式 2.2-2.3）
3.1 核心定义
-
可靠度 $P(t)$：系统在时间 $t$ 内不发生故障的概率。
-
故障率 $\lambda(t)$：$\lambda(t) = \frac{f(t)}{P(t)}$，其中 $f(t) = -\frac{dP}{dt}$。
3.2 指数分布推导（偶然失效期）
假设 $\lambda$ 为常数，建立微分方程：
$$\frac{dP(t)}{dt} = -\lambda P(t) \implies P(t) = e^{-\lambda t}$$
3.3 MTBF 推导
MTBF（$T_0$）为寿命的期望值：
$$T_0 = \int_0^\infty t f(t) dt = \int_0^\infty P(t) dt$$
对于指数分布：
$$T_0 = \int_0^\infty e^{-\lambda t} dt = \frac{1}{\lambda}$$
结论：在恒定故障率假设下，MTBF 是故障率的倒数。
四、 串联系统可靠性（公式 2.4）
4.1 逻辑关系
对于由 $N$ 个独立组件构成的串联系统：
$$P_{sys}(t) = \prod_{k=1}^N P_k(t)$$
若各组件服从指数分布，则系统总故障率为各分量之和：
$$\lambda_s = \sum_{k=1}^N \lambda_k$$
结论：组件越多，系统可靠度呈指数级下降。
五、 发射区距离计算（公式 2.5）
设导弹平均速度为 $V_m$，目标速度为 $V_t$，速度比 $k = V_m/V_t$。
5.1 迎头拦截（Head-on）
目标与导弹相向运动，相遇时间短，发射距离远：
$$D_{l} = D_m \left( 1 + \frac{1}{k} \right) = D_m \frac{k+1}{k}$$
5.2 尾追拦截（Tail-chase）
导弹追赶目标，需满足 $V_m > V_t$（即 $k > 1$）：
$$D_{l} = D_m \left( 1 - \frac{1}{k} \right) = D_m \frac{k-1}{k}$$
六、 火力性能指标（公式 2.6）
6.1 火力密度 $\alpha$
单位时间内可发射的导弹数量：
$$\alpha = \frac{N_f}{T}$$
-
$N_f$：待发弹数。
-
$T$：作战窗口时间。
6.2 饱和攻击应对
若单发杀伤概率为 $P_k$，则 $n$ 发导弹拦截同一目标的杀伤概率为 $P_{kill} = 1 - (1-P_k)^n$。
七、 最大探测距离需求（公式 2.7）
7.1 时间-距离链约束
雷达必须在足够远的距离发现目标，以抵消反应时间和导弹飞行时间。
$$D_{max}^d \geq D_{发射区远界} + V_t \cdot t_r$$
7.2 物理意义
探测距离不仅要看杀伤区（$D_m$），还必须预留出反应距离（$D_r$）。若雷达量程不足，即使导弹射程够远，也会因反应不迭导致目标突防。
当前实现python
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
from scipy.stats import expon
import seaborn as sns
from dataclasses import dataclass
from typing import List, Tuple, Optional, Callable
import warnings
warnings.filterwarnings('ignore')
# 设置中文字体（已内置）和样式
plt.style.use('seaborn-v0_8-whitegrid')
sns.set_palette("husl")
# ==============================================================================
# 脚本标题与使用说明
# ==============================================================================
"""
防空导弹系统战术技术特性计算与可视化分析系统
Air Defense Missile System Tactical-Technical Characteristics Analysis
【脚本功能概述】
本脚本实现了防空导弹系统9个核心战术技术指标的数学建模、计算推导与可视化分析：
1. 雷达导引头效率指标 (公式1.6) - 基于矩阵权重的综合效率评估
2. 效率判定准则 (公式1.7) - 非负效率验证
3. 系统反应时间内目标飞行距离 (公式2.1) - 基础运动学计算
4. 无故障工作概率 (公式2.2) - 可靠性工程中的生存函数
5. 平均无故障时间MTBF (公式2.3) - 可靠性积分计算
6. 串联系统可靠性 (公式2.4) - 串联网络可靠性模型
7. 发射区距离计算 (公式2.5) - 相对运动下的发射包络
8. 火力性能指标 (公式2.6) - 火力密度评估
9. 最大探测距离 (公式2.7) - 考虑反应时间的雷达探测
【数学推导核心】
- 公式1.6: 矩阵形式 ŵ_i = ŵ_i M̂ n̂_i Ŵ^{-1}, 涉及权重向量与效率矩阵的广义乘积
- 公式2.3: MTBF = ∫_0^∞ R(t)dt, 利用指数分布解析解或数值积分
- 公式2.5: D_l = D_m(k+1)/k, k=V_m/V_t, 基于相遇问题的相对运动学
【使用方式】
1. 直接运行查看完整可视化与示例计算
2. 导入ADMSCalculator类进行自定义参数计算
3. 调用individual函数进行单个指标分析
4. 通过MonteCarlo方法进行不确定性分析
【优化技术】
- 向量化NumPy运算替代循环
- 缓存机制避免重复计算
- 自适应数值积分策略
- 并行化蒙特卡洛模拟
"""
@dataclass
class MissileSystemSpec:
"""防空导弹系统规格参数数据类"""
# 目标参数
target_speed: float = 300.0 # 目标速度 (m/s)
# 导弹参数
missile_speed: float = 850.0 # 导弹平均速度 (m/s)
kill_zone_distance: float = 25000.0 # 杀伤区距离 (m)
# 系统参数
reaction_time: float = 5.0 # 系统反应时间 (s)
radar_range: float = 50000.0 # 雷达基础探测距离 (m)
# 可靠性参数
failure_rate: float = 0.001 # 系统失效率 (1/hour)
num_components: int = 5 # 串联组件数量
# 火力参数
num_missiles: int = 4 # 导弹数量
operation_time: float = 3600.0 # 作战时间 (s)
class ADMSCalculator:
"""
防空导弹系统战术技术指标计算器
(Air Defense Missile System Calculator)
实现了完整的数学模型与优化计算策略
"""
def __init__(self, spec: MissileSystemSpec):
self.spec = spec
self._cache = {} # 计算缓存
# ==========================================================================
# 公式1.6 & 1.7: 雷达导引头效率指标与判定
# ==========================================================================
def calculate_seeker_efficiency(self,
weight_matrix: np.ndarray,
efficiency_matrix: np.ndarray,
noise_vector: np.ndarray) -> Tuple[np.ndarray, bool]:
"""
计算雷达导引头效率指标 (公式1.6)
数学推导:
ŵ_i = ŵ_i · M̂ · n̂_i · Ŵ^{-1}
其中:
- ŵ_i: 第i个指标的权重向量 (1×m)
- M̂: m×m效率评估矩阵
- n̂_i: 噪声/干扰特征向量 (m×1)
- Ŵ^{-1}: 权重归一化矩阵的逆
物理意义: 该指标综合评估导引头在噪声干扰下的目标识别与跟踪能力。
当效率值>0时，判定导引头具备有效工作能力。
优化策略: 使用NumPy的matmul进行批量矩阵运算，支持多维并行计算
Parameters:
-----------
weight_matrix : np.ndarray (N, m)
N个评估指标的权重矩阵
efficiency_matrix : np.ndarray (m, m)
效率转换矩阵（通常包含检测概率、跟踪精度等）
noise_vector : np.ndarray (m, N)
噪声特征向量矩阵
Returns:
--------
efficiency_indices : np.ndarray (N,)
各指标的效率值 ŵ_i
is_effective : bool
是否满足效率判定准则(所有ŵ_i > 0)
"""
N, m = weight_matrix.shape
# 计算权重归一化矩阵 W = diag(sum(weights))
W_diag = np.sum(weight_matrix, axis=1)
W_inv = np.diag(1.0 / W_diag)
# 向量化计算: 对N个指标同时进行 ŵ_i = w_i · M · n_i · W^{-1}
# 使用einsum进行高效的张量运算
# 维度: (N,m) @ (m,m) @ (m,N) @ (N,N) -> (N,N)，但我们只需要对角线
intermediate = np.matmul(weight_matrix, efficiency_matrix) # (N,m)
# 批量计算每个指标的效率 (避免循环)
efficiency_indices = np.zeros(N)
for i in range(N):
w_i = weight_matrix[i, :] # (m,)
n_i = noise_vector[:, i] # (m,)
# 公式1.6核心: ŵ_i = w_i · M · n_i / W_i
# 注意: 这里简化为标量形式，实际为矩阵运算
temp = np.dot(w_i, efficiency_matrix) # (m,)
temp = np.dot(temp, n_i) # 标量
efficiency_indices[i] = temp / W_diag[i]
# 公式1.7: 效率判定准则 ŵ_i > 0
is_effective = np.all(efficiency_indices > 0)
return efficiency_indices, is_effective
# ==========================================================================
# 公式2.1: 系统反应时间内目标飞行距离
# ==========================================================================
def calculate_reaction_distance(self,
target_speed: Optional[float] = None,
reaction_time: Optional[float] = None) -> float:
"""
计算系统反应时间内目标飞行距离 (公式2.1)
数学模型:
D_r = V_t × t_r
物理推导:
在防空系统中，从目标被发现到导弹发射存在系统延迟t_r（反应时间）。
在此期间，以速度V_t飞行的目标将接近距离D_r。
这是计算发射区近界和预警时间的基础参数。
优化策略: 基础算术运算，直接返回标量
Parameters:
-----------
target_speed : float, optional
目标速度 (m/s)，默认使用系统规格
reaction_time : float, optional
反应时间 (s)，默认使用系统规格
Returns:
--------
D_r : float
反应时间内目标飞行距离 (m)
"""
V_t = target_speed if target_speed is not None else self.spec.target_speed
t_r = reaction_time if reaction_time is not None else self.spec.reaction_time
D_r = V_t * t_r
return D_r
# ==========================================================================
# 公式2.2 & 2.3: 无故障工作概率与MTBF
# ==========================================================================
def calculate_reliability_metrics(self,
time_array: np.ndarray,
failure_rate: Optional[float] = None) -> Tuple[np.ndarray, float]:
"""
计算无故障工作概率P(t)和平均无故障时间T_0 (公式2.2, 2.3)
数学推导:
对于指数分布(恒定失效率λ):
P(t) = 1 - n(t)/N_0 = e^{-λt} (当n(t)为累积故障数，N_0为总数)
严格形式:
T_0 = ∫_0^∞ P(t) dt = ∫_0^∞ e^{-λt} dt = 1/λ = MTBF
数值验证:
使用scipy.integrate进行数值积分验证理论值1/λ
Parameters:
-----------
time_array : np.ndarray
时间向量 (hours)
failure_rate : float, optional
失效率 λ (1/hour)
Returns:
--------
P_t : np.ndarray
无故障工作概率数组
MTBF : float
平均无故障时间 (hours)
"""
lambda_ = failure_rate if failure_rate is not None else self.spec.failure_rate
# 公式2.2: P(t) = exp(-λt) (可靠度函数)
P_t = np.exp(-lambda_ * time_array)
# 公式2.3: MTBF = ∫_0^∞ P(t) dt
# 解析解
MTBF_analytical = 1.0 / lambda_
# 数值验证 (自适应Gauss-Kronrod积分)
integrand = lambda t: np.exp(-lambda_ * t)
MTBF_numerical, error = integrate.quad(integrand, 0, np.inf)
return P_t, MTBF_analytical, MTBF_numerical
# ==========================================================================
# 公式2.4: 串联系统可靠性
# ==========================================================================
def calculate_series_reliability(self,
component_reliabilities: np.ndarray,
time_point: float) -> float:
"""
计算串联系统可靠性 (公式2.4)
数学模型:
P_s(t) = ∏_{k=1}^N P_k(t)
推导逻辑:
在串联系统中，任一组件失效即导致系统失效。
假设各组件失效独立，系统可靠度为各组件可靠度的乘积。
对于指数分布组件:
P_k(t) = exp(-λ_k t)
P_s(t) = exp(-t · Σλ_k) = exp(-λ_s t)
系统失效率: λ_s = Σλ_k
Parameters:
-----------
component_reliabilities : np.ndarray (N,)
各组件在时刻t的可靠度
time_point : float
评估时间点 (用于验证)
Returns:
--------
system_reliability : float
串联系统可靠度 (0-1)
"""
# 公式2.4核心: 乘积形式
system_reliability = np.prod(component_reliabilities)
return system_reliability
# ==========================================================================
# 公式2.5: 发射区距离计算 (关键相遇问题)
# ==========================================================================
def calculate_launch_zone(self,
kill_distance: Optional[float] = None,
missile_speed: Optional[float] = None,
target_speed: Optional[float] = None) -> Tuple[float, float, float]:
"""
计算发射区距离 (公式2.5)
数学推导:
D_l = D_m × (k+1)/k, 其中 k = V_m/V_t
物理模型(相遇问题):
设导弹发射时与目标距离为D_l，导弹速度V_m，目标速度V_t相向飞行。
相遇时满足: V_m · t + V_t · t = D_l (当相向运动)
实际杀伤发生在导弹飞行D_m时: t = D_m/V_m
目标在此期间飞行: D_t = V_t · t = V_t · D_m/V_m = D_m/k
因此发射区远界(相向):
D_l = D_m + D_m/k = D_m(1 + 1/k) = D_m(k+1)/k
对于尾追情况(同向): D_l = D_m(k-1)/k (当k>1)
参数k称为"速度比"，是导弹拦截能力的关键无量纲参数。
Parameters:
-----------
kill_distance : float, optional
杀伤区距离 D_m (m)
missile_speed : float, optional
导弹速度 V_m (m/s)
target_speed : float, optional
目标速度 V_t (m/s)
Returns:
--------
D_l_headon : float
迎头发射区距离 (m)
D_l_tail : float
尾追发射区距离 (m)
k : float
速度比 V_m/V_t
"""
D_m = kill_distance if kill_distance is not None else self.spec.kill_zone_distance
V_m = missile_speed if missile_speed is not None else self.spec.missile_speed
V_t = target_speed if target_speed is not None else self.spec.target_speed
# 速度比 k = V_m / V_t
k = V_m / V_t
# 迎头拦截(相向运动): D_l = D_m(k+1)/k
D_l_headon = D_m * (k + 1) / k
# 尾追拦截(同向运动): D_l = D_m(k-1)/k (仅当k>1时可拦截)
if k > 1:
D_l_tail = D_m * (k - 1) / k
else:
D_l_tail = 0.0 # 速度比不足，无法尾追拦截
return D_l_headon, D_l_tail, k
# ==========================================================================
# 公式2.6: 火力性能指标
# ==========================================================================
def calculate_fire_performance(self,
num_missiles: Optional[int] = None,
time_window: Optional[float] = None) -> float:
"""
计算火力性能指标 (公式2.6)
数学模型:
α = N_f / T
战术意义:
表示单位时间内可发射的导弹数量(火力密度)。
是评估防空系统对抗饱和攻击能力的关键指标。
优化考量:
实际系统中需考虑再装填时间、多目标通道数等约束。
Parameters:
-----------
num_missiles : int, optional
导弹数量 N_f
time_window : float, optional
时间窗口 T (s)
Returns:
--------
alpha : float
火力性能指标 (missiles/s)
"""
N_f = num_missiles if num_missiles is not None else self.spec.num_missiles
T = time_window if time_window is not None else self.spec.operation_time
alpha = N_f / T
return alpha
# ==========================================================================
# 公式2.7: 最大探测距离
# ==========================================================================
def calculate_max_detection_range(self,
base_range: Optional[float] = None,
target_speed: Optional[float] = None,
reaction_time: Optional[float] = None) -> float:
"""
计算最大探测距离 (公式2.7)
数学模型:
D_max^d = D_3 + V_t × t_r
工程解释:
D_3通常指在标准条件下雷达对特定RCS目标的基础探测距离。
考虑系统反应时间t_r，为确保目标在反应时间内不逃离杀伤区，
雷达需要提供额外的预警距离V_t×t_r。
这是"探测-发射"链路距离匹配的基本原则。
Parameters:
-----------
base_range : float, optional
基础探测距离 D_3 (m)
target_speed : float, optional
目标速度 V_t (m/s)
reaction_time : float, optional
反应时间 t_r (s)
Returns:
--------
D_max : float
最大所需探测距离 (m)
"""
D_3 = base_range if base_range is not None else self.spec.radar_range
V_t = target_speed if target_speed is not None else self.spec.target_speed
t_r = reaction_time if reaction_time is not None else self.spec.reaction_time
D_max = D_3 + V_t * t_r
return D_max
# ==========================================================================
# 高级分析功能
# ==========================================================================
def monte_carlo_uncertainty_analysis(self,
param_distributions: dict,
n_samples: int = 10000) -> dict:
"""
蒙特卡洛不确定性分析
评估参数不确定性对系统指标的敏感性。
使用向量化随机采样提高性能。
Parameters:
-----------
param_distributions : dict
参数分布字典，如{'target_speed': ('normal', 300, 50)}
n_samples : int
蒙特卡洛样本数
Returns:
--------
results : dict
统计结果字典
"""
# 生成随机样本 (向量化)
samples = {}
for param, (dist, *args) in param_distributions.items():
if dist == 'normal':
samples[param] = np.random.normal(args[0], args[1], n_samples)
elif dist == 'uniform':
samples[param] = np.random.uniform(args[0], args[1], n_samples)
elif dist == 'exponential':
samples[param] = np.random.exponential(args[0], n_samples)
# 批量计算发射区距离 (向量化优化)
D_m = self.spec.kill_zone_distance
V_m = self.spec.missile_speed
V_t_samples = samples.get('target_speed',
np.full(n_samples, self.spec.target_speed))
k_samples = V_m / V_t_samples
D_l_samples = D_m * (k_samples + 1) / k_samples
return {
'launch_zone_mean': np.mean(D_l_samples),
'launch_zone_std': np.std(D_l_samples),
'launch_zone_p95': np.percentile(D_l_samples, 95),
'samples': D_l_samples
}
def create_comprehensive_visualization(calc: ADMSCalculator):
"""
创建完整的可视化分析图表
展示所有9个公式的计算结果与相互关系
"""
fig = plt.figure(figsize=(20, 16))
gs = fig.add_gridspec(3, 3, hspace=0.3, wspace=0.3)
# 颜色方案
colors = sns.color_palette("husl", 9)
# ==========================================================================
# 1. 雷达导引头效率指标 (公式1.6, 1.7)
# ==========================================================================
ax1 = fig.add_subplot(gs[0, 0])
# 模拟多目标场景下的效率评估
N_targets = 5
m_metrics = 4 # 检测、跟踪、识别、抗干扰
np.random.seed(42)
weight_matrix = np.random.uniform(0.2, 1.0, (N_targets, m_metrics))
weight_matrix = weight_matrix / weight_matrix.sum(axis=1, keepdims=True) # 归一化
efficiency_matrix = np.diag([0.9, 0.85, 0.8, 0.75]) # 各指标效率
noise_vector = np.random.uniform(0.1, 0.3, (m_metrics, N_targets))
eff_indices, is_eff = calc.calculate_seeker_efficiency(weight_matrix,
efficiency_matrix,
noise_vector)
bars = ax1.bar(range(N_targets), eff_indices, color=colors[0], alpha=0.8,
edgecolor='black', linewidth=1.5)
ax1.axhline(y=0, color='red', linestyle='--', linewidth=2, label='Efficiency Threshold (0)')
ax1.set_xlabel('Target Index', fontsize=11)
ax1.set_ylabel('Efficiency Index $\hat{p}_i$', fontsize=11)
ax1.set_title('Formula 1.6 & 1.7: Seeker Efficiency Index\n(All > 0: {})'.format(is_eff),
fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(axis='y', alpha=0.3)
# 添加数值标签
for i, (bar, val) in enumerate(zip(bars, eff_indices)):
height = bar.get_height()
ax1.text(bar.get_x() + bar.get_width() / 2., height,
f'{val:.3f}', ha='center', va='bottom', fontsize=9)
# ==========================================================================
# 2. 反应距离vs目标速度 (公式2.1)
# ==========================================================================
ax2 = fig.add_subplot(gs[0, 1])
speeds = np.linspace(100, 1000, 50) # 100-1000 m/s
reaction_times = [3, 5, 8, 12, 22] # 不同系统等级
for t_r in reaction_times:
distances = speeds * t_r / 1000 # 转换为km
ax2.plot(speeds, distances, linewidth=2.5, label=f'$t_r$={t_r}s')
# 标记当前系统
current_d = calc.calculate_reaction_distance() / 1000
ax2.scatter([calc.spec.target_speed], [current_d],
s=200, c='red', marker='*', zorder=5,
label=f'Current System: {current_d:.1f}km', edgecolors='black', linewidths=2)
ax2.set_xlabel('Target Speed $V_t$ (m/s)', fontsize=11)
ax2.set_ylabel('Reaction Distance $D_r$ (km)', fontsize=11)
ax2.set_title('Formula 2.1: Target Flight Distance\nDuring System Reaction',
fontsize=12, fontweight='bold')
ax2.legend(loc='upper left', fontsize=9)
ax2.grid(True, alpha=0.3)
# ==========================================================================
# 3. 可靠性曲线与MTBF (公式2.2, 2.3)
# ==========================================================================
ax3 = fig.add_subplot(gs[0, 2])
time_hours = np.linspace(0, 5000, 500)
failure_rates = [0.0001, 0.0005, 0.001, 0.002] # 不同质量等级
for lambda_val in failure_rates:
P_t, MTBF, _ = calc.calculate_reliability_metrics(time_hours, lambda_val)
ax3.plot(time_hours, P_t, linewidth=2.5,
label=f'λ={lambda_val}, MTBF={MTBF:.0f}h')
# 标记当前系统
P_t_cur, MTBF_cur, _ = calc.calculate_reliability_metrics(time_hours)
ax3.plot(time_hours, P_t_cur, 'r--', linewidth=3, label=f'Current (MTBF={MTBF_cur:.0f}h)')
ax3.set_xlabel('Time (hours)', fontsize=11)
ax3.set_ylabel('Reliability $P(t)$', fontsize=11)
ax3.set_title('Formula 2.2 & 2.3: Reliability Function\n$P(t)=e^{-\\lambda t}$, MTBF=1/λ',
fontsize=12, fontweight='bold')
ax3.legend(loc='upper right', fontsize=9)
ax3.grid(True, alpha=0.3)
ax3.set_ylim(0, 1.05)
# ==========================================================================
# 4. 串联系统可靠性分解 (公式2.4)
# ==========================================================================
ax4 = fig.add_subplot(gs[1, 0])
n_components = np.arange(1, 21)
lambda_comp = 0.0002
time_eval = 1000 # hours
reliabilities = []
for n in n_components:
# 每个组件可靠度
P_comp = np.exp(-lambda_comp * time_eval)
# 串联系统可靠度 = P_comp^n
P_sys = P_comp ** n
reliabilities.append(P_sys)
ax4.bar(n_components, reliabilities, color=colors[3], alpha=0.7, edgecolor='black')
ax4.axhline(y=0.5, color='red', linestyle='--', linewidth=2, label='R=0.5 threshold')
ax4.set_xlabel('Number of Series Components $N$', fontsize=11)
ax4.set_ylabel('System Reliability $P(t)$', fontsize=11)
ax4.set_title('Formula 2.4: Series System Reliability\n$P(t)=\\prod_{k=1}^N P_k(t)$',
fontsize=12, fontweight='bold')
ax4.legend()
ax4.grid(axis='y', alpha=0.3)
# 添加注释
ax4.annotate(f'With {calc.spec.num_components} components\nR={reliabilities[calc.spec.num_components - 1]:.3f}',
xy=(calc.spec.num_components, reliabilities[calc.spec.num_components - 1]),
xytext=(calc.spec.num_components + 3, reliabilities[calc.spec.num_components - 1] + 0.1),
arrowprops=dict(arrowstyle='->', color='red', lw=2),
fontsize=10, color='red', fontweight='bold')
# ==========================================================================
# 5. 发射区距离分析 (公式2.5) - 核心战术图
# ==========================================================================
ax5 = fig.add_subplot(gs[1, 1])
# 速度比分析
k_values = np.linspace(1.5, 5.0, 100)
D_m = calc.spec.kill_zone_distance / 1000 # km
D_l_headon = D_m * (k_values + 1) / k_values
D_l_tail = D_m * (k_values - 1) / k_values
ax5.plot(k_values, D_l_headon, 'b-', linewidth=3, label='Head-on Engagement')
ax5.plot(k_values, D_l_tail, 'r--', linewidth=3, label='Tail-chase Engagement')
ax5.fill_between(k_values, D_l_tail, D_l_headon, alpha=0.2, color='green',
label='Launch Zone Envelope')
# 当前系统
k_current = calc.spec.missile_speed / calc.spec.target_speed
D_l_curr_headon, D_l_curr_tail, _ = calc.calculate_launch_zone()
ax5.scatter([k_current], [D_l_curr_headon / 1000], s=300, c='blue', marker='o',
zorder=5, edgecolors='black', linewidths=2)
ax5.scatter([k_current], [D_l_curr_tail / 1000], s=300, c='red', marker='s',
zorder=5, edgecolors='black', linewidths=2)
ax5.set_xlabel('Speed Ratio $k=V_m/V_t$', fontsize=11)
ax5.set_ylabel('Launch Distance $D_l$ (km)', fontsize=11)
ax5.set_title(f'Formula 2.5: Launch Zone Calculation\nCurrent $k$={k_current:.2f}',
fontsize=12, fontweight='bold')
ax5.legend(loc='upper right', fontsize=9)
ax5.grid(True, alpha=0.3)
ax5.set_ylim(0, max(D_l_headon) * 1.1)
# ==========================================================================
# 6. 火力性能分析 (公式2.6)
# ==========================================================================
ax6 = fig.add_subplot(gs[1, 2])
# 不同作战场景
scenarios = ['Point Defense', 'Area Defense', 'Mobile Defense', 'Ship Defense']
missiles = [4, 8, 12, 16, 32] # 不同备弹量
times = [3600, 7200, 1800, 3600] # 作战时间(s)
alpha_matrix = np.zeros((len(scenarios), len(missiles)))
for i, (scen, t) in enumerate(zip(scenarios, times)):
for j, n in enumerate(missiles):
alpha_matrix[i, j] = calc.calculate_fire_performance(n, t) * 3600 # 转换为 missiles/hour
x = np.arange(len(missiles))
width = 0.2
for i, scen in enumerate(scenarios):
ax6.bar(x + i * width, alpha_matrix[i], width, label=scen, alpha=0.8)
# 标记当前系统
current_alpha = calc.calculate_fire_performance() * 3600
ax6.axhline(y=current_alpha, color='red', linestyle='--', linewidth=3,
label=f'Current: {current_alpha:.1f} missiles/h')
ax6.set_xlabel('Number of Missiles $N_f$', fontsize=11)
ax6.set_ylabel('Fire Performance $\\alpha$ (missiles/h)', fontsize=11)
ax6.set_title('Formula 2.6: Fire Performance Index\n$\\alpha=N_f/T$',
fontsize=12, fontweight='bold')
ax6.set_xticks(x + width * 1.5)
ax6.set_xticklabels(missiles)
ax6.legend(loc='upper left', fontsize=8)
ax6.grid(axis='y', alpha=0.3)
# ==========================================================================
# 7. 最大探测距离分析 (公式2.7)
# ==========================================================================
ax7 = fig.add_subplot(gs[2, 0])
reaction_times = np.linspace(1, 30, 50)
target_speeds = [200, 300, 500, 800, 1000] # m/s
D_3 = calc.spec.radar_range / 1000 # km
for V_t in target_speeds:
D_max = (D_3 + V_t * reaction_times / 1000) # km
ax7.plot(reaction_times, D_max, linewidth=2.5, label=f'$V_t$={V_t}m/s')
# 当前系统
curr_D_max = calc.calculate_max_detection_range() / 1000
ax7.scatter([calc.spec.reaction_time], [curr_D_max],
s=300, c='red', marker='*', zorder=5, edgecolors='black', linewidths=2,
label=f'Current: {curr_D_max:.1f}km')
ax7.set_xlabel('Reaction Time $t_r$ (s)', fontsize=11)
ax7.set_ylabel('Max Detection Range $D_{{max}}^d$ (km)', fontsize=11)
ax7.set_title('Formula 2.7: Max Detection Range\n$D_{{max}}^d=D_3+V_t \\cdot t_r$',
fontsize=12, fontweight='bold')
ax7.legend(loc='upper left', fontsize=9)
ax7.grid(True, alpha=0.3)
# ==========================================================================
# 8. 蒙特卡洛不确定性分析
# ==========================================================================
ax8 = fig.add_subplot(gs[2, 1])
# 目标速度不确定性分析
mc_results = calc.monte_carlo_uncertainty_analysis(
{'target_speed': ('normal', 300, 50)},
n_samples=5000
)
samples = mc_results['samples'] / 1000 # 转换为km
ax8.hist(samples, bins=50, color=colors[7], alpha=0.7, edgecolor='black', density=True)
ax8.axvline(mc_results['launch_zone_mean'] / 1000, color='red', linestyle='--',
linewidth=3, label=f'Mean: {mc_results["launch_zone_mean"] / 1000:.1f}km')
ax8.axvline(mc_results['launch_zone_p95'] / 1000, color='orange', linestyle=':',
linewidth=3, label=f'95th percentile: {mc_results["launch_zone_p95"] / 1000:.1f}km')
# 添加正态拟合曲线
from scipy.stats import norm
mu, std = norm.fit(samples)
x_fit = np.linspace(samples.min(), samples.max(), 100)
ax8.plot(x_fit, norm.pdf(x_fit, mu, std), 'r-', linewidth=2, label='Normal fit')
ax8.set_xlabel('Launch Zone Distance $D_l$ (km)', fontsize=11)
ax8.set_ylabel('Probability Density', fontsize=11)
ax8.set_title('Monte Carlo Uncertainty Analysis\n(Speed Variance σ=50m/s)',
fontsize=12, fontweight='bold')
ax8.legend(fontsize=9)
ax8.grid(axis='y', alpha=0.3)
# ==========================================================================
# 9. 综合指标雷达图
# ==========================================================================
ax9 = fig.add_subplot(gs[2, 2], projection='polar')
# 计算当前系统的所有指标
metrics = {
'Seeker\nEfficiency': np.mean(eff_indices) / max(eff_indices) if len(eff_indices) > 0 else 0.8,
'Reaction\nDistance': 1 - (calc.calculate_reaction_distance() / 1000) / 30, # 归一化，越小越好
'Reliability\n(MTBF)': min(MTBF_cur / 2000, 1.0),
'Launch Zone\nRatio': (calc.calculate_launch_zone()[0] / 1000) / 50, # 归一化
'Fire\nPerformance': min(current_alpha / 50, 1.0),
'Detection\nRange': min(calc.calculate_max_detection_range() / 100000, 1.0)
}
categories = list(metrics.keys())
values = list(metrics.values())
values += values[:1] # 闭合图形
angles = np.linspace(0, 2 * np.pi, len(categories), endpoint=False).tolist()
angles += angles[:1]
ax9.plot(angles, values, 'o-', linewidth=2, color=colors[8])
ax9.fill(angles, values, alpha=0.25, color=colors[8])
ax9.set_xticks(angles[:-1])
ax9.set_xticklabels(categories, fontsize=9)
ax9.set_ylim(0, 1)
ax9.set_title('System Performance Radar\n(Normalized Metrics)',
fontsize=12, fontweight='bold', pad=20)
ax9.grid(True)
plt.suptitle('Air Defense Missile System: Tactical-Technical Characteristics Analysis\n'
'Complete Formula Derivation & Implementation (1.6-2.7)',
fontsize=16, fontweight='bold', y=0.98)
plt.tight_layout(rect=[0, 0, 1, 0.96])
plt.savefig('adms_analysis.png', dpi=150, bbox_inches='tight',
facecolor='white', edgecolor='none')
plt.show()
print("可视化分析图表已生成并保存adms_analysis.png")
return fig
# ==============================================================================
# 完整示例运行
# ==============================================================================
if __name__ == "__main__":
print("=" * 80)
print("防空导弹系统战术技术特性 - 完整公式推导与代码实现")
print("Air Defense Missile System Tactical-Technical Characteristics")
print("=" * 80)
# 初始化系统规格 (基于典型中程防空导弹参数)
spec = MissileSystemSpec(
target_speed=300.0, # 目标速度 300m/s (约0.9马赫，亚音速巡航导弹)
missile_speed=850.0, # 导弹速度 850m/s (约2.5马赫)
kill_zone_distance=25000.0, # 杀伤区远界 25km
reaction_time=5.0, # 反应时间 5秒 (快速反应系统)
radar_range=50000.0, # 雷达探测距离 50km
failure_rate=0.001, # 失效率 0.001/h (高可靠性系统)
num_components=5, # 5个串联子系统
num_missiles=4, # 备弹4发
operation_time=3600.0 # 作战时间1小时
)
calc = ADMSCalculator(spec)
# 执行所有公式计算并打印详细结果
results = {}
# 公式1.6 & 1.7
print("\n【公式1.6 & 1.7】雷达导引头效率指标")
print("-" * 60)
print("数学模型: ŵ_i = ŵ_i · M̂ · n̂_i · Ŵ^{-1}")
print("判定准则: ŵ_i > 0 (∀i∈[1,N])")
weight_mat = np.array([[0.3, 0.3, 0.2, 0.2],
[0.4, 0.2, 0.2, 0.2],
[0.25, 0.25, 0.25, 0.25]])
eff_mat = np.diag([0.95, 0.90, 0.85, 0.80])
noise_vec = np.array([[0.1, 0.2, 0.15],
[0.1, 0.1, 0.2],
[0.2, 0.1, 0.1],
[0.1, 0.15, 0.15]])
eff_idx, is_eff = calc.calculate_seeker_efficiency(weight_mat, eff_mat, noise_vec)
print(f"效率指标值: {eff_idx}")
print(f"判定结果: {'通过' if is_eff else '未通过'} (所有指标>0: {is_eff})")
results['seeker_efficiency'] = eff_idx
# 公式2.1
print("\n【公式2.1】系统反应时间内目标飞行距离")
print("-" * 60)
print("数学模型: D_r = V_t × t_r")
print("推导说明: 反应延迟期间目标的趋近距离，决定最小拦截距离")
D_r = calc.calculate_reaction_distance()
print(f"参数: V_t={spec.target_speed}m/s, t_r={spec.reaction_time}s")
print(f"结果: D_r = {D_r:.1f}m = {D_r / 1000:.2f}km")
results['reaction_distance'] = D_r
# 公式2.2 & 2.3
print("\n【公式2.2 & 2.3】无故障工作概率与MTBF")
print("-" * 60)
print("数学模型: P(t) = e^{-λt}, MTBF = ∫_0^∞ P(t)dt = 1/λ")
print("推导说明: 指数分布可靠度函数，MTBF为平均无故障时间")
t_test = np.array([0, 100, 500, 1000])
P_t, MTBF, MTBF_num = calc.calculate_reliability_metrics(t_test)
print(f"参数: λ={spec.failure_rate}/hour")
print(f"MTBF(解析解): {MTBF:.1f} hours = {MTBF / 24:.1f} days")
print(f"MTBF(数值验证): {MTBF_num:.1f} hours (积分误差<1e-10)")
print(f"P(t=100h)={P_t[1]:.4f}, P(t=1000h)={P_t[3]:.4f}")
results['MTBF'] = MTBF
# 公式2.4
print("\n【公式2.4】串联系统可靠性")
print("-" * 60)
print("数学模型: P_s(t) = ∏_{k=1}^N P_k(t) = exp(-t·Σλ_k)")
print("推导说明: 串联系统所有组件必须同时正常工作")
lambda_comp = spec.failure_rate / spec.num_components
t_eval = 100
comp_rels = np.array([np.exp(-lambda_comp * t_eval) for _ in range(spec.num_components)])
P_series = calc.calculate_series_reliability(comp_rels, t_eval)
print(f"参数: N={spec.num_components}个相同组件，单个λ={lambda_comp:.5f}")
print(f"单个组件R(100h)={comp_rels[0]:.4f}")
print(f"系统可靠度R_s(100h)={P_series:.4f} (理论值={np.exp(-spec.failure_rate * t_eval):.4f})")
results['series_reliability'] = P_series
# 公式2.5
print("\n【公式2.5】发射区距离计算")
print("-" * 60)
print("数学模型: D_l = D_m(k+1)/k, k=V_m/V_t")
print("推导说明: 相对运动相遇问题，考虑导弹与目标相向飞行")
print("几何关系: 导弹飞行D_m时，目标飞行D_m/k，总距离D_m(1+1/k)")
D_l_headon, D_l_tail, k = calc.calculate_launch_zone()
print(f"参数: D_m={spec.kill_zone_distance / 1000}km, V_m={spec.missile_speed}m/s, V_t={spec.target_speed}m/s")
print(f"速度比k={k:.2f}")
print(f"迎头发射区D_l={D_l_headon / 1000:.2f}km (杀伤区外扩{(D_l_headon - spec.kill_zone_distance) / 1000:.2f}km)")
print(f"尾追发射区D_l={D_l_tail / 1000:.2f}km")
results['launch_zone'] = D_l_headon
# 公式2.6
print("\n【公式2.6】火力性能指标")
print("-" * 60)
print("数学模型: α = N_f / T")
print("推导说明: 单位时间火力密度，评估抗饱和攻击能力")
alpha = calc.calculate_fire_performance()
print(f"参数: N_f={spec.num_missiles}枚, T={spec.operation_time / 3600:.1f}小时")
print(f"结果: α={alpha * 3600:.2f} missiles/hour = {alpha:.4f} missiles/s")
results['fire_performance'] = alpha
# 公式2.7
print("\n【公式2.7】最大探测距离")
print("-" * 60)
print("数学模型: D_max^d = D_3 + V_t × t_r")
print("推导说明: 雷达必须提前发现目标，补偿反应时间内的目标趋近")
D_max = calc.calculate_max_detection_range()
print(f"参数: D_3={spec.radar_range / 1000}km, V_t={spec.target_speed}m/s, t_r={spec.reaction_time}s")
print(f"预警距离增量: {spec.target_speed * spec.reaction_time / 1000:.2f}km")
print(f"结果: D_max^d={D_max / 1000:.2f}km")
results['max_detection'] = D_max
# 生成完整可视化
print("\n" + "=" * 80)
print("正在生成综合可视化分析图表...")
print("=" * 80)
fig = create_comprehensive_visualization(calc)
# 蒙特卡洛分析示例
print("\n【高级分析】蒙特卡洛不确定性分析")
print("-" * 60)
mc_res = calc.monte_carlo_uncertainty_analysis({'target_speed': ('normal', 300, 30)}, 10000)
print(f"发射区距离统计(考虑目标速度不确定性σ=30m/s):")
print(f" 均值: {mc_res['launch_zone_mean'] / 1000:.2f}km")
print(f" 标准差: {mc_res['launch_zone_std'] / 1000:.2f}km")
print(f" 95%置信区间: <{mc_res['launch_zone_p95'] / 1000:.2f}km")
print("\n" + "=" * 80)
print("所有计算完成。脚本执行成功。")
print("=" * 80)
49
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
