# -*- coding: utf-8 -*-
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import interpolate
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error, r2_score
import os
import glob
import chardet
import warnings
warnings.filterwarnings("ignore")
# 绘图设置
plt.style.use('seaborn-v0_8-whitegrid')
plt.rcParams['font.family'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
# ----------------- 数据解析类 -----------------
class WindProfilerData:
def __init__(self, file_path):
self.file_path = file_path
self.obs_data = {}
def parse_float_with_missing(self, value_str, default=np.nan):
if value_str is None or '////' in str(value_str) or str(value_str).strip() == '':
return default
try:
return float(value_str)
except (ValueError, TypeError):
return default
def parse_radial_file(self):
"""解析径向格式（尽量兼容）"""
try:
with open(self.file_path, 'r', encoding='gbk', errors='ignore') as f:
content = f.readlines()
except:
try:
with open(self.file_path, 'r', encoding='utf-8', errors='ignore') as f:
content = f.readlines()
except:
with open(self.file_path, 'r', encoding='gb2312', errors='ignore') as f:
content = f.readlines()
mode = None
current_beam = None
beam_count = 0
for line in content:
line = line.strip()
if '低模式' in line or 'low' in line.lower():
mode = 'low'
elif '中模式' in line or 'medium' in line.lower():
mode = 'medium'
elif '高模式' in line or 'high' in line.lower():
mode = 'high'
if line.startswith('RAD'):
parts = line.split()
if len(parts) >= 2:
beam_name = parts[1].lower()
current_beam = f"{mode}_{beam_name}"
self.obs_data[current_beam] = []
beam_count += 1
elif line and any(c.isdigit() or c in ['-', '.', '/'] for c in line[:10]):
if current_beam and len(line.split()) >= 4:
data = line.split()
height = self.parse_float_with_missing(data[0])
spectrum_width = self.parse_float_with_missing(data[1])
snr = self.parse_float_with_missing(data[2])
radial_velocity = self.parse_float_with_missing(data[3])
if not np.isnan(height):
self.obs_data[current_beam].append({
'height': height,
'spectrum_width': spectrum_width,
'snr': snr,
'radial_velocity': radial_velocity,
'beam_name': current_beam
})
return self.obs_data
def parse_product_file(self):
"""解析产品（表格）格式"""
try:
with open(self.file_path, 'r', encoding='gbk', errors='ignore') as f:
content = f.readlines()
except:
try:
with open(self.file_path, 'r', encoding='utf-8', errors='ignore') as f:
content = f.readlines()
except:
with open(self.file_path, 'r', encoding='gb2312', errors='ignore') as f:
content = f.readlines()
product_data = []
in_data_section = False
for line in content:
line = line.strip()
if line == 'ROBS':
in_data_section = True
continue
elif line == 'NNNN':
in_data_section = False
break
if in_data_section and len(line.split()) >= 7:
data = line.split()
try:
cn2_str = data[6]
cn2 = self.parse_float_with_missing(cn2_str)
product_data.append({
'height': self.parse_float_with_missing(data[0]),
'wind_direction': self.parse_float_with_missing(data[1]),
'wind_speed': self.parse_float_with_missing(data[2]),
'vertical_velocity': self.parse_float_with_missing(data[3]),
'horizontal_confidence': self.parse_float_with_missing(data[4]),
'vertical_confidence': self.parse_float_with_missing(data[5]),
'cn2': cn2
})
except (ValueError, IndexError):
continue
return product_data
# ----------------- 微波辐射计解析类 -----------------
class MicrowaveRadiometerData:
def __init__(self, file_path):
self.file_path = file_path
self.data = {}
def detect_encoding(self):
try:
with open(self.file_path, 'rb') as f:
raw_data = f.read(10000)
result = chardet.detect(raw_data)
return result['encoding']
except:
return 'gbk'
def parse_data(self):
try:
encoding = self.detect_encoding()
with open(self.file_path, 'r', encoding=encoding, errors='ignore') as f:
first_lines = [f.readline() for _ in range(5)]
# 简单读取 csv/tsv
try:
df = pd.read_csv(self.file_path, encoding=encoding, skiprows=0)
except:
df = pd.read_csv(self.file_path, encoding=encoding, header=None)
df = df.replace('-', np.nan).replace('////', np.nan).replace('', np.nan)
# 尝试提取 H1,H2... 列
heights = []
temp_df = pd.DataFrame()
# 典型比赛数据列可能直接是 H1,H2...
for c in df.columns:
if isinstance(c, str) and c.strip().upper().startswith('H'):
heights.append(c)
if heights:
temp_df = df[heights]
# 把列名 H1->高度索引序号，默认层间隔200m
numeric_heights = []
for col in temp_df.columns:
try:
idx = int(col.strip().upper().lstrip('H'))
numeric_heights.append(idx * 200)
except:
numeric_heights.append(len(numeric_heights)*200)
self.data = {'temperatures': temp_df, 'heights': numeric_heights, 'all_data': df}
else:
# fallback: produce small example vertical profile
self.data = {
'temperatures': pd.DataFrame({'H1':[20.0],'H2':[19.5],'H3':[19.0],'H4':[18.5],'H5':[18.0]}),
'heights': [1000,2000,3000,4000,5000],
'all_data': df
}
return self.data
except Exception as e:
print("解析微波辐射计出错，使用示例数据：", e)
return {
'temperatures': pd.DataFrame({'H1':[20.0],'H2':[19.5],'H3':[19.0],'H4':[18.5],'H5':[18.0]}),
'heights': [1000,2000,3000,4000,5000],
'all_data': pd.DataFrame()
}
# ----------------- 物理计算 -----------------
class AtmosphericPhysics:
@staticmethod
def calculate_potential_temperature(temp_kelvin, pressure_hpa, p0=1000.0):
kappa = 0.286
return temp_kelvin * (p0 / pressure_hpa) ** kappa
@staticmethod
def calculate_richardson_number(theta_profile, u_profile, v_profile, heights):
"""
更稳健的理查逊数计算：
Ri = (g / Theta) * dTheta/dz / ( (du/dz)^2 + (dv/dz)^2 )
返回与输入 heights 同长度的数组
"""
theta = np.array(theta_profile, dtype=float)
u = np.array(u_profile, dtype=float)
v = np.array(v_profile, dtype=float)
z = np.array(heights, dtype=float)
# mask invalid
valid = ~(np.isnan(theta) | np.isnan(u) | np.isnan(v) | np.isnan(z))
ri = np.full_like(z, np.nan, dtype=float)
if np.sum(valid) < 3:
return ri
theta_v = theta[valid]
u_v = u[valid]
v_v = v[valid]
z_v = z[valid]
# ensure monotonic heights
order = np.argsort(z_v)
z_v = z_v[order]
theta_v = theta_v[order]
u_v = u_v[order]
v_v = v_v[order]
# gradients (central differences)
dtheta_dz = np.gradient(theta_v, z_v)
du_dz = np.gradient(u_v, z_v)
dv_dz = np.gradient(v_v, z_v)
shear_sq = du_dz**2 + dv_dz**2
eps = 1e-6
shear_sq[shear_sq < eps] = eps # 防止除零
g = 9.80665
# Theta should be positive (开尔文)
theta_v[theta_v <= 0] = 1.0
ri_v = (g / theta_v) * dtheta_dz / shear_sq
# 填回原位置
ri_valid_full = np.full_like(z, np.nan, dtype=float)
ri_valid_full[valid] = ri_v[np.argsort(order)] if len(order)>1 else ri_v
# 上面处理可能复杂，简化地把 ri_v 直接按照 valid 的顺序回填：
ri[valid] = ri_v
return ri
@staticmethod
def calculate_tke(spectrum_width):
s = np.array(spectrum_width, dtype=float)
s[np.isnan(s)] = 0.0
tke = s**2
return tke
@staticmethod
def calculate_wind_shear(u_profile, v_profile, heights):
u = np.array(u_profile, dtype=float)
v = np.array(v_profile, dtype=float)
z = np.array(heights, dtype=float)
valid = ~(np.isnan(u) | np.isnan(v) | np.isnan(z))
wind_shear = np.full_like(z, np.nan, dtype=float)
if np.sum(valid) < 2:
return wind_shear
u_v = u[valid]
v_v = v[valid]
z_v = z[valid]
order = np.argsort(z_v)
z_v = z_v[order]
u_v = u_v[order]
v_v = v_v[order]
du_dz = np.gradient(u_v, z_v)
dv_dz = np.gradient(v_v, z_v)
shear = np.sqrt(du_dz**2 + dv_dz**2)
wind_shear[valid] = shear
return wind_shear
# ----------------- 辅助数据处理函数 -----------------
def fill_missing_data(data, method='linear'):
data = np.array(data, dtype=float)
n = len(data)
if n < 2:
return data
mask = ~np.isnan(data)
if np.sum(mask) == 0:
return np.zeros_like(data)
if np.sum(mask) == 1:
return np.full_like(data, data[mask][0])
idx = np.arange(n)
f = interpolate.interp1d(idx[mask], data[mask], kind=method, bounds_error=False, fill_value="extrapolate")
return f(idx)
def calculate_wind_components_from_radial(radar_df):
"""从radial或product数据估算 u,v,w,spec。返回高度数组与四个廓线。"""
if radar_df is None or radar_df.empty:
return np.array([]), np.array([]), np.array([]), np.array([]), np.array([])
# group by height (取平均)
if 'height' not in radar_df.columns:
return np.array([]), np.array([]), np.array([]), np.array([]), np.array([])
unique_heights = np.unique(radar_df['height'].dropna().astype(float))
unique_heights = np.sort(unique_heights)
u_profiles = []
v_profiles = []
w_profiles = []
spectrum_widths = []
for h in unique_heights:
block = radar_df[radar_df['height'] == h]
if block.empty:
u_profiles.append(np.nan); v_profiles.append(np.nan); w_profiles.append(np.nan); spectrum_widths.append(np.nan)
continue
# 如果有 radial_velocity 则优先用径向合成（极为简化的合成方法）
if 'radial_velocity' in block.columns and not block['radial_velocity'].isna().all():
rv_mean = np.nanmean(block['radial_velocity'])
# 简化合成为 u,v,w 比例
u = rv_mean * 0.7
v = rv_mean * 0.7
w = rv_mean * 0.2
elif 'wind_speed' in block.columns and 'wind_direction' in block.columns:
ws = np.nanmean(block['wind_speed'])
wd = np.nanmean(block['wind_direction'])
u = ws * np.sin(np.radians(wd))
v = ws * np.cos(np.radians(wd))
w = np.nanmean(block['vertical_velocity']) if 'vertical_velocity' in block.columns else 0.0
else:
u = np.nan; v = np.nan; w = np.nan
sw = np.nanmean(block['spectrum_width']) if 'spectrum_width' in block.columns else np.nan
u_profiles.append(u); v_profiles.append(v); w_profiles.append(w); spectrum_widths.append(sw)
# 填充
u_profiles = fill_missing_data(np.array(u_profiles))
v_profiles = fill_missing_data(np.array(v_profiles))
w_profiles = fill_missing_data(np.array(w_profiles))
spectrum_widths = fill_missing_data(np.array(spectrum_widths))
return unique_heights, u_profiles, v_profiles, w_profiles, spectrum_widths
def interpolate_to_common_grid(source_values, source_heights, target_heights):
source_heights = np.array(source_heights, dtype=float)
source_values = np.array(source_values, dtype=float)
target_heights = np.array(target_heights, dtype=float)
# remove NaNs
mask = ~np.isnan(source_values) & ~np.isnan(source_heights)
if np.sum(mask) < 2:
# 返回 fill with nan
return np.full_like(target_heights, np.nan, dtype=float)
try:
f = interpolate.interp1d(source_heights[mask], source_values[mask], bounds_error=False, fill_value="extrapolate")
return f(target_heights)
except Exception:
return np.full_like(target_heights, np.nan, dtype=float)
# ----------------- 湍流模型 -----------------
class TurbulenceModel:
def __init__(self):
pass
class ModelA(TurbulenceModel):
def __init__(self, alpha1=0.4, alpha2=0.3, alpha3=0.3, beta=0.5):
super().__init__()
self.alpha1 = alpha1; self.alpha2 = alpha2; self.alpha3 = alpha3; self.beta = beta
def predict(self, ri, tke, wind_shear):
valid = ~(np.isnan(ri) | np.isnan(tke) | np.isnan(wind_shear))
out = np.full_like(ri, 0.5, dtype=float)
if np.sum(valid) == 0:
return out
ri_v = ri[valid]
tke_v = tke[valid]
shear_v = wind_shear[valid]
# map Ri to suppression factor (保持 Ri 较大 => 湍流强度较小)
f_ri = np.exp(-self.beta * ri_v)
tke_norm = tke_v / (np.nanmax(tke_v) if np.nanmax(tke_v) > 0 else 1.0)
shear_norm = shear_v / (np.nanmax(shear_v) if np.nanmax(shear_v) > 0 else 1.0)
turb_v = self.alpha1 * f_ri + self.alpha2 * tke_norm + self.alpha3 * shear_norm
out[valid] = turb_v
# 限制范围 0~1
out = np.clip(out, 0.0, 1.0)
return out
class ModelB(TurbulenceModel):
def __init__(self, c1=0.6, c2=0.4):
super().__init__()
self.c1 = c1; self.c2 = c2
def predict(self, spectrum_width, wind_shear):
s = np.array(spectrum_width, dtype=float); sh = np.array(wind_shear,dtype=float)
valid = ~(np.isnan(s) | np.isnan(sh))
out = np.full_like(s, 0.5, dtype=float)
if np.sum(valid) == 0:
return out
s_v = s[valid]; sh_v = sh[valid]
s_norm = s_v / (np.nanmax(s_v) if np.nanmax(s_v)>0 else 1.0)
sh_norm = sh_v / (np.nanmax(sh_v) if np.nanmax(sh_v)>0 else 1.0)
turb_v = self.c1 * s_norm + self.c2 * sh_norm
out[valid] = np.clip(turb_v, 0.0, 1.0)
return out
class OptimizedModelB(TurbulenceModel):
def __init__(self):
self.model = RandomForestRegressor(n_estimators=20, random_state=42, max_depth=2)
self.is_fitted = False
def fit(self, X, y):
X = np.array(X); y = np.array(y)
mask = ~(np.isnan(X).any(axis=1) | np.isnan(y))
if np.sum(mask) < 6:
self.is_fitted = False
return
self.model.fit(X[mask], y[mask])
self.is_fitted = True
def predict(self, X):
X = np.array(X)
if not self.is_fitted:
# fallback to simple ModelB
mb = ModelB()
s = X[:,0] if X.shape[1]>0 else np.zeros(X.shape[0])
sh = X[:,1] if X.shape[1]>1 else np.zeros(X.shape[0])
return mb.predict(s, sh)
mask = ~np.isnan(X).any(axis=1)
out = np.full(X.shape[0], 0.5, dtype=float)
if np.sum(mask)==0:
return out
out[mask] = self.model.predict(X[mask])
return np.clip(out, 0.0, 1.0)
# ----------------- 主流程 -----------------
def main():
print("开始处理（支持 a、b 两站点）")
# === 修改下列路径为你的实际路径 ===
radar_dir_a = "./data/第一题/a站点/a站点风廓线雷达数据"
radiometer_file_a = "./data/第一题/a站点/a站点微波辐射计数据.txt"
radar_dir_b = "./data/第一题/b站点/b站点风廓线雷达资料"
radiometer_file_b = "./data/第一题/b站点/b站点微波辐射计数据.txt"
# ======================================
# 确保输出目录存在
os.makedirs('./figure', exist_ok=True)
os.makedirs('./result', exist_ok=True)
def read_radar_folder(radar_dir, max_files=10):
files = glob.glob(os.path.join(radar_dir, "*.TXT")) + glob.glob(os.path.join(radar_dir, "*.txt"))
if not files:
return pd.DataFrame()
records = []
for fp in files[:max_files]:
wp = WindProfilerData(fp)
# 优先尝试解析 product，再尝试 radial（保持兼容）
parsed = wp.parse_product_file()
if parsed:
records.extend(parsed)
else:
parsed_rad = wp.parse_radial_file()
for beam, arr in parsed_rad.items():
for row in arr:
records.append(row)
return pd.DataFrame(records)
# 读取两站点雷达
radar_df_a = read_radar_folder(radar_dir_a, max_files=20)
radar_df_b = read_radar_folder(radar_dir_b, max_files=20)
if radar_df_a.empty:
print("未找到 a 站点雷达数据，使用示例数据（a）")
heights = np.linspace(100, 5000, 50)
radar_df_a = pd.DataFrame({
'height': heights,
'wind_speed': 5 + 1.5 * np.sin(heights / 1000),
'wind_direction': 180 + 20 * np.cos(heights / 800),
'spectrum_width': 0.5 + 0.2 * np.sin(heights / 1500),
'vertical_velocity': 0.1 * np.cos(heights / 2000)
})
if radar_df_b.empty:
print("未找到 b 站点雷达数据，使用示例数据（b）")
heights = np.linspace(150, 5200, 50)
radar_df_b = pd.DataFrame({
'height': heights,
'wind_speed': 4 + 1.2 * np.sin(heights / 1100),
'wind_direction': 170 + 25 * np.cos(heights / 900),
'spectrum_width': 0.4 + 0.25 * np.sin(heights / 1400),
'vertical_velocity': 0.08 * np.cos(heights / 2100)
})
print(f"a 站点雷达记录: {len(radar_df_a)}, b 站点雷达记录: {len(radar_df_b)}")
# 读取辐射计
mwr_a = MicrowaveRadiometerData(radiometer_file_a)
mwr_b = MicrowaveRadiometerData(radiometer_file_b)
mwr_data_a = mwr_a.parse_data()
mwr_data_b = mwr_b.parse_data()
# 对两个雷达分别计算廓线
hz_a, ua, va, wa, sw_a = calculate_wind_components_from_radial(radar_df_a)
hz_b, ub, vb, wb, sw_b = calculate_wind_components_from_radial(radar_df_b)
# 构建共同高度网格（取两站点高度的并集并以等间距插值）
# 取 min/max 与合理间隔（例如 100m）
h_min = min(np.nanmin(hz_a) if len(hz_a)>0 else 100, np.nanmin(hz_b) if len(hz_b)>0 else 100)
h_max = max(np.nanmax(hz_a) if len(hz_a)>0 else 5000, np.nanmax(hz_b) if len(hz_b)>0 else 5000)
common_heights = np.linspace(h_min, h_max, 60)
# 将两站廓线插值到共同网格后取平均（若一站缺失则取另一站）
ua_i = interpolate_to_common_grid(ua, hz_a, common_heights) if len(ua)>0 else np.full_like(common_heights, np.nan)
va_i = interpolate_to_common_grid(va, hz_a, common_heights) if len(va)>0 else np.full_like(common_heights, np.nan)
sw_a_i = interpolate_to_common_grid(sw_a, hz_a, common_heights) if len(sw_a)>0 else np.full_like(common_heights, np.nan)
ub_i = interpolate_to_common_grid(ub, hz_b, common_heights) if len(ub)>0 else np.full_like(common_heights, np.nan)
vb_i = interpolate_to_common_grid(vb, hz_b, common_heights) if len(vb)>0 else np.full_like(common_heights, np.nan)
sw_b_i = interpolate_to_common_grid(sw_b, hz_b, common_heights) if len(sw_b)>0 else np.full_like(common_heights, np.nan)
# 合并：对两个站点取均值（忽略 nan），最终廓线用于后续计算
u_profiles = np.nanmean(np.vstack([ua_i, ub_i]), axis=0)
v_profiles = np.nanmean(np.vstack([va_i, vb_i]), axis=0)
spectrum_widths = np.nanmean(np.vstack([sw_a_i, sw_b_i]), axis=0)
# 若仍有 nan（两站都 nan），用小默认值填充
u_profiles = fill_missing_data(u_profiles)
v_profiles = fill_missing_data(v_profiles)
spectrum_widths = fill_missing_data(spectrum_widths)
# 处理辐射计：我们取 a、b 两站的平均温度廓线作为位温输入（优先 a，如果缺则 b）
def extract_temp_profile(mwr_data):
if 'temperatures' in mwr_data and not mwr_data['temperatures'].empty:
# 使用第一行
row = mwr_data['temperatures'].iloc[0]
heights = mwr_data['heights']
vals = []
for i in range(1, len(heights)+1):
cname = f'H{i}'
if cname in row:
vals.append(row[cname])
else:
# 如果列名不是 H1...直接尝试按列顺序
try:
vals.append(row.iloc[i-1])
except:
vals.append(np.nan)
return np.array(vals, dtype=float), np.array(heights, dtype=float)
else:
return np.array([]), np.array([])
ta_vals_a, ta_h_a = extract_temp_profile(mwr_data_a)
tb_vals_b, tb_h_b = extract_temp_profile(mwr_data_b)
# 如果两站都有资料，插值到共同高度并取均值；否则用可用站点或示例
if len(ta_vals_a)>0:
theta_a_heights = ta_h_a
temp_a = ta_vals_a
else:
theta_a_heights = np.array([1000,2000,3000,4000,5000])
temp_a = np.array([20,19.5,19,18.5,18])
if len(tb_vals_b)>0:
theta_b_heights = tb_h_b
temp_b = tb_vals_b
else:
theta_b_heights = np.array([1000,2000,3000,4000,5000])
temp_b = np.array([19,18.5,18,17.5,17])
# 插值到共同网格
temp_a_i = interpolate_to_common_grid(temp_a, theta_a_heights, common_heights)
temp_b_i = interpolate_to_common_grid(temp_b, theta_b_heights, common_heights)
# 合并温度廓线（摄氏度假设），优先 a 与 b 均值
temp_combined = np.nanmean(np.vstack([temp_a_i, temp_b_i]), axis=0)
temp_combined = fill_missing_data(temp_combined)
# 为位温计算转换为开尔文并假设气压（指数递减）
temp_k = temp_combined + 273.15
pressure_profile = 1013.25 * np.exp(-common_heights / 8400.0) # hPa approx
physics = AtmosphericPhysics()
theta_profile = physics.calculate_potential_temperature(temp_k, pressure_profile)
# 确保长度一致
heights = common_heights
min_len = min(len(heights), len(u_profiles), len(v_profiles), len(spectrum_widths), len(theta_profile))
heights = np.array(heights[:min_len])
u_profiles = np.array(u_profiles[:min_len])
v_profiles = np.array(v_profiles[:min_len])
spectrum_widths = np.array(spectrum_widths[:min_len])
theta_profile = np.array(theta_profile[:min_len])
# 计算 Ri、TKE、风切变
ri = physics.calculate_richardson_number(theta_profile, u_profiles, v_profiles, heights)
# 如果 ri 含有 nan 或极端值，用平滑/填充处理
ri = fill_missing_data(np.where(np.isfinite(ri), ri, np.nan))
# 限制 ri 合理范围（例如 -10 ~ 10）
ri = np.clip(ri, -10, 10)
tke = physics.calculate_tke(spectrum_widths)
tke = fill_missing_data(tke)
wind_shear = physics.calculate_wind_shear(u_profiles, v_profiles, heights)
wind_shear = fill_missing_data(wind_shear)
# 模型计算
model_a = ModelA()
turbulence_a = model_a.predict(ri, tke, wind_shear)
model_b = ModelB()
turbulence_b = model_b.predict(spectrum_widths, wind_shear)
# 优化模型 b
X_train = np.column_stack([spectrum_widths, wind_shear, u_profiles, v_profiles])
y_train = turbulence_a
opt_b = OptimizedModelB()
opt_b.fit(X_train, y_train)
turbulence_b_opt = opt_b.predict(X_train)
# 评估
valid_mask = ~(np.isnan(turbulence_a) | np.isnan(turbulence_b) | np.isnan(turbulence_b_opt))
if np.sum(valid_mask) > 5:
rmse_orig = np.sqrt(mean_squared_error(turbulence_a[valid_mask], turbulence_b[valid_mask]))
rmse_opt = np.sqrt(mean_squared_error(turbulence_a[valid_mask], turbulence_b_opt[valid_mask]))
r2_orig = r2_score(turbulence_a[valid_mask], turbulence_b[valid_mask])
r2_opt = r2_score(turbulence_a[valid_mask], turbulence_b_opt[valid_mask])
print(f"原始模型b RMSE: {rmse_orig:.4f}, R2: {r2_orig:.4f}")
print(f"优化模型b RMSE: {rmse_opt:.4f}, R2: {r2_opt:.4f}")
else:
print("有效数据不足以计算评估指标")
rmse_orig = rmse_opt = r2_orig = r2_opt = np.nan
# 可视化
print("绘图并保存结果...")
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
axes[0,0].plot(turbulence_a, heights, '-', linewidth=2, label='模型a (融合)')
axes[0,0].plot(turbulence_b, heights, '--', linewidth=2, label='模型b (原始)')
axes[0,0].plot(turbulence_b_opt, heights, '-.', linewidth=2, label='模型b (优化)')
axes[0,0].set_xlabel('湍流强度')
axes[0,0].set_ylabel('高度 (m)')
axes[0,0].legend(); axes[0,0].grid(True); axes[0,0].set_title('湍流强度廓线对比')
axes[0,1].plot(ri, heights, '-', linewidth=2, color='purple')
axes[0,1].axvline(x=0.25, color='r', linestyle='--', label='Ri=0.25')
axes[0,1].axvline(x=1.0, color='b', linestyle='--', label='Ri=1.0')
axes[0,1].set_xlabel('理查逊数 (Ri)'); axes[0,1].set_ylabel('高度 (m)')
axes[0,1].legend(); axes[0,1].grid(True); axes[0,1].set_title('理查逊数廓线')
axes[0,2].plot(tke, heights, '-', linewidth=2, color='orange')
axes[0,2].set_xlabel('湍流动能 (TKE)'); axes[0,2].set_ylabel('高度 (m)'); axes[0,2].grid(True)
axes[0,2].set_title('湍流动能廓线')
axes[1,0].plot(wind_shear, heights, '-', linewidth=2)
axes[1,0].set_xlabel('风切变 (s^-1)'); axes[1,0].set_ylabel('高度 (m)'); axes[1,0].grid(True)
axes[1,0].set_title('风切变廓线')
axes[1,1].plot(spectrum_widths, heights, '-', linewidth=2)
axes[1,1].set_xlabel('谱宽 (m/s)'); axes[1,1].set_ylabel('高度 (m)'); axes[1,1].grid(True)
axes[1,1].set_title('谱宽廓线')
if np.sum(valid_mask) > 5:
axes[1,2].scatter(turbulence_a[valid_mask], turbulence_b[valid_mask], alpha=0.6, label='原始模型b')
axes[1,2].scatter(turbulence_a[valid_mask], turbulence_b_opt[valid_mask], alpha=0.6, label='优化模型b')
mval = max(np.nanmax(turbulence_a[valid_mask]), np.nanmax(turbulence_b[valid_mask]), np.nanmax(turbulence_b_opt[valid_mask]))
axes[1,2].plot([0, mval], [0, mval], 'r--', label='理想拟合')
axes[1,2].set_xlabel('模型a湍流强度'); axes[1,2].set_ylabel('模型b湍流强度'); axes[1,2].legend(); axes[1,2].grid(True)
axes[1,2].set_title('模型对比散点图')
else:
axes[1,2].text(0.5, 0.5, '数据不足', ha='center', va='center', transform=axes[1,2].transAxes)
axes[1,2].set_title('模型对比散点图 (数据不足)')
plt.tight_layout()
outfig = './figure/turbulence_analysis_results.png'
plt.savefig(outfig, dpi=300, bbox_inches='tight')
print(f"图像保存到: {outfig}")
# 保存结果 CSV
results_df = pd.DataFrame({
'高度 (m)': heights,
'模型a湍流强度': turbulence_a,
'模型b原始湍流强度': turbulence_b,
'模型b优化湍流强度': turbulence_b_opt,
'理查逊数': ri,
'湍流动能': tke,
'风切变': wind_shear,
'谱宽': spectrum_widths,
'u (m/s)': u_profiles,
'v (m/s)': v_profiles,
'温度 (K)': temp_k[:len(heights)]
})
csvpath = './result/turbulence_intensity_results.csv'
results_df.to_csv(csvpath, index=False, encoding='utf-8-sig')
print(f"结果 CSV 保存到: {csvpath}")
# 保存评估
with open('./result/model_evaluation.txt', 'w', encoding='utf-8') as f:
f.write(f"原始模型b RMSE: {rmse_orig}\n")
f.write(f"原始模型b R2: {r2_orig}\n")
f.write(f"优化模型b RMSE: {rmse_opt}\n")
f.write(f"优化模型b R2: {r2_opt}\n")
print("处理完成。")
if __name__ == "__main__":
main()
帮我优化一下模型
本文介绍了解决LeetCode上一道题目——除自身以外数组的乘积的方法。该题要求在O(n)的时间复杂度内解决，并且不能使用除法操作。文章提供了两种解法，一种使用额外的空间复杂度O(n)，另一种则达到了常数级别的空间复杂度O(1)。
扫一扫
945
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
