火箭仿真系列-详细的大气与风场模型模块介绍

详细的大气与风场模型模块详解

大气与风场模型模块是其高精度弹道仿真的核心基础，提供了从标准大气到真实气象数据的全面环境模拟能力。以下是该模块的深入解析：

一、模块架构概览

1.1 核心类层次结构

Environment (主类)
├── AtmosphereModel (大气模型)
│ ├── StandardAtmosphere (标准大气模型)
│ ├── CustomAtmosphere (自定义大气模型)
│ ├── COSPAAtmosphere (COSPA大气模型)
│ └── RealAtmosphere (真实大气模型)
├── WindModel (风场模型)
│ ├── ConstantWind (恒定风)
│ ├── WindShear (风切变)
│ ├── WindGust (阵风)
│ ├── TurbulenceModel (湍流模型)
│ └── RealWind (真实风场)
├── EarthModel (地球模型)
│ └── WGS84 (WGS84地球模型)
└── WeatherForecast (天气预报)
├── NOAA_GFS (NOAA GFS)
├── ECMWF (欧洲中期天气预报)
├── RAP (快速更新)
└── CustomWeather (自定义天气)

二、大气模型详解

2.1 StandardAtmosphere (标准大气模型)

US Standard Atmosphere 1976 实现：

class USStandardAtmosphere1976:
"""美国标准大气 1976"""

def __init__(self):
# 大气分层定义 (海拔高度, 温度梯度) (m, K/m)
self.layers = [
(0, -0.0065), # 对流层: 0-11km
(11000, 0.0), # 对流层顶: 11-20km
(20000, 0.001), # 平流层: 20-32km
(32000, 0.0028), # 平流层: 32-47km
(47000, 0.0), # 平流层顶: 47-51km
(51000, -0.0028), # 中间层: 51-71km
(71000, -0.002), # 中间层: 71-86km
]

# 海平面参考值
self.sea_level = {
'temperature': 288.15, # K
'pressure': 101325.0, # Pa
'density': 1.225, # kg/m³
'speed_of_sound': 340.294, # m/s
}

def compute_at_altitude(self, altitude):
"""计算指定高度的大气参数"""

# 确保高度非负
h = max(altitude, 0)

# 确定所在层
layer_idx = 0
for i, (base_height, _) in enumerate(self.layers):
if h >= base_height:
layer_idx = i
else:
break

# 获取层参数
h_base, lapse_rate = self.layers[layer_idx]

if layer_idx > 0:
# 计算从海平面到该层底部的累积效应
T, P = self._integrate_to_base(layer_idx, h_base)
else:
T = self.sea_level['temperature']
P = self.sea_level['pressure']

# 计算该层内的温度
if lapse_rate == 0:
# 等温层
T_layer = T
P_layer = P * np.exp(-9.80665 * (h – h_base) /
(287.05 * T_layer))
else:
# 变温层
T_layer = T + lapse_rate * (h – h_base)
P_layer = P * (T_layer / T) ** (-9.80665 /
(lapse_rate * 287.05))

# 计算密度
rho = P_layer / (287.05 * T_layer)

# 计算声速
a = np.sqrt(1.4 * 287.05 * T_layer)

return {
'temperature': T_layer, # K
'pressure': P_layer, # Pa
'density': rho, # kg/m³
'speed_of_sound': a, # m/s
'dynamic_viscosity': self.compute_viscosity(T_layer),
'kinematic_viscosity': self.compute_kinematic_viscosity(T_layer, rho)
}

def compute_viscosity(self, T):
"""计算动力粘度 (Sutherland公式)"""
T0 = 273.15
mu0 = 1.716e-5
S = 110.4
mu = mu0 * (T/T0)**1.5 * (T0 + S)/(T + S)
return mu

def compute_kinematic_viscosity(self, T, rho):
"""计算运动粘度"""
mu = self.compute_viscosity(T)
nu = mu / rho
return nu

COSPA 大气模型：

class COSPAAtmosphere:
"""委员会标准大气模型 (COSPAR Reference Atmosphere)"""

def __init__(self, season='annual', latitude=45.0):
"""
参数:
season: 'annual', 'january', 'july', 'equinox', 'solstice'
latitude: 纬度 (度)
"""
self.season = season
self.latitude = latitude

# 加载COSPA数据
self.data = self.load_cospa_data()

def load_cospa_data(self):
"""加载COSPA数据表"""
# 这里简化表示，实际应从文件加载
# COSPA提供不同季节、纬度的详细大气剖面
data = {
'altitude': np.array([0, 1000, 2000, 5000, 10000, 15000,
20000, 30000, 40000, 50000, 60000]),
'temperature': np.array([288.15, 281.65, 275.15, 255.65, 223.15,
216.65, 216.65, 228.65, 270.65, 270.65,
247.02]),
'pressure': np.array([101325, 89876, 79495, 54020, 26436,
12045, 5475, 1197, 287, 80, 22]),
'density': np.array([1.225, 1.111, 1.006, 0.736, 0.412,
0.194, 0.088, 0.018, 0.004, 0.001, 0.0003])
}
return data

def compute_at_altitude(self, altitude):
"""通过插值计算"""
from scipy.interpolate import interp1d

f_T = interp1d(self.data['altitude'], self.data['temperature'],
kind='cubic', fill_value='extrapolate')
f_P = interp1d(self.data['altitude'], self.data['pressure'],
kind='cubic', fill_value='extrapolate')
f_rho = interp1d(self.data['altitude'], self.data['density'],
kind='cubic', fill_value='extrapolate')

T = float(f_T(altitude))
P = float(f_P(altitude))
rho = float(f_rho(altitude))
a = np.sqrt(1.4 * 287.05 * T)

return {
'temperature': T,
'pressure': P,
'density': rho,
'speed_of_sound': a
}

2.2 真实大气数据集成

class RealAtmosphere:
"""真实大气数据集成"""

def __init__(self, source='NOAA', date=None, location=None):
self.source = source
self.date = date or datetime.now()
self.location = location # (latitude, longitude)

# 数据源映射
self.sources = {
'NOAA': self.load_noaa_data,
'ECMWF': self.load_ecmwf_data,
'GFS': self.load_gfs_data,
'MERRA2': self.load_merra2_data,
'NRLMSISE': self.load_nrlmsise_data
}

# 加载数据
self.data = self.load_data()

def load_noaa_data(self):
"""加载NOAA数据"""
import requests

# NOAA API端点
lat, lon = self.location
date_str = self.date.strftime('%Y-%m-%d')

# 获取探空数据
url = f"https://www.ncei.noaa.gov/access/services/data/v1"
params = {
'dataset': 'global-gridded-sounding',
'dataTypes': 'TEMP,PRES,HGHT',
'stations': f'{lat},{lon}',
'startDate': date_str,
'endDate': date_str,
'format': 'json'
}

response = requests.get(url, params=params)
data = response.json()

# 解析数据
altitudes = []
temperatures = []
pressures = []
densities = []

for record in data:
alt = float(record['height'])
temp = float(record['temperature'])
pres = float(record['pressure'])

# 计算密度: ρ = P / (R * T)
R = 287.05
density = pres / (R * temp)

altitudes.append(alt)
temperatures.append(temp)
pressures.append(pres)
densities.append(density)

return {
'altitude': np.array(altitudes),
'temperature': np.array(temperatures),
'pressure': np.array(pressures),
'density': np.array(densities)
}

def load_nrlmsise_data(self):
"""加载NRLMSISE-00模型数据"""
# NRLMSISE-00是NASA的标准大气模型，考虑太阳和地磁活动

from datetime import datetime
import spaceweather

# 获取空间天气数据
sw_data = spaceweather.get_spaceweather_data(self.date)

# NRLMSISE-00输入参数
# 这里简化表示，实际实现需要调用NRLMSISE-00 Fortran代码
f107 = sw_data['f107'] # 10.7cm太阳通量
f107a = sw_data['f107a'] # 81天平均f107
ap = sw_data['ap'] # 地磁指数

# 计算大气参数
# 实际实现会调用NRLMSISE-00模型
altitudes = np.linspace(0, 100000, 1000)
temperatures = np.zeros_like(altitudes)
densities = np.zeros_like(altitudes)

for i, alt in enumerate(altitudes):
# 简化计算，实际应使用完整模型
T, rho = self.nrlmsise00_compute(alt, f107, f107a, ap)
temperatures[i] = T
densities[i] = rho

return {
'altitude': altitudes,
'temperature': temperatures,
'density': densities
}

def get_profile(self, altitude, parameter='all'):
"""获取大气剖面"""
from scipy.interpolate import interp1d

if parameter == 'all':
f_T = interp1d(self.data['altitude'], self.data['temperature'],
kind='cubic', fill_value='extrapolate')
f_P = interp1d(self.data['altitude'], self.data['pressure'],
kind='cubic', fill_value='extrapolate')
f_rho = interp1d(self.data['altitude'], self.data['density'],
kind='cubic', fill_value='extrapolate')

T = float(f_T(altitude))
P = float(f_P(altitude))
rho = float(f_rho(altitude))
a = np.sqrt(1.4 * 287.05 * T)

return {
'temperature': T,
'pressure': P,
'density': rho,
'speed_of_sound': a
}
else:
f_param = interp1d(self.data['altitude'], self.data[parameter],
kind='cubic', fill_value='extrapolate')
return float(f_param(altitude))

三、风场模型详解

3.1 基本风场模型

class WindModel:
"""风场模型基类"""

def __init__(self, reference_height=10.0):
self.reference_height = reference_height # 参考高度 (m)

def get_wind_vector(self, altitude, time=None):
"""获取风矢量 (东向, 北向, 垂直) (m/s)"""
raise NotImplementedError

def get_wind_speed(self, altitude, time=None):
"""获取风速"""
wind = self.get_wind_vector(altitude, time)
return np.linalg.norm(wind[:2]) # 只考虑水平风

def get_wind_direction(self, altitude, time=None):
"""获取风向 (度, 从北顺时针)"""
wind = self.get_wind_vector(altitude, time)
u, v, _ = wind

# 计算风向
direction = np.degrees(np.arctan2(u, v)) # atan2(东向, 北向)
direction = (direction + 360) % 360 # 转换为0-360度

return direction

3.2 恒定风场

class ConstantWind(WindModel):
"""恒定风场模型"""

def __init__(self, wind_speed=0.0, wind_direction=0.0,
reference_height=10.0):
super().__init__(reference_height)

# 风向转换为弧度
direction_rad = np.radians(wind_direction)

# 计算风矢量分量 (东向, 北向)
# 风向定义: 0度=北风, 90度=东风
self.u = -wind_speed * np.sin(direction_rad) # 东向分量
self.v = -wind_speed * np.cos(direction_rad) # 北向分量

def get_wind_vector(self, altitude, time=None):
"""恒定风，不随高度变化"""
return np.array([self.u, self.v, 0.0])

3.3 风切变模型

class WindShear(WindModel):
"""风切变模型"""

def __init__(self, wind_speed_10m=5.0, wind_direction=0.0,
shear_exponent=0.14, reference_height=10.0):
"""
参数:
wind_speed_10m: 10米高度风速 (m/s)
wind_direction: 风向 (度)
shear_exponent: 风切变指数
reference_height: 参考高度 (m)
"""
super().__init__(reference_height)

self.wind_speed_ref = wind_speed_10m
self.direction = wind_direction
self.shear_exponent = shear_exponent

# 计算参考高度处的风矢量
direction_rad = np.radians(wind_direction)
self.u_ref = -wind_speed_10m * np.sin(direction_rad)
self.v_ref = -wind_speed_10m * np.cos(direction_rad)

def get_wind_vector(self, altitude, time=None):
"""使用指数律计算风切变"""
if altitude <= 0:
# 地面风速为0
return np.array([0.0, 0.0, 0.0])

# 指数律: U(z) = U_ref * (z/z_ref)^alpha
z = max(altitude, 0.1) # 避免除零
wind_factor = (z / self.reference_height) ** self.shear_exponent

u = self.u_ref * wind_factor
v = self.v_ref * wind_factor

return np.array([u, v, 0.0])

def get_logarithmic_profile(self, altitude, roughness_length=0.03):
"""使用对数律计算风切变"""
if altitude <= roughness_length:
return np.array([0.0, 0.0, 0.0])

# 对数律: U(z) = (u*/κ) * ln(z/z0)
# 其中u*是摩擦速度，κ=0.4是冯卡门常数

# 计算摩擦速度
u_star = self.wind_speed_ref * 0.4 / np.log(self.reference_height/roughness_length)

# 计算风速
wind_speed = (u_star / 0.4) * np.log(max(altitude, roughness_length)/roughness_length)

# 转换为矢量
direction_rad = np.radians(self.direction)
u = -wind_speed * np.sin(direction_rad)
v = -wind_speed * np.cos(direction_rad)

return np.array([u, v, 0.0])

3.4 阵风模型

class WindGust(WindModel):
"""阵风模型"""

def __init__(self, base_wind_speed=5.0, base_direction=0.0,
gust_amplitude=10.0, gust_duration=10.0,
gust_start_time=30.0, reference_height=10.0):
super().__init__(reference_height)

self.base_wind_speed = base_wind_speed
self.base_direction = base_direction
self.gust_amplitude = gust_amplitude
self.gust_duration = gust_duration
self.gust_start_time = gust_start_time

# 基流风矢量
direction_rad = np.radians(base_direction)
self.u_base = -base_wind_speed * np.sin(direction_rad)
self.v_base = -base_wind_speed * np.cos(direction_rad)

def get_wind_vector(self, altitude, time=None):
"""计算包含阵风的风矢量"""
if time is None:
return np.array([self.u_base, self.v_base, 0.0])

# 计算阵风因子
gust_factor = self.compute_gust_factor(time)

# 总风速
wind_speed = self.base_wind_speed + self.gust_amplitude * gust_factor

# 风向可能有轻微变化（可选）
direction = self.base_direction

direction_rad = np.radians(direction)
u = -wind_speed * np.sin(direction_rad)
v = -wind_speed * np.cos(direction_rad)

return np.array([u, v, 0.0])

def compute_gust_factor(self, time):
"""计算阵风因子 (0到1)"""
t_rel = time – self.gust_start_time

if t_rel < 0 or t_rel > self.gust_duration:
return 0.0

# 使用余弦形状的阵风
return 0.5 * (1 – np.cos(2 * np.pi * t_rel / self.gust_duration))

def get_1_minus_cosine_gust(self, time, amplitude, duration, start_time):
"""1-cosine阵风模型 (EASA标准)"""
t_rel = time – start_time

if t_rel < 0 or t_rel > duration:
return 0.0

return 0.5 * amplitude * (1 – np.cos(2 * np.pi * t_rel / duration))

3.5 湍流模型

class TurbulenceModel(WindModel):
"""大气湍流模型"""

def __init__(self, mean_wind_speed=5.0, mean_direction=0.0,
turbulence_intensity=0.1, scale_length=100.0,
reference_height=10.0, seed=None):
"""
参数:
turbulence_intensity: 湍流强度 (σ_u / U_mean)
scale_length: 湍流尺度 (m)
"""
super().__init__(reference_height)

self.mean_wind_speed = mean_wind_speed
self.mean_direction = mean_direction
self.turbulence_intensity = turbulence_intensity
self.scale_length = scale_length

# 随机种子
self.rng = np.random.RandomState(seed)

# 生成湍流序列
self.turbulence_cache = {}

def get_wind_vector(self, altitude, time=None):
"""获取包含湍流的风矢量"""

# 平均风
direction_rad = np.radians(self.mean_direction)
u_mean = -self.mean_wind_speed * np.sin(direction_rad)
v_mean = -self.mean_wind_speed * np.cos(direction_rad)

if time is None:
return np.array([u_mean, v_mean, 0.0])

# 添加湍流
u_turb, v_turb = self.compute_turbulence(altitude, time)

u_total = u_mean + u_turb
v_total = v_mean + v_turb

return np.array([u_total, v_total, 0.0])

def compute_turbulence(self, altitude, time):
"""计算湍流分量 (Dryden模型)"""

# 湍流标准差
sigma_u = self.mean_wind_speed * self.turbulence_intensity

# 湍流尺度
L_u = self.scale_length

# Dryden谱滤波器参数
# 传递函数: H(s) = σ_u * √(2L_u/πU) * 1/(1 + L_u*s/U)

# 时间常数
tau = L_u / self.mean_wind_speed

# 生成白噪声
if time not in self.turbulence_cache:
# 使用一阶滤波器模拟湍流
dt = 0.1 # 时间步长
n_steps = int(time / dt) + 1

# 生成白噪声序列
white_noise = self.rng.randn(n_steps) * np.sqrt(dt)

# 滤波器响应
turbulence = np.zeros(n_steps)
for i in range(1, n_steps):
# 一阶低通滤波器
turbulence[i] = (1 – dt/tau) * turbulence[i-1] + \\
sigma_u * np.sqrt(2*dt/tau) * white_noise[i]

self.turbulence_cache[time] = turbulence[-1]

turbulence_value = self.turbulence_cache[time]

# 纵向和横向湍流 (简化假设相同)
u_turb = turbulence_value
v_turb = turbulence_value * 0.8 # 横向湍流通常较小

return u_turb, v_turb

def get_von_karman_turbulence(self, frequency):
"""Von Karman湍流谱"""
# Von Karman谱密度
L = self.scale_length
sigma = self.mean_wind_speed * self.turbulence_intensity

# 纵向谱
f = frequency
phi_u = (sigma**2 * 2*L) / (np.pi * (1 + (1.339*L*f)**2)**(5/6))

# 横向谱
phi_v = (sigma**2 * 2*L) / (np.pi * (1 + (1.339*L*f)**2)**(5/6))

return phi_u, phi_v

四、真实风场数据集成

4.1 NOAA GFS 数据

class NOAA_GFS_Wind:
"""NOAA GFS (全球预报系统) 风场数据"""

def __init__(self, date, location, forecast_hour=0):
self.date = date
self.location = location # (lat, lon)
self.forecast_hour = forecast_hour

# GFS参数
self.levels = [
1000, 975, 950, 925, 900, 875, 850, 825, 800, 775, 750,
725, 700, 675, 650, 625, 600, 575, 550, 525, 500, 475,
450, 425, 400, 375, 350, 325, 300, 275, 250, 225, 200,
175, 150, 125, 100, 70, 50, 30, 20, 10, 7, 5, 3, 2, 1
] # hPa

def download_gfs_data(self):
"""下载GFS数据"""
import xarray as xr
import cfgrib

lat, lon = self.location
date_str = self.date.strftime('%Y%m%d')

# GFS数据URL模板
base_url = "https://nomads.ncep.noaa.gov/dods/gfs_0p25"
url = f"{base_url}/gfs{date_str}/gfs_0p25_{self.forecast_hour:02d}z"

try:
# 打开数据集
ds = xr.open_dataset(url, engine='pydap')

# 提取风场数据
u_wind = ds['ugrdprs'] # 东向风
v_wind = ds['vgrdprs'] # 北向风

# 插值到指定位置
u_interp = u_wind.interp(lat=lat, lon=lon, method='linear')
v_interp = v_wind.interp(lat=lat, lon=lon, method='linear')

return {
'u_wind': u_interp.values,
'v_wind': v_interp.values,
'levels': self.levels,
'lat': lat,
'lon': lon,
'time': self.date
}

except Exception as e:
print(f"下载GFS数据失败: {e}")
return None

def get_wind_profile(self, altitude):
"""获取指定高度的风剖面"""
from scipy.interpolate import interp1d

# 将气压高度转换为几何高度 (简化)
# 实际应该使用大气模型进行转换
geometric_heights = self.pressure_to_height(self.levels)

# 插值
f_u = interp1d(geometric_heights, self.data['u_wind'],
kind='cubic', fill_value='extrapolate')
f_v = interp1d(geometric_heights, self.data['v_wind'],
kind='cubic', fill_value='extrapolate')

u = float(f_u(altitude))
v = float(f_v(altitude))

return np.array([u, v, 0.0])

def pressure_to_height(self, pressure_hpa):
"""气压高度转换为几何高度 (简化)"""
# 使用标准大气公式
heights = []
for p in pressure_hpa:
# 压高公式: h = 44330 * (1 – (P/P0)^(1/5.255))
P0 = 1013.25 # 海平面标准气压 (hPa)
h = 44330 * (1 – (p/P0)**(1/5.255))
heights.append(h)

return np.array(heights)

4.2 ECMWF ERA5 数据

class ECMWF_ERA5_Wind:
"""ECMWF ERA5 再分析数据"""

def __init__(self, date, location, api_key=None):
self.date = date
self.location = location
self.api_key = api_key

# ERA5压力层
self.levels = [
1, 2, 3, 5, 7, 10, 20, 30, 50, 70, 100, 125, 150, 175, 200,
225, 250, 300, 350, 400, 450, 500, 550, 600, 650, 700, 750,
775, 800, 825, 850, 875, 900, 925, 950, 975, 1000
] # hPa

def download_era5_data(self):
"""通过CDS API下载ERA5数据"""
import cdsapi

if not self.api_key:
raise ValueError("需要ECMWF CDS API密钥")

c = cdsapi.Client(key=self.api_key)

lat, lon = self.location
date_str = self.date.strftime('%Y-%m-%d')

request = {
'product_type': 'reanalysis',
'format': 'netcdf',
'variable': ['u_component_of_wind', 'v_component_of_wind'],
'pressure_level': self.levels,
'year': self.date.year,
'month': self.date.month,
'day': self.date.day,
'time': '00:00',
'area': [lat+1, lon-1, lat-1, lon+1], # 小区域
}

# 下载数据
filename = f'era5_wind_{date_str}.nc'
c.retrieve('reanalysis-era5-pressure-levels', request, filename)

# 加载数据
import xarray as xr
ds = xr.open_dataset(filename)

return ds

五、综合环境类

5.1 Environment 主类

class Environment:
"""综合环境类"""

def __init__(self, latitude=0.0, longitude=0.0, elevation=0.0,
date=None, atmosphere_model='standard',
wind_model='constant', earth_model='WGS84'):

# 地理位置
self.latitude = latitude
self.longitude = longitude
self.elevation = elevation

# 日期时间
self.date = date or datetime.now()

# 地球模型
self.earth_model = self.create_earth_model(earth_model)

# 大气模型
self.atmosphere = self.create_atmosphere_model(atmosphere_model)

# 风场模型
self.wind = self.create_wind_model(wind_model)

# 重力模型
self.gravity_model = self.create_gravity_model()

# 缓存
self.cache = {}

def create_atmosphere_model(self, model_type):
"""创建大气模型"""
if model_type == 'standard':
return USStandardAtmosphere1976()
elif model_type == 'cospar':
return COSPAAtmosphere()
elif model_type == 'custom':
return CustomAtmosphere()
elif model_type == 'real':
return RealAtmosphere(source='NOAA', date=self.date,
location=(self.latitude, self.longitude))
else:
raise ValueError(f"未知的大气模型: {model_type}")

def create_wind_model(self, model_type):
"""创建风场模型"""
if model_type == 'constant':
return ConstantWind(wind_speed=0.0, wind_direction=0.0)
elif model_type == 'shear':
return WindShear(wind_speed_10m=5.0, wind_direction=0.0)
elif model_type == 'gust':
return WindGust(base_wind_speed=5.0, base_direction=0.0)
elif model_type == 'turbulence':
return TurbulenceModel(mean_wind_speed=5.0, mean_direction=0.0)
elif model_type == 'real':
return RealWind(source='GFS', date=self.date,
location=(self.latitude, self.longitude))
else:
raise ValueError(f"未知的风场模型: {model_type}")

def create_earth_model(self, model_type):
"""创建地球模型"""
if model_type == 'WGS84':
return WGS84EarthModel()
else:
raise ValueError(f"未知的地球模型: {model_type}")

def create_gravity_model(self):
"""创建重力模型"""
return GravityModel(self.latitude, self.elevation)

def get_conditions(self, altitude, time=None):
"""获取指定高度和时间的综合环境条件"""

cache_key = (altitude, time)
if cache_key in self.cache:
return self.cache[cache_key]

# 大气条件
atmosphere = self.atmosphere.compute_at_altitude(altitude)

# 风场条件
wind_vector = self.wind.get_wind_vector(altitude, time)

# 重力
gravity = self.gravity_model.get_gravity(altitude)

# 科里奥利参数
coriolis = self.earth_model.coriolis_parameter(self.latitude)

conditions = {
'altitude': altitude,
'time': time,
'temperature': atmosphere['temperature'], # K
'pressure': atmosphere['pressure'], # Pa
'density': atmosphere['density'], # kg/m³
'speed_of_sound': atmosphere['speed_of_sound'], # m/s
'wind_u': wind_vector[0], # 东向风 (m/s)
'wind_v': wind_vector[1], # 北向风 (m/s)
'wind_w': wind_vector[2], # 垂直风 (m/s)
'wind_speed': np.linalg.norm(wind_vector[:2]),
'wind_direction': self.compute_wind_direction(wind_vector[0], wind_vector[1]),
'gravity': gravity, # m/s²
'coriolis_parameter': coriolis, # 1/s
'dynamic_pressure': 0.5 * atmosphere['density'] *
np.linalg.norm(wind_vector[:2])**2
}

self.cache[cache_key] = conditions
return conditions

def compute_wind_direction(self, u, v):
"""计算风向 (度，从北顺时针)"""
if u == 0 and v == 0:
return 0.0

# 风向: atan2(东向, 北向)
direction_rad = np.arctan2(u, v)
direction_deg = np.degrees(direction_rad)
direction_deg = (direction_deg + 360) % 360

return direction_deg

def get_density(self, altitude):
"""获取密度"""
return self.atmosphere.compute_at_altitude(altitude)['density']

def get_pressure(self, altitude):
"""获取压力"""
return self.atmosphere.compute_at_altitude(altitude)['pressure']

def get_temperature(self, altitude):
"""获取温度"""
return self.atmosphere.compute_at_altitude(altitude)['temperature']

def get_speed_of_sound(self, altitude):
"""获取声速"""
return self.atmosphere.compute_at_altitude(altitude)['speed_of_sound']

def get_wind_components(self, altitude, time=None):
"""获取风分量"""
wind = self.wind.get_wind_vector(altitude, time)
return {
'u': wind[0], # 东向
'v': wind[1], # 北向
'w': wind[2] # 垂直
}

def set_date(self, year, month, day, hour=0, minute=0, second=0):
"""设置日期时间"""
self.date = datetime(year, month, day, hour, minute, second)
self.clear_cache()

def set_location(self, latitude, longitude, elevation=0.0):
"""设置地理位置"""
self.latitude = latitude
self.longitude = longitude
self.elevation = elevation

# 更新重力模型
self.gravity_model = self.create_gravity_model()
self.clear_cache()

def set_atmosphere_model(self, model_type, **kwargs):
"""设置大气模型"""
if model_type == 'standard':
self.atmosphere = USStandardAtmosphere1976()
elif model_type == 'custom':
self.atmosphere = CustomAtmosphere(**kwargs)
elif model_type == 'real':
self.atmosphere = RealAtmosphere(**kwargs)

self.clear_cache()

def set_wind_model(self, model_type, **kwargs):
"""设置风场模型"""
if model_type == 'constant':
self.wind = ConstantWind(**kwargs)
elif model_type == 'shear':
self.wind = WindShear(**kwargs)
elif model_type == 'gust':
self.wind = WindGust(**kwargs)
elif model_type == 'turbulence':
self.wind = TurbulenceModel(**kwargs)
elif model_type == 'real':
self.wind = RealWind(**kwargs)

self.clear_cache()

def clear_cache(self):
"""清除缓存"""
self.cache.clear()

六、高级功能

6.1 风场可视化

class WindVisualizer:
"""风场可视化工具"""

def __init__(self, environment):
self.env = environment

def plot_wind_profile(self, max_altitude=30000, n_points=100):
"""绘制风剖面图"""
import matplotlib.pyplot as plt

altitudes = np.linspace(0, max_altitude, n_points)
wind_speeds = []
wind_directions = []

for alt in altitudes:
conditions = self.env.get_conditions(alt)
wind_speeds.append(conditions['wind_speed'])
wind_directions.append(conditions['wind_direction'])

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))

# 风速剖面
ax1.plot(wind_speeds, altitudes)
ax1.set_xlabel('风速 (m/s)')
ax1.set_ylabel('高度 (m)')
ax1.set_title('风速剖面')
ax1.grid(True)

# 风向剖面
ax2.plot(wind_directions, altitudes)
ax2.set_xlabel('风向 (度)')
ax2.set_ylabel('高度 (m)')
ax2.set_title('风向剖面')
ax2.grid(True)
ax2.set_xlim(0, 360)

plt.tight_layout()
return fig

def plot_wind_rose(self, altitude=1000, n_bins=16):
"""绘制风玫瑰图"""
import matplotlib.pyplot as plt
from windrose import WindroseAxes

# 获取风数据
conditions = self.env.get_conditions(altitude)
wind_speed = conditions['wind_speed']
wind_direction = conditions['wind_direction']

fig = plt.figure(figsize=(8, 8))
ax = WindroseAxes.from_ax(fig=fig)

# 绘制风玫瑰
ax.bar(wind_direction, wind_speed, normed=True,
opening=0.8, edgecolor='white', bins=np.arange(0, 21, 5))
ax.set_legend(title='风速 (m/s)')
ax.set_title(f'高度 {altitude}m 风玫瑰图')

return fig

def plot_3d_wind_field(self, x_range, y_range, z_range, resolution=20):
"""绘制3D风场"""
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

x = np.linspace(x_range[0], x_range[1], resolution)
y = np.linspace(y_range[0], y_range[1], resolution)
z = np.linspace(z_range[0], z_range[1], resolution)

X, Y, Z = np.meshgrid(x, y, z, indexing='ij')

U = np.zeros_like(X)
V = np.zeros_like(Y)
W = np.zeros_like(Z)

for i in range(resolution):
for j in range(resolution):
for k in range(resolution):
# 计算高度
alt = Z[i, j, k]
wind = self.env.get_wind_components(alt)
U[i, j, k] = wind['u']
V[i, j, k] = wind['v']
W[i, j, k] = wind['w']

fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')

# 绘制流线
ax.streamplot(X[:, :, resolution//2], Y[:, :, resolution//2],
U[:, :, resolution//2], V[:, :, resolution//2],
color='blue', linewidth=0.5, arrowsize=0.5)

ax.set_xlabel('东向 (m)')
ax.set_ylabel('北向 (m)')
ax.set_zlabel('高度 (m)')
ax.set_title('3D风场')

return fig

6.2 环境条件对弹道的影响分析

class EnvironmentalSensitivity:
"""环境条件敏感性分析"""

def __init__(self, rocket, environment):
self.rocket = rocket
self.env = environment

def analyze_wind_sensitivity(self, wind_speeds, wind_directions):
"""分析风敏感性"""
results = []

for ws in wind_speeds:
for wd in wind_directions:
# 设置风场
self.env.set_wind_model('constant',
wind_speed=ws,
wind_direction=wd)

# 运行仿真
flight = Flight(rocket=self.rocket, environment=self.env)
flight.run()

results.append({
'wind_speed': ws,
'wind_direction': wd,
'apogee': flight.apogee,
'impact_range': flight.x_impact,
'drift': flight.y_impact
})

return pd.DataFrame(results)

def analyze_density_sensitivity(self, density_factors):
"""分析密度敏感性"""
results = []

for factor in density_factors:
# 创建自定义大气模型
class ModifiedAtmosphere(USStandardAtmosphere1976):
def compute_at_altitude(self, altitude):
base = super().compute_at_altitude(altitude)
base['density'] *= factor
# 保持压力温度关系
base['pressure'] = base['density'] * 287.05 * base['temperature']
return base

self.env.atmosphere = ModifiedAtmosphere()

# 运行仿真
flight = Flight(rocket=self.rocket, environment=self.env)
flight.run()

results.append({
'density_factor': factor,
'apogee': flight.apogee,
'max_velocity': flight.max_speed,
'impact_range': flight.x_impact
})

return pd.DataFrame(results)

def monte_carlo_environment(self, n_simulations=1000):
"""环境条件蒙特卡洛分析"""
results = []

for i in range(n_simulations):
# 随机采样环境参数
wind_speed = np.random.normal(5, 2) # 平均5m/s，标准差2m/s
wind_direction = np.random.uniform(0, 360)
density_factor = np.random.normal(1.0, 0.05) # 密度变化±5%

# 设置风场
self.env.set_wind_model('constant',
wind_speed=max(wind_speed, 0),
wind_direction=wind_direction)

# 设置大气
class RandomAtmosphere(USStandardAtmosphere1976):
def compute_at_altitude(self, altitude):
base = super().compute_at_altitude(altitude)
base['density'] *= density_factor
base['pressure'] = base['density'] * 287.05 * base['temperature']
return base

self.env.atmosphere = RandomAtmosphere()

# 运行仿真
flight = Flight(rocket=self.rocket, environment=self.env)
flight.run()

results.append({
'simulation': i,
'wind_speed': wind_speed,
'wind_direction': wind_direction,
'density_factor': density_factor,
'apogee': flight.apogee,
'impact_range': flight.x_impact,
'drift': flight.y_impact
})

df = pd.DataFrame(results)

# 统计分析
stats = {
'apogee_mean': df['apogee'].mean(),
'apogee_std': df['apogee'].std(),
'apogee_5th': df['apogee'].quantile(0.05),
'apogee_95th': df['apogee'].quantile(0.95),
'range_mean': df['impact_range'].mean(),
'range_std': df['impact_range'].std(),
'drift_mean': df['drift'].mean(),
'drift_std': df['drift'].std()
}

return df, stats

七、使用示例

7.1 复杂环境设置

# 创建综合环境
from datetime import datetime
from rocketpy import Environment, Flight, Rocket
import numpy as np

# 1. 创建真实环境
env = Environment(
latitude=28.5721, # 肯尼迪航天中心
longitude=-80.6480,
elevation=3.0, # 海拔3米

# 设置发射时间
date=datetime(2024, 6, 15, 12, 0, 0), # 2024年6月15日中午

# 使用真实大气和风场
atmosphere_model='real',
wind_model='real',

# 地球模型
earth_model='WGS84'
)

# 2. 自定义风场叠加
from rocketpy import WindShear, WindGust, TurbulenceModel

# 添加风切变
wind_shear = WindShear(
wind_speed_10m=8.0, # 10米高度风速8m/s
wind_direction=45.0, # 东北风
shear_exponent=0.15
)

# 添加阵风
wind_gust = WindGust(
base_wind_speed=8.0,
base_direction=45.0,
gust_amplitude=12.0, # 12m/s阵风
gust_duration=15.0, # 持续15秒
gust_start_time=30.0 # 发射后30秒开始
)

# 添加湍流
turbulence = TurbulenceModel(
mean_wind_speed=8.0,
mean_direction=45.0,
turbulence_intensity=0.12, # 12%湍流强度
scale_length=150.0, # 湍流尺度150m
seed=42 # 随机种子
)

# 组合风场
class CombinedWind:
def __init__(self, *wind_models):
self.wind_models = wind_models

def get_wind_vector(self, altitude, time=None):
total_wind = np.zeros(3)
for model in self.wind_models:
total_wind += model.get_wind_vector(altitude, time)
return total_wind

combined_wind = CombinedWind(wind_shear, wind_gust, turbulence)
env.wind = combined_wind

# 3. 创建火箭并运行仿真
calisto = Rocket(
radius=0.0635,
mass=14.426,
inertia=(6.321, 6.321, 0.034)
)

# … 添加发动机、气动面等组件 …

# 4. 运行飞行仿真
flight = Flight(
rocket=calisto,
environment=env,
rail_length=5.0,
inclination=85,
heading=90
)

flight.run()

# 5. 分析环境影响
print("环境条件分析:")
print(f"发射场温度: {env.get_temperature(3):.1f} K")
print(f"发射场风速: {env.get_conditions(3)['wind_speed']:.1f} m/s")
print(f"发射场风向: {env.get_conditions(3)['wind_direction']:.0f}°")
print(f"10km高度风速: {env.get_conditions(10000)['wind_speed']:.1f} m/s")
print(f"顶点高度密度: {env.get_density(flight.apogee):.4f} kg/m³")

# 6. 风场影响分析
wind_drift = flight.y_impact
wind_correction_angle = np.degrees(np.arctan2(wind_drift, flight.x_impact))
print(f"\\n风致漂移: {wind_drift:.1f} m")
print(f"风修正角: {wind_correction_angle:.1f}°")

7.2 季节性环境变化分析

def seasonal_analysis(latitude, longitude, year=2024):
"""分析季节对环境的影响"""

seasons = {
'Winter': (datetime(year, 1, 15, 12, 0), 'january'),
'Spring': (datetime(year, 4, 15, 12, 0), 'april'),
'Summer': (datetime(year, 7, 15, 12, 0), 'july'),
'Fall': (datetime(year, 10, 15, 12, 0), 'october')
}

results = {}

for season, (date, month) in seasons.items():
# 创建季节环境
env = Environment(
latitude=latitude,
longitude=longitude,
date=date,
atmosphere_model='cospar', # 使用COSPA季节模型
wind_model='real'
)

# 设置COSPA季节
env.atmosphere.season = month

# 分析不同高度的条件
altitudes = [0, 5000, 10000, 20000, 30000]

season_data = {}
for alt in altitudes:
conditions = env.get_conditions(alt)
season_data[alt] = {
'temperature': conditions['temperature'],
'density': conditions['density'],
'wind_speed': conditions['wind_speed']
}

results[season] = season_data

return results

# 运行季节分析
lat, lon = 40.0, -100.0 # 美国中部
seasonal_results = seasonal_analysis(lat, lon, 2024)

# 打印结果
for season, data in seasonal_results.items():
print(f"\\n{season}:")
for alt in [0, 10000, 30000]:
print(f" {alt/1000:.0f}km: T={data[alt]['temperature']:.1f}K, "
f"ρ={data[alt]['density']:.4f}kg/m³, "
f"W={data[alt]['wind_speed']:.1f}m/s")

八、验证与测试

8.1 大气模型验证

def validate_atmosphere_models():
"""验证不同大气模型的一致性"""

altitudes = np.linspace(0, 80000, 100)

models = {
'US76': USStandardAtmosphere1976(),
'COSPA': COSPAAtmosphere(season='annual', latitude=45.0),
'NRLMSISE': NRLMSISE00Atmosphere()
}

results = {}

for name, model in models.items():
temperatures = []
densities = []

for alt in altitudes:
atm = model.compute_at_altitude(alt)
temperatures.append(atm['temperature'])
densities.append(atm['density'])

results[name] = {
'altitude': altitudes,
'temperature': np.array(temperatures),
'density': np.array(densities)
}

# 计算差异
print("大气模型比较 (相对US76):")
for name in ['COSPA', 'NRLMSISE']:
t_diff = np.mean(np.abs(results[name]['temperature'] –
results['US76']['temperature']))
rho_diff = np.mean(np.abs(results[name]['density'] –
results['US76']['density']))

print(f"{name}: 平均温度差 {t_diff:.2f}K, 平均密度差 {rho_diff*100:.2f}%")

return results

8.2 风场模型验证

def validate_wind_models(measured_data):
"""验证风场模型与实际测量数据的吻合度"""

models = {
'Constant': ConstantWind(wind_speed=5.0, wind_direction=0.0),
'Shear': WindShear(wind_speed_10m=5.0, wind_direction=0.0),
'Real': RealWind(source='GFS', date=datetime.now())
}

errors = {}

for name, model in models.items():
predicted = []
measured = []

for data_point in measured_data:
alt = data_point['altitude']
time = data_point.get('time')

# 模型预测
wind_pred = model.get_wind_vector(alt, time)
predicted.append(wind_pred[:2]) # 只考虑水平分量

# 实测数据
wind_meas = np.array([data_point['u'], data_point['v']])
measured.append(wind_meas)

predicted = np.array(predicted)
measured = np.array(measured)

# 计算误差
u_error = np.mean(np.abs(predicted[:,0] – measured[:,0]))
v_error = np.mean(np.abs(predicted[:,1] – measured[:,1]))
speed_error = np.mean(np.abs(np.linalg.norm(predicted, axis=1) –
np.linalg.norm(measured, axis=1)))

errors[name] = {
'u_error': u_error,
'v_error': v_error,
'speed_error': speed_error,
'rmse': np.sqrt(np.mean((predicted – measured)**2))
}

return errors

总结

大气与风场模型模块提供了：

该模块为火箭发射提供高精度的环境条件模拟，显著提高弹道预测的准确性。无论是业余火箭爱好者还是专业航天工程师，都能通过这个模块获得可靠的仿真环境。