主题013：非平稳风荷载模拟方法

摘要

实际风荷载往往具有显著的非平稳特性，特别是在台风、雷暴等极端天气条件下。本教程系统介绍非平稳随机过程的理论基础、时变功率谱密度、演变谱理论以及非平稳风荷载的模拟方法，包括乘积模型、滤波器模型、小波变换方法和时频分析方法。通过详细的Python数值算例，演示非平稳风荷载的生成、时频特性分析和结构响应计算，为极端风况下的结构风振分析提供完整的技术指南。

关键词：非平稳随机过程；演变谱；时变功率谱；乘积模型；小波变换；时频分析；台风风场

---

1. 引言

1.1 非平稳风荷载的工程背景

在前面的章节中，我们主要讨论了平稳随机风荷载的模拟和分析方法。然而，实际工程中的风荷载往往具有显著的非平稳特性：

（1）台风风场

台风是结构抗风设计中最具挑战性的荷载条件之一。台风风场具有以下非平稳特征：

• 眼墙过境：风速和风向在短时间内剧烈变化
• 强度演变：台风从生成、发展到消亡，风速强度持续变化
• 路径移动：台风中心移动导致风场空间分布时变

（2）雷暴大风

雷暴产生的下击暴流具有极强的非平稳性：

• 突发性强：风速在数秒内急剧增大
• 持续时间短：典型持续时间5-15分钟
• 空间移动：雷暴单体移动导致风场时空变化

（3）阵风边界层

大气边界层中的阵风也具有非平稳特性：

• 湍流强度变化：随大气稳定度变化
• 平均风速变化：日变化和天气系统影响
• 相干结构：大尺度湍流结构的通过

1.2 非平稳性的数学描述

与平稳随机过程不同，非平稳随机过程的统计特性随时间变化：

时变均值：${\mu }_X(t)=E[X(t)]$

时变方差：${\sigma }_X^2(t)=E[(X(t)-{\mu }_X(t){)}^2]$

时变自相关函数：$R_X(t_1,t_2)=E[X(t_1)X(t_2)]$

时变功率谱密度：$S_X(t,\omega )$（演变谱）

1.3 非平稳风荷载模拟的意义

准确模拟非平稳风荷载对于以下工程问题至关重要：

1. 极端风况下的结构安全评估
2. 疲劳寿命的精确预测
3. 风敏感结构的动力响应分析
4. 结构控制系统的优化设计
5. 风能资源的准确评估

---

2. 非平稳随机过程理论基础

2.1 非平稳随机过程的分类

根据非平稳特性的表现形式，非平稳随机过程可分为以下几类：

2.1.1 一阶非平稳（均值非平稳）

过程的均值随时间变化，但协方差结构保持平稳：

$$X(t)=\mu (t)+Y(t)$$

其中$Y(t)$为零均值平稳过程。

2.1.2 二阶非平稳（方差非平稳）

过程的方差或强度随时间变化：

$$X(t)=A(t)\cdot Y(t)$$

其中$A(t)$为确定性包络函数，$Y(t)$为平稳过程。

2.1.3 全非平稳

均值、方差和相关结构都随时间变化：

$$X(t)=\mu (t)+A(t)\cdot Y(g(t))$$

其中$g(t)$为时间变换函数。

2.2 演变谱理论

2.2.1 Priestley演变谱

Priestley（1965）提出了演变谱（Evolutionary Spectrum）的概念，用于描述非平稳随机过程的时频特性。

对于非平稳过程$X(t)$，若可以表示为：

$$X(t)={\int }_{-\infty }^{\infty }A(t,\omega )e^{i\omega t}dZ(\omega )$$

其中$A(t,\omega )$为时频调制函数，$Z(\omega )$为正交增量过程，则演变谱定义为：

$$S(t,\omega )=|A(t,\omega ){|}^{2}S(\omega )$$

其中$S(\omega )$为平稳过程的功率谱密度。

2.2.2 物理意义

演变谱$S(t,\omega )$表示在时刻$t$、频率$\omega$处的能量密度。与平稳过程的功率谱不同，演变谱是时间和频率的二元函数。

2.2.3 时变均方值

由演变谱可计算时变均方值：

$${\sigma }^2(t)={\int }_{-\infty }^{\infty }S(t,\omega )d\omega$$

2.3 时频分析方法

2.3.1 短时傅里叶变换（STFT）

STFT通过加窗傅里叶变换分析信号的局部频谱特性：

$$STFT(t,\omega )={\int }_{-\infty }^{\infty }X(\tau )w(\tau -t)e^{-i\omega \tau }d\tau$$

其中$w(\tau )$为窗函数（如汉宁窗、高斯窗）。

特点：

• 时间分辨率和频率分辨率存在权衡（不确定性原理）
• 适用于缓慢变化的非平稳过程
• 窗函数选择影响分析结果

2.3.2 小波变换

小波变换采用可变带宽的基函数，具有多分辨率特性：

$$WT(t,a)=\frac{1}{\sqrt{a}}{\int }_{-\infty }^{\infty }X(\tau ){\psi }^{\ast }\left(\frac{\tau -t}{a}\right)d\tau$$

其中$\psi$为小波母函数，$a$为尺度参数。

常用小波函数：

• Morlet小波：复指数调制的高斯函数
• Mexican Hat小波：高斯函数的二阶导数
• Daubechies小波：紧支集正交小波

特点：

• 高频处时间分辨率高，低频处频率分辨率高
• 适用于具有瞬态特征的非平稳过程
• 小波基选择影响分析结果

2.3.3 Wigner-Ville分布

Wigner-Ville分布是一种二次型时频分布：

$$W(t,\omega )={\int }_{-\infty }^{\infty }X(t+\tau /2)X^{\ast }(t-\tau /2)e^{-i\omega \tau }d\tau$$

特点：

• 时频分辨率高
• 存在交叉项干扰（多分量信号）
• 可能出现负值

---

3. 非平稳风荷载模拟方法

3.1 乘积模型（Separable Model）

3.1.1 模型形式

乘积模型是最简单的非平稳模型，假设时频调制函数可分离：

$$X(t)=A(t)\cdot Y(t)$$

其中：

• $A(t)$：确定性时变包络函数（调制函数）
• $Y(t)$：零均值平稳随机过程

演变谱为：

$$S(t,\omega )=A^2(t)\cdot S_Y(\omega )$$

3.1.2 包络函数的形式

常用的包络函数形式包括：

（1）指数型

$$A(t)=A_0{\left(\frac{t}{t_{max}}\right)}^{\alpha }\exp \left[-\alpha \left(\frac{t}{t_{max}}-1\right)\right]$$

适用于描述台风眼墙过境时的风速变化。

（2）高斯型

$$A(t)=A_0\exp \left[-\frac{(t-t_c{)}^2}{2{\sigma }_t^2}\right]$$

适用于描述阵风或下击暴流。

（3）分段函数

$$A(t)=\{\begin{array}{ll}{A}_0{\left(\frac{t}{t_1}\right)}^{{\alpha }_1}& 0\le t<t_1\\ A_0& t_1\le t<t_2\\ A_0\exp \left[-\beta (t-t_2)\right]& t\ge t_2\end{array}$$

适用于描述完整的台风过程（增强-稳定-衰减）。

3.1.3 模拟步骤

步骤1：根据工程经验或实测数据确定包络函数$A(t)$

步骤2：生成平稳随机过程$Y(t)$（如采用谐波叠加法或AR模型）

步骤3：计算非平稳过程$X(t)=A(t)\cdot Y(t)$

3.2 滤波器模型（Filtered Model）

3.2.1 模型原理

滤波器模型通过时变线性滤波器对白噪声进行滤波，生成非平稳过程：

$$X(t)={\int }_{-\infty }^{t}h(t,\tau )W(\tau )d\tau$$

其中：

• $h(t,\tau )$：时变脉冲响应函数
• $W(t)$：白噪声过程

3.2.2 时变AR模型

时变自回归（TV-AR）模型：

$$X(t)=-\sum _{k=1}^p{a}_k(t)X(t-k\Delta t)+\sigma (t)W(t)$$

其中$a_k(t)$和$\sigma (t)$为时变参数。

参数识别：

• 滑动窗口法：在每个时间窗口内假设局部平稳
• 自适应算法：如最小均方（LMS）算法
• 基函数展开：$a_k(t)={\sum }_j{c}_{kj}{\varphi }_j(t)$

3.3 谱表示法的扩展

3.3.1 非均匀调制谱表示

将谐波叠加法扩展到非平稳情况：

$$X(t)=\sum _{k=1}^{N}A(t,{\omega }_k)\sqrt{2S({\omega }_k)\Delta \omega }\cos ({\omega }_{k}t+{\varphi }_k)$$

其中$A(t,{\omega }_k)$为时频调制函数。

3.3.2 多点非平稳风场

对于空间多点非平稳风场：

$$U_i(t)={\bar{U}}_i(t)+u_i(t)$$

其中脉动分量：

$$u_i(t)=\sum _{k=1}^N{A}_i(t,{\omega }_k)\sqrt{2S_{ii}({\omega }_k)\Delta \omega }\cos ({\omega }_{k}t+{\theta }_{ik})$$

相位角${\theta }_{ik}$需满足空间相干条件。

3.4 基于小波变换的模拟方法

3.4.1 小波重构法

步骤1：对目标演变谱进行小波分解

步骤2：生成各尺度上的小波系数（随机过程）

步骤3：通过小波逆变换重构非平稳过程

3.4.2 小波包方法

小波包提供更精细的时频分解，适用于复杂非平稳过程：

$$X(t)=\sum _{j,k}{d}_{j,k}{\psi }_{j,k}(t)$$

其中$d_{j,k}$为小波包系数，可根据目标时频特性生成。

---

4. Python数值算例

以下Python代码演示了非平稳风荷载模拟方法的实现，包括乘积模型、时频分析、小波变换等。

"""
非平稳风荷载模拟方法
主题013：非平稳风荷载模拟方法
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.interpolate import interp1d
import pywt
import warnings
warnings.filterwarnings('ignore')

# 设置matplotlib后端为Agg，不弹出窗口
plt.switch_backend('Agg')
plt.rcParams['font.size'] = 10
plt.rcParams['figure.dpi'] = 100

# ==================== 1. 包络函数定义 ====================

def envelope_exponential(t, t_max, alpha, A0=1.0):
    """
    指数型包络函数
    
    参数:
        t: 时间数组
        t_max: 最大强度时刻
        alpha: 形状参数
        A0: 最大强度
    
    返回:
        A: 包络函数值
    """
    A = A0 * (t / t_max)**alpha * np.exp(-alpha * (t / t_max - 1))
    A[t < 0] = 0
    return A

def envelope_gaussian(t, tc, sigma_t, A0=1.0):
    """
    高斯型包络函数
    
    参数:
        t: 时间数组
        tc: 中心时刻
        sigma_t: 时间宽度
        A0: 最大强度
    
    返回:
        A: 包络函数值
    """
    return A0 * np.exp(-(t - tc)**2 / (2 * sigma_t**2))

def envelope_piecewise(t, t1, t2, alpha1, beta, A0=1.0):
    """
    分段型包络函数（台风模型）
    
    参数:
        t: 时间数组
        t1: 增强阶段结束时刻
        t2: 稳定阶段结束时刻
        alpha1: 增强阶段指数
        beta: 衰减系数
        A0: 最大强度
    
    返回:
        A: 包络函数值
    """
    A = np.zeros_like(t)
    # 增强阶段
    mask1 = (t >= 0) & (t < t1)
    A[mask1] = A0 * (t[mask1] / t1)**alpha1
    # 稳定阶段
    mask2 = (t >= t1) & (t < t2)
    A[mask2] = A0
    # 衰减阶段
    mask3 = t >= t2
    A[mask3] = A0 * np.exp(-beta * (t[mask3] - t2))
    return A

def envelope_typhoon_eyewall(t, t_center, width, A0=1.0):
    """
    台风眼墙过境包络函数
    
    参数:
        t: 时间数组
        t_center: 眼墙中心时刻
        width: 眼墙宽度
        A0: 最大强度
    
    返回:
        A: 包络函数值
    """
    # 使用双曲正割函数模拟眼墙结构
    A = A0 / np.cosh((t - t_center) / width)**2
    return A

# ==================== 2. 平稳风场生成 ====================

def generate_stationary_wind(T, dt, f_cutoff, spectrum_func, seed=None):
    """
    生成平稳脉动风（谐波叠加法）
    
    参数:
        T: 总时长
        dt: 时间步长
        f_cutoff: 截止频率
        spectrum_func: 谱函数
        seed: 随机种子
    
    返回:
        t: 时间数组
        u: 脉动风速
    """
    if seed is not None:
        np.random.seed(seed)
    
    t = np.arange(0, T, dt)
    N = len(t)
    
    # 频率离散
    df = 1 / T
    f = np.arange(df, f_cutoff, df)
    
    # 生成随机相位
    phi = np.random.uniform(0, 2*np.pi, len(f))
    
    # 谐波叠加
    u = np.zeros(N)
    for i, f_i in enumerate(f):
        S_i = spectrum_func(f_i)
        amplitude = np.sqrt(2 * S_i * df)
        u += amplitude * np.cos(2 * np.pi * f_i * t + phi[i])
    
    return t, u

def kaimal_spectrum(f, z=50, U_mean=30, u_star=2.0):
    """Kaimal顺风向脉动风速谱"""
    x = f * z / U_mean
    S_u = u_star**2 * 200 * x / (f * (1 + 50 * x)**(5/3))
    return S_u

# ==================== 3. 非平稳风场生成 ====================

def generate_nonstationary_wind_product(T, dt, envelope_func, spectrum_func, seed=None):
    """
    乘积模型生成非平稳风场
    
    参数:
        T: 总时长
        dt: 时间步长
        envelope_func: 包络函数
        spectrum_func: 谱函数
        seed: 随机种子
    
    返回:
        t: 时间数组
        U: 非平稳风速
        A: 包络函数
        u_stationary: 平稳分量
    """
    t = np.arange(0, T, dt)
    
    # 计算包络
    A = envelope_func(t)
    
    # 生成平稳过程
    _, u_stationary = generate_stationary_wind(T, dt, 5.0, spectrum_func, seed)
    
    # 调制
    U = A * u_stationary
    
    return t, U, A, u_stationary

def generate_typhoon_wind(T, dt, t_center, max_wind, radius_max_wind, translation_speed, 
                         spectrum_func, seed=None):
    """
    生成台风风场（简化模型）
    
    参数:
        T: 总时长
        dt: 时间步长
        t_center: 台风中心过境时刻
        max_wind: 最大风速 (m/s)
        radius_max_wind: 最大风速半径 (km)
        translation_speed: 移动速度 (m/s)
        spectrum_func: 脉动风谱函数
        seed: 随机种子
    
    返回:
        t: 时间数组
        U_mean: 时变平均风速
        u_fluctuating: 脉动风速
        U_total: 总风速
    """
    t = np.arange(0, T, dt)
    
    # 计算结构到台风中心的距离随时间变化
    # 假设结构位于台风路径右侧，台风从左向右移动
    distance = np.abs((t - t_center) * translation_speed)  # m
    distance_km = distance / 1000  # km
    
    # Holland风场模型计算平均风速
    B = 1.5  # Holland参数
    rho = 1.225  # 空气密度
    dp = 50e3  # 中心气压差 (Pa)
    
    # 简化Holland模型
    r_max = radius_max_wind  # km
    U_mean = np.zeros_like(t)
    
    for i, r in enumerate(distance_km):
        if r < 1:  # 台风眼内
            U_mean[i] = max_wind * (r / r_max)
        else:
            # Holland风场
            U_mean[i] = max_wind * (r_max / r)**B * np.exp(1 - (r_max / r)**B)
    
    # 生成脉动风（使用乘积模型）
    # 湍流强度随平均风速变化
    I_u = 0.1 + 0.05 * (U_mean / 30)  # 简化模型
    
    _, u_stationary = generate_stationary_wind(T, dt, 5.0, spectrum_func, seed)
    
    # 调制脉动风
    u_fluctuating = I_u * U_mean * u_stationary / np.std(u_stationary)
    
    U_total = U_mean + u_fluctuating
    
    return t, U_mean, u_fluctuating, U_total

# ==================== 4. 时频分析 ====================

def stft_analysis(signal, fs, window='hann', nperseg=256, noverlap=None):
    """
    短时傅里叶变换分析
    
    参数:
        signal: 信号
        fs: 采样频率
        window: 窗函数
        nperseg: 每段长度
        noverlap: 重叠长度
    
    返回:
        f: 频率数组
        t: 时间数组
        Zxx: STFT结果
    """
    if noverlap is None:
        noverlap = nperseg // 2
    
    f, t, Zxx = signal.stft(signal, fs, window=window, nperseg=nperseg, noverlap=noverlap)
    return f, t, Zxx

def wavelet_analysis(signal, scales, wavelet='cmor1.5-1.0', fs=1.0):
    """
    连续小波变换分析
    
    参数:
        signal: 信号
        scales: 尺度数组
        wavelet: 小波函数
        fs: 采样频率
    
    返回:
        coef: 小波系数
        freqs: 频率数组
    """
    coef, freqs = pywt.cwt(signal, scales, wavelet, sampling_period=1/fs)
    return coef, freqs

def compute_instantaneous_frequency(signal, fs):
    """
    计算瞬时频率（通过Hilbert变换）
    
    参数:
        signal: 信号
        fs: 采样频率
    
    返回:
        t: 时间数组
        inst_freq: 瞬时频率
        amplitude_envelope: 包络
    """
    # Hilbert变换
    analytic_signal = signal.hilbert(signal)
    
    # 瞬时幅值（包络）
    amplitude_envelope = np.abs(analytic_signal)
    
    # 瞬时相位
    instantaneous_phase = np.unwrap(np.angle(analytic_signal))
    
    # 瞬时频率
    inst_freq = np.diff(instantaneous_phase) / (2.0 * np.pi) * fs
    
    t = np.arange(len(signal)) / fs
    
    return t[:-1], inst_freq, amplitude_envelope

# ==================== 5. 算例1：不同包络函数对比 ====================

print("="*60)
print("算例1：不同包络函数对比")
print("="*60)

T = 600  # 总时长 (s)
dt = 0.1  # 时间步长
t = np.arange(0, T, dt)

# 定义不同的包络函数
A_exp = envelope_exponential(t, 300, 2.0, 1.0)
A_gauss = envelope_gaussian(t, 300, 100, 1.0)
A_piece = envelope_piecewise(t, 200, 400, 1.5, 0.01, 1.0)
A_typhoon = envelope_typhoon_eyewall(t, 300, 30, 1.0)

print("包络函数类型:")
print("  1. 指数型包络")
print("  2. 高斯型包络")
print("  3. 分段型包络（台风模型）")
print("  4. 台风眼墙包络")

# 绘制包络函数
fig, axes = plt.subplots(2, 2, figsize=(12, 8))

axes[0, 0].plot(t, A_exp, 'b-', linewidth=1.5)
axes[0, 0].set_xlabel('时间 (s)')
axes[0, 0].set_ylabel('包络 A(t)')
axes[0, 0].set_title('指数型包络')
axes[0, 0].grid(True, alpha=0.3)

axes[0, 1].plot(t, A_gauss, 'r-', linewidth=1.5)
axes[0, 1].set_xlabel('时间 (s)')
axes[0, 1].set_ylabel('包络 A(t)')
axes[0, 1].set_title('高斯型包络')
axes[0, 1].grid(True, alpha=0.3)

axes[1, 0].plot(t, A_piece, 'g-', linewidth=1.5)
axes[1, 0].set_xlabel('时间 (s)')
axes[1, 0].set_ylabel('包络 A(t)')
axes[1, 0].set_title('分段型包络（台风模型）')
axes[1, 0].grid(True, alpha=0.3)

axes[1, 1].plot(t, A_typhoon, 'm-', linewidth=1.5)
axes[1, 1].set_xlabel('时间 (s)')
axes[1, 1].set_ylabel('包络 A(t)')
axes[1, 1].set_title('台风眼墙包络')
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('envelope_functions.png', dpi=150, bbox_inches='tight')
print("\n图形已保存: envelope_functions.png")
plt.close()

# ==================== 6. 算例2：非平稳风场生成（乘积模型） ====================

print("\n" + "="*60)
print("算例2：非平稳风场生成（乘积模型）")
print("="*60)

# 使用台风眼墙包络生成非平稳风场
envelope_func = lambda t: envelope_typhoon_eyewall(t, 300, 30, 5.0)

t, U_nonstat, A, u_stat = generate_nonstationary_wind_product(
    T, dt, envelope_func, lambda f: kaimal_spectrum(f, z=50, U_mean=30, u_star=2.0), seed=42)

print(f"模拟参数:")
print(f"  总时长: {T} s")
print(f"  时间步长: {dt} s")
print(f"  数据点数: {len(t)}")
print(f"\n统计特性:")
print(f"  平稳分量标准差: {np.std(u_stat):.4f} m/s")
print(f"  非平稳过程标准差: {np.std(U_nonstat):.4f} m/s")
print(f"  最大包络值: {np.max(A):.4f}")
print(f"  最大风速: {np.max(U_nonstat):.4f} m/s")

# 计算时变标准差（滑动窗口）
window_size = int(50 / dt)  # 50秒窗口
sigma_t = np.array([np.std(U_nonstat[max(0, i-window_size//2):min(len(U_nonstat), i+window_size//2)]) 
                    for i in range(len(U_nonstat))])

# 绘制结果
fig, axes = plt.subplots(3, 1, figsize=(12, 10))

# 包络函数
axes[0].plot(t, A, 'b-', linewidth=1.5, label='包络函数')
axes[0].set_ylabel('A(t)')
axes[0].set_title('包络函数')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# 平稳分量
axes[1].plot(t, u_stat, 'g-', linewidth=0.5, alpha=0.7, label='平稳分量')
axes[1].set_ylabel('u(t) (m/s)')
axes[1].set_title('平稳脉动风分量')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

# 非平稳过程
axes[2].plot(t, U_nonstat, 'r-', linewidth=0.5, alpha=0.7, label='非平稳风速')
axes[2].plot(t, A * np.std(u_stat), 'k--', linewidth=2, label='理论包络×σ')
axes[2].plot(t, sigma_t, 'm-', linewidth=1.5, label='滑动窗口σ')
axes[2].set_xlabel('时间 (s)')
axes[2].set_ylabel('U(t) (m/s)')
axes[2].set_title('非平稳脉动风')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('nonstationary_wind_generation.png', dpi=150, bbox_inches='tight')
print("\n图形已保存: nonstationary_wind_generation.png")
plt.close()

# ==================== 7. 算例3：时频分析 ====================

print("\n" + "="*60)
print("算例3：非平稳风场时频分析")
print("="*60)

fs = 1 / dt  # 采样频率

# STFT分析
print("进行STFT分析...")
f_stft, t_stft, Zxx = stft_analysis(U_nonstat, fs, nperseg=256)

# 小波分析
print("进行小波分析...")
scales = np.arange(1, 128)
coef, freqs_wavelet = wavelet_analysis(U_nonstat, scales, 'cmor1.5-1.0', fs)

# 瞬时频率分析
print("计算瞬时频率...")
t_inst, inst_freq, amp_envelope = compute_instantaneous_frequency(U_nonstat, fs)

print("\n时频分析完成")
print(f"  STFT频率范围: {f_stft[0]:.3f} - {f_stft[-1]:.3f} Hz")
print(f"  小波频率范围: {freqs_wavelet[-1]:.3f} - {freqs_wavelet[0]:.3f} Hz")

# 绘制时频分析结果
fig = plt.figure(figsize=(14, 10))

# 原始信号
ax1 = plt.subplot(4, 1, 1)
ax1.plot(t, U_nonstat, 'b-', linewidth=0.5)
ax1.set_ylabel('U(t) (m/s)')
ax1.set_title('非平稳风速时程')
ax1.grid(True, alpha=0.3)

# STFT时频图
ax2 = plt.subplot(4, 1, 2)
im2 = ax2.pcolormesh(t_stft, f_stft, np.abs(Zxx), shading='gouraud', cmap='jet')
ax2.set_ylabel('频率 (Hz)')
ax2.set_title('短时傅里叶变换 (STFT)')
ax2.set_ylim([0, 2])
plt.colorbar(im2, ax=ax2, label='幅值')

# 小波时频图
ax3 = plt.subplot(4, 1, 3)
im3 = ax3.pcolormesh(t, freqs_wavelet, np.abs(coef), shading='gouraud', cmap='jet')
ax3.set_ylabel('频率 (Hz)')
ax3.set_title('连续小波变换 (CWT)')
ax3.set_ylim([0, 2])
plt.colorbar(im3, ax=ax3, label='幅值')

# 瞬时频率
ax4 = plt.subplot(4, 1, 4)
ax4.plot(t_inst, inst_freq, 'g-', linewidth=0.5)
ax4.set_xlabel('时间 (s)')
ax4.set_ylabel('瞬时频率 (Hz)')
ax4.set_title('瞬时频率（Hilbert变换）')
ax4.set_ylim([0, 2])
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('time_frequency_analysis.png', dpi=150, bbox_inches='tight')
print("\n图形已保存: time_frequency_analysis.png")
plt.close()

# ==================== 8. 算例4：台风风场模拟 ====================

print("\n" + "="*60)
print("算例4：台风风场模拟")
print("="*60)

T_typhoon = 3600  # 1小时
dt_typhoon = 1.0  # 1秒时间步长

# 台风参数
t_center = 1800  # 台风中心过境时刻 (s)
max_wind = 50  # 最大风速 (m/s)
radius_max_wind = 30  # 最大风速半径 (km)
translation_speed = 15  # 台风移动速度 (m/s)

t_typhoon, U_mean, u_fluct, U_total = generate_typhoon_wind(
    T_typhoon, dt_typhoon, t_center, max_wind, radius_max_wind, 
    translation_speed, lambda f: kaimal_spectrum(f, z=50, U_mean=30, u_star=2.0), seed=123)

print(f"台风参数:")
print(f"  中心过境时刻: {t_center} s")
print(f"  最大风速: {max_wind} m/s")
print(f"  最大风速半径: {radius_max_wind} km")
print(f"  移动速度: {translation_speed} m/s")
print(f"\n风速统计:")
print(f"  平均风速范围: {np.min(U_mean):.2f} - {np.max(U_mean):.2f} m/s")
print(f"  总风速范围: {np.min(U_total):.2f} - {np.max(U_total):.2f} m/s")
print(f"  平均风速标准差: {np.std(U_mean):.2f} m/s")
print(f"  总风速标准差: {np.std(U_total):.2f} m/s")

# 绘制台风风场
fig, axes = plt.subplots(3, 1, figsize=(12, 10))

# 平均风速
axes[0].plot(t_typhoon, U_mean, 'b-', linewidth=1.5)
axes[0].axvline(x=t_center, color='r', linestyle='--', alpha=0.5, label='中心过境')
axes[0].set_ylabel('平均风速 (m/s)')
axes[0].set_title('台风平均风速（Holland模型）')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# 脉动风速
axes[1].plot(t_typhoon, u_fluct, 'g-', linewidth=0.3, alpha=0.7)
axes[1].axvline(x=t_center, color='r', linestyle='--', alpha=0.5)
axes[1].set_ylabel('脉动风速 (m/s)')
axes[1].set_title('台风脉动风速分量')
axes[1].grid(True, alpha=0.3)

# 总风速
axes[2].plot(t_typhoon, U_total, 'r-', linewidth=0.3, alpha=0.7)
axes[2].plot(t_typhoon, U_mean, 'b-', linewidth=1.5, label='平均风速')
axes[2].axvline(x=t_center, color='r', linestyle='--', alpha=0.5, label='中心过境')
axes[2].set_xlabel('时间 (s)')
axes[2].set_ylabel('总风速 (m/s)')
axes[2].set_title('台风总风速')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('typhoon_wind_simulation.png', dpi=150, bbox_inches='tight')
print("\n图形已保存: typhoon_wind_simulation.png")
plt.close()

# ==================== 9. 算例5：非平稳结构响应分析 ====================

print("\n" + "="*60)
print("算例5：非平稳风荷载下的结构响应")
print("="*60)

# 单自由度结构参数
m = 1000.0  # 质量 (kg)
f_n = 1.0  # 固有频率 (Hz)
omega_n = 2 * np.pi * f_n
k = m * omega_n**2
c = 2 * m * omega_n * 0.05  # 5%阻尼比

# 使用台风风场计算结构响应（简化计算）
# 假设结构位于50m高度，迎风面积10m²
rho = 1.225
A = 10.0
C_d = 1.2

# 计算风荷载
F_wind = 0.5 * rho * U_mean**2 * A * C_d + rho * U_mean * A * C_d * u_fluct

# 数值积分计算响应（Newmark-beta法）
def newmark_beta(M, C, K, F, dt, beta=0.25, gamma=0.5):
    """Newmark-beta法计算动力响应"""
    N = len(F)
    x = np.zeros(N)
    v = np.zeros(N)
    a = np.zeros(N)
    
    # 初始条件
    a[0] = (F[0] - C * v[0] - K * x[0]) / M
    
    for i in range(N - 1):
        # 预测
        x_pred = x[i] + dt * v[i] + (0.5 - beta) * dt**2 * a[i]
        v_pred = v[i] + (1 - gamma) * dt * a[i]
        
        # 求解加速度
        a[i+1] = (F[i+1] - C * v_pred - K * x_pred) / (M + gamma * dt * C + beta * dt**2 * K)
        
        # 校正
        v[i+1] = v_pred + gamma * dt * a[i+1]
        x[i+1] = x_pred + beta * dt**2 * a[i+1]
    
    return x, v, a

# 计算响应
x_disp, x_vel, x_acc = newmark_beta(m, c, k, F_wind, dt_typhoon)

print(f"结构参数:")
print(f"  质量: {m} kg")
print(f"  刚度: {k:.2f} N/m")
print(f"  固有频率: {f_n} Hz")
print(f"  阻尼比: 5%")
print(f"\n响应统计:")
print(f"  最大位移: {np.max(np.abs(x_disp)):.4f} m")
print(f"  最大速度: {np.max(np.abs(x_vel)):.4f} m/s")
print(f"  最大加速度: {np.max(np.abs(x_acc)):.4f} m/s²")
print(f"  位移标准差: {np.std(x_disp):.4f} m")

# 绘制响应
fig, axes = plt.subplots(3, 1, figsize=(12, 10))

# 风荷载
axes[0].plot(t_typhoon, F_wind / 1000, 'b-', linewidth=0.3, alpha=0.7)
axes[0].axvline(x=t_center, color='r', linestyle='--', alpha=0.5)
axes[0].set_ylabel('风荷载 (kN)')
axes[0].set_title('台风风荷载')
axes[0].grid(True, alpha=0.3)

# 位移响应
axes[1].plot(t_typhoon, x_disp * 1000, 'g-', linewidth=0.3, alpha=0.7)
axes[1].axvline(x=t_center, color='r', linestyle='--', alpha=0.5)
axes[1].set_ylabel('位移 (mm)')
axes[1].set_title('结构位移响应')
axes[1].grid(True, alpha=0.3)

# 加速度响应
axes[2].plot(t_typhoon, x_acc, 'r-', linewidth=0.3, alpha=0.7)
axes[2].axvline(x=t_center, color='r', linestyle='--', alpha=0.5)
axes[2].set_xlabel('时间 (s)')
axes[2].set_ylabel('加速度 (m/s²)')
axes[2].set_title('结构加速度响应')
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('structural_response_nonstationary.png', dpi=150, bbox_inches='tight')
print("\n图形已保存: structural_response_nonstationary.png")
plt.close()

print("\n" + "="*60)
print("非平稳风荷载模拟完成！")
print("="*60)

---