平稳高斯随机过程：Rice 水平穿越 + 峰值期望全套推导

平稳高斯随机过程：Rice 水平穿越 + 峰值期望全套推导

适用：随机振动、汽车耐久载荷谱、疲劳极值预估；所有变量沿用前文定义
\(\(X(t)\)\)：零均值平稳高斯过程；\(\dot \(X(t)\)=\frac{dX}{dt}\)；\(\ddot \(X(t)\)=\frac{d^2X}{dt^2}\)
\(\sigma_X^2=E[X^2],\ \sigma_{\dot X}^2=E[\dot X^2],\ \sigma_{\ddot X}^2=E[\ddot X^2]\)

一、回顾：向上水平穿越率（Upcrossing rate）

单位时间内，信号从下往上跨过阈值\(\alpha\) 的平均次数：

\[N^+(\alpha)=\int_{0}^{\infty}\dot{x}\,p(\alpha,\dot{x})\,d\dot{x} \]

零均值平稳高斯下，\(X\)与\(\dot X\)独立，联合密度：

\[p(x,\dot x)=\frac{1}{2\pi\sigma_X\sigma_{\dot X}}\exp\left[-\frac{x^2}{2\sigma_X^2}-\frac{\dot x^2}{2\sigma_{\dot X}^2}\right] \]

代入\(x=\alpha\)积分，得到Rice upcrossing公式：

\[\boldsymbol{N^+(\alpha)=\frac{\sigma_{\dot X}}{2\pi\sigma_X}\exp\left(-\frac{\alpha^2}{2\sigma_X^2}\right)} \]

• \(\boldsymbol{\alpha=0}\)：零点向上穿越率

\[N_0=N^+(0)=\frac{\sigma_{\dot X}}{2\pi\sigma_X} \]

\(N_0\)：平均零穿越频率（单位 Hz），工程上常用来表征随机振动的平均振荡频率。

---

二、峰值期望率（Peak rate，单位时间平均峰值个数）

峰值定义：在时刻\(t\)满足 \(\dot \(X(t)\)=0,\ \ddot \(X(t)\)<0\)（导数为0，二阶导负，局部极大值）。

峰值期望率公式定义

\[N_p=\int_{-\infty}^{0}|\ddot{x}|\,p(0,\ddot{x})\,d\ddot{x} \]

\(p(\dot x,\ddot x)\)：\(\dot \(X(t)\)\)和\(\ddot \(X(t)\)\)的联合概率密度；积分只取\(\ddot x<0\)（极大峰，排除谷值）。

对于零均值平稳高斯过程：
\(\dot \(X(t)\)\)同样是零均值高斯；\(\dot X\)与\(\ddot X\)相互独立（和前面\(X,\dot X\)独立道理完全一样）。

\[E[\dot \(X(t)\)\ddot \(X(t)\)]=\frac12\frac{d}{dt}E\left[\dot X^2(t)\right]=0 \]

平稳，方差\(E[\dot X^2]\)不随时间变化，时间导数为0 → \(\dot X,\ddot X\)不相关，高斯变量不相关即独立。

联合密度：

\[p(\dot x,\ddot x)=\frac{1}{2\pi\sigma_{\dot X}\sigma_{\ddot X}}\exp\left[-\frac{\dot x^2}{2\sigma_{\dot X}^2}-\frac{\ddot x^2}{2\sigma_{\ddot X}^2}\right] \]

峰值条件\(\dot x=0\)，代入：

\[p(0,\ddot{x})=\frac{1}{2\pi\sigma_{\dot X}\sigma_{\ddot X}}\exp\left[-\frac{0}{2\sigma_{\dot X}^2}-\frac{\ddot x^2}{2\sigma_{\ddot X}^2}\right] \]

代入峰值率积分：

\[N_p=\int_{-\infty}^{0} |\ddot x| \cdot \frac{1}{2\pi\sigma_{\dot X}\sigma_{\ddot X}} \exp\left(-\frac{\ddot x^2}{2\sigma_{\ddot X}^2}\right) d\ddot x \]

令 \(\ddot x=-u,\ u>0,\ d\ddot{x}=-du\)：

\[N_p=\int_{0}^{\infty} u \cdot \frac{1}{2\pi\sigma_{\dot X}\sigma_{\ddot X}} \exp\left(-\frac{u^2}{2\sigma_{\ddot X}^2}\right) du \]

积分 \(\displaystyle\int_{0}^{\infty}u\exp\left(-\frac{u^2}{2\sigma_{\ddot X}^2}\right)du=\sigma_{\ddot X}^2\)

\[N_p=\frac{\sigma_{\ddot X}^2}{2\pi\sigma_{\dot X}\sigma_{\ddot X}}=\boldsymbol{\frac{\sigma_{\ddot X}}{2\pi\sigma_{\dot X}}} \]

\(N_p\)：单位时间内局部极大峰值的平均数量。

三、不规则因子（Irregularity factor，非常关键！）

\[\boldsymbol{\gamma=\frac{N_0}{N_p}=\frac{\sigma_{\dot X}^2}{\sigma_X\sigma_{\ddot X}}} \]

• \(\gamma\to1\)：窄带随机过程；峰值几乎全部在零点穿越附近，峰和零点一一对应；峰值服从瑞利分布；
• \(\gamma\ll1\)：宽带随机过程；一个周期内出现很多小峰，峰数量 > 零点穿越数量；疲劳计算必须考虑宽带修正（Dirlik模型就是基于这个）。

窄带：\(\gamma\approx1\)，\(N_p\approx N_0\)；宽带：\(\gamma<1\)，\(N_p>N_0\)

---

四、用功率谱密度PSD表达所有量（工程最常用形式）

单边功率谱 \(S_X(f),\ f\ge0\)，单位 \(\text{X}^2/\text{Hz}\)

\[\begin{aligned} \sigma_X^2 &= \int_{0}^{\infty} S_X(f) df \\ \sigma_{\dot X}^2 &= \int_{0}^{\infty}(2\pi f)^2 S_X(f) df \\ \sigma_{\ddot X}^2 &= \int_{0}^{\infty}(2\pi f)^4 S_X(f) df \end{aligned} \]

记谱矩（Spectral moments）：

\[m_k=\int_{0}^{\infty} (2\pi f)^k S_X(f) df \]

则：

\[\sigma_X^2=m_0,\quad \sigma_{\dot X}^2=m_2,\quad \sigma_{\ddot X}^2=m_4 \]

把上面全部公式换成谱矩形式：

\[\begin{aligned} N_0 &= \frac{1}{2\pi}\sqrt{\frac{m_2}{m_0}} \\ N_p &= \frac{1}{2\pi}\sqrt{\frac{m_4}{m_2}} \\ \gamma &= \frac{N_0}{N_p}=\sqrt{\frac{m_2^2}{m_0 m_4}} \end{aligned} \]

Rice向上穿越率写成谱矩：

\[N^+(\alpha)=N_0 \exp\left(-\frac{\alpha^2}{2m_0}\right) \]

谱矩 \(m_0,m_2,m_4\) 是Dirlik、Wirsching-Light等随机疲劳应力谱模型的输入参数，耐久载荷谱里必用。

---

五、峰值幅值的概率分布（Rice峰值分布）

单位时间内，幅值大于\(\alpha\)的峰值期望数目：

\[N_p(\alpha)=N_p\cdot \frac{1}{\sqrt{2\pi}\gamma}\int_{\alpha/\sqrt{m_0}}^{\infty} \exp\left(-\frac{x^2}{2\gamma^2}\right)\int_{0}^{x/\gamma}\exp\left(-\frac{y^2}{2}\right)dy \ dx \]

窄带极限 \(\gamma=1\)：Rice分布退化为瑞利分布：

\[p_A(a)=\frac{a}{m_0}\exp\left(-\frac{a^2}{2m_0}\right),\quad a\ge0 \]

窄带下，超过\(\alpha\)峰值个数：

\[N_p(\alpha)=N_p\exp\left(-\frac{\alpha^2}{2m_0}\right)=N_0\exp\left(-\frac{\alpha^2}{2m_0}\right)=N^+(\alpha) \]

物理含义：窄带，每一次向上穿越\(\alpha\)，对应一个高度大于\(\alpha\)的峰值。

宽带情况：\(\gamma<1\)，不能用瑞利；要用Dirlik近似，Dirlik把峰值分布拆成瑞利+两个高斯分量，是汽车振动疲劳载荷谱工业标准。

---

六、Python代码实现：PSD计算谱矩、\(N_0,N_p,\gamma\)，Rice穿越率

import numpy as np

def spectral_moments(f, Sx, k):
    """
    计算单边PSD的谱矩 m_k = ∫ (2π f)^k * Sx(f) df
    f:频率向量 Hz
    Sx:单边功率谱密度向量
    k:谱矩阶数
    """
    df = np.diff(f, prepend=f[0])
    integrand = ((2*np.pi*f)**k) * Sx
    mk = np.sum(integrand * df)
    return mk

def rice_upcrossing(alpha, m0, N0):
    """Rice向上穿越率 N+(\(\alpha\))"""
    return N0 * np.exp(-alpha**2/(2*m0))

# ========== Demo：窄带PSD示例 ==========
f = np.linspace(0, 200, 2000)
# 窄带：中心频率50Hz
Sx = np.zeros_like(f)
mask = np.logical_and(f>48, f<52)
Sx[mask] = 1.0

m0 = spectral_moments(f, Sx, 0)
m2 = spectral_moments(f, Sx, 2)
m4 = spectral_moments(f, Sx, 4)

N0 = 1/(2*np.pi)*np.sqrt(m2/m0)
Np = 1/(2*np.pi)*np.sqrt(m4/m2)
gamma = N0/Np

print(f"m0={m0:.3f}, m2={m2:.3f}, m4={m4:.3f}")
print(f"N0零点穿越率={N0:.2f} Hz")
print(f"Np峰值率={Np:.2f} Hz")
print(f"不规则因子\(\gamma\)={gamma:.3f}")

alpha = np.sqrt(m0)
Na = rice_upcrossing(alpha, m0, N0)
print(f"阈值 \(\alpha\)={alpha:.2f}，向上穿越率 N+(\(\alpha\))={Na:.3f} /s")

---

七、扩展：Dirlik峰值分布（宽带疲劳，载荷谱雨流）

Dirlik PDF：

\[p(a)=\frac{1}{\sqrt{m_0}}\left[ \frac{D_1}{Q}\exp\left(-\frac{a}{Q\sqrt{m_0}}\right)+ \frac{D_2 a}{R^2\sqrt{m_0}}\exp\left(-\frac{a^2}{2R^2 m_0}\right)+ D_3\frac{a}{\sqrt{m_0}}\exp\left(-\frac{a^2}{2m_0}\right) \right] \]

系数由\(\gamma\)和\(\lambda_2=\sqrt{(1-\gamma^2)}\)算出，输出雨流幅值分布，直接用于Miner线性累积损伤计算。

如果你需要，我可以：

1. 单独推导Dirlik模型完整公式，含系数推导；
2. 扩展代码：Dirlik PDF + Miner损伤计算；
3. 或者做NCode/LMS TestLab对应的谱矩输入对照。你想继续哪一部分？