穿越分析的详解11
穿越分析的基本结论
正穿越 \(y=a\)的平均频率一般公式为：
\[\nu_a^+ = \int_0^\infty p(a,\dot y)\dot y\,d\dot y
\]
对高斯过程，代入联合密度后得：
\[\nu_a^+ = \frac{1}{2\pi}\frac{\sigma_{\dot y}}{\sigma_y} e^{-a^2/2\sigma_y^2}
\]
特别地，正穿越零的频率（即过程平均频率）为：
\[\nu_0^+ = \frac{1}{2\pi}\frac{\sigma_{\dot y}}{\sigma_y}
\]
5.1 穿越分析
5.1.1 正穿越概率→莱斯公式
无论是时域疲劳还是频域疲劳，都可以通过求解幅值分布来计算，不同的是计算方式。
而幅值分布可以看成如下图所示，求取峰值的分布。为了进行一般性的推导，我们先求解y>a的分布。
窄带过程的典型样本
上图中，y(t)有三次穿越y=a。
在dt时间内，y穿越y=a的条件：
\[\begin{aligned}
\text{i.}\quad & \dot{y} > 0 \\
\text{ii.}\quad & y < a \\
\text{iii.}\quad & y > a - \dot{y}dt
\end{aligned}\]
第一个条件保证是正向穿越；第二体条件保证y在a之下；
第三个是穿越速度的要求，y在t+dt时刻能够穿越a，需要保证y在t时刻的斜率大于 \(（a-y）/dt\)，
即
\[\dot{y}>(a-y)/dt
\]
因此
\[y<a-\dot{y}dt
\]
合并条件ii和iii，得到
\[\begin{aligned}
\text{i.}\quad & \dot{y} > 0 \\
\text{ii.}\quad & a - \dot{y}dt<y < a \\\
\end{aligned}\]
基于上述条件，发生正穿越a的概率为y和\(\dot{y}\)的联合分布中如图楔形部分所示的体积。
积分区域：
$$
\begin{aligned}
\left(\begin{aligned}\text{在 } dt \text{ 时间内}\\\text{发生正穿越}\text{a}\end{aligned}\right)\text{的概率}
&= \iint p(y, \dot{y}) \, \mathrm{d}y \mathrm{d}\dot{y}, \\
&= \int_{0}^{\infty} \mathrm{d}\dot{y} \int_{a-\dot{y}dt}^{a} \mathrm{d}y \, p(y, \dot{y})
\end{aligned}
$$
当时$dt→0，p(y)=p(a)$
\[p(y,\dot{y})=p(y=a,\dot{y})
\]
所以，对dy积分
\[\begin{aligned}
\int_{a-\dot{y}dt}^{a} p(y, \dot{y}) \, \mathrm{d}y &\approx \int_{a-\dot{y}dt}^{a} p(a, \dot{y}) \, \mathrm{d}y \\
&= p(a, \dot{y}) \int_{a-\dot{y}dt}^{a} 1 \, \mathrm{d}y \\
&= p(a, \dot{y}) [a - (a - \dot{y}dt)] \\
&= p(a, \dot{y})\dot{y}dt
\end{aligned}
\]
\[\begin{aligned}
\left( \begin{array}{c} \text{在 } dt \text{ 时间内} \\ \text{发生正穿越 } a \end{array} \right) \text{的概率} &= \int_{0}^{\infty} p(y = a, \dot{y})\dot{y}dt\mathrm{d}\dot{y} \\
&= dt \int_{0}^{\infty} p(a, \dot{y})\dot{y}\mathrm{d}\dot{y}
\end{aligned}
\]
正穿越频率记为\(\nu_a^+\)，在T时刻正穿越a的次数是 \(\nu_a^+ T\), 在dt事件内正穿越a的次数是 $\nu_a^+ dt $。联系起宏观频率和微观概率的一个重要结论：
而在dt事件正穿越a的平均次数等于dt时间正穿越a的概率。
数学上描述：
\[\nu_a^+ dt = dt \int_{0}^{\infty} p(a, \dot{y}) \dot{y} d\dot{y}
\]
正穿越频率：
\[\nu_a^+ = \int_{0}^{\infty} p(a, \dot{y})\dot{y} \, d\dot{y}
\]
这是一般意义的莱斯公式，和是不是高斯过程无关。
尝试用例子来解释上面那个结论：
假设的例子：
假设地铁闸机在极短的 1 秒内最多只能通过 1 个人（这意味着要么过 1 个，要么过 0 个），且闸机不能保证每次都能让人通过。
- 1. 时间维度（看单一个体）：经过统计，1个闸机在 1000秒内通过了 20 个人，那么它每秒通过的频率就是 20/1000=0.0220/1000=0.02 人/秒。
- 2. 空间维度（看集合群体）：如果全世界有 1000 个类似的闸机，你同时观察它们 1 秒钟。因为大家概率一样，平均每个闸机通过的人数依然是 0.02 人（这 1000 个闸机加起来总共过了 20 人）。
因此：
前者，可以说每秒通过闸机的平均次数是 0.02 次/秒；
后者，可以说单台闸机发生通过的概率是 0.02。
（因为在 1 秒内最多只过 1 个人，所以“平均次数”在数值上刚好等于“概率”）。
实际上用了随机过程的遍历性和平稳性。实际中，我们往往会用时间的上的统计量去估计空间的统计量，这里面是反过来的。
5.2.2 窄带高斯过程的莱斯公式
y为窄带高斯过程，因为y和\(\dot{y}\) 独立：
\[p(y,\dot{y})=p(y)p(\dot{y})
\]
有因为p(a) 是常数，所以
\[\begin{aligned}
v_a^+ &= \int_{0}^{\infty}\dot y \cdot p(a,\dot y)\,\mathrm{d}\dot y \\
&=\int_{0}^{\infty}\dot y \cdot p(a)\,p(\dot y)\,\mathrm{d}\dot y \\
&= p(a)\int_{0}^{\infty}\dot y \,p(\dot y)\,\mathrm{d}\dot y
\end{aligned}
\]
物理意义如下：
\[\underbrace{v_a^+}_{\text{总穿越频率}} = \underbrace{p(a)}_{\text{位置概率密度（常数）}} \cdot \underbrace{\int_{0}^{\infty} \underbrace{p(\dot{y})}_{\text{速度概率密度}} \cdot \underbrace{\dot{y}}_{\text{速度}} \, \underbrace{d\dot{y}}_{\text{速度微元}}}_{\text{所有正向速度的加权平均值（期望）}}
\]
两个高斯分布的概率密度如下：
\[p(a)=\frac{1}{\sqrt{2\pi}\,\sigma_y}\exp\left(-\frac{a^2}{2\sigma_y^2}\right)\\
p(\dot y)=\frac{1}{\sqrt{2\pi}\,\sigma_{\dot y}}\exp\left(-\frac{\dot y^2}{2\sigma_{\dot y}^2}\right)
\]
速度积分积分可以算出来（计算过程见附录）
\[\int_{0}^{\infty}\dot y\,p(\dot y)\,\mathrm{d}\dot y=\frac{\sigma_{\dot y}}{\sqrt{2\pi}}
\]
因此，向上穿越率为
\[\begin{aligned}
v_a^+
&={\frac{\sigma_{\dot y}}{2\pi\sigma_y}\,e^{-\frac{a^2}{2\sigma_y^2}}}
\end{aligned}
\]
这是经典的莱斯公式（Rice's formula）在窄带高斯过程的形式。
5.2.3 峰值分布→瑞利分布
xxxx
附录：
推导
Rice向上穿越率完整推导
目标：\(v_a^+=\int_{0}^{+\infty}\dot y\, p(a,\dot y)\,\mathrm{d}\dot y\)
①代入联合密度，常数提出积分外
\[\begin{aligned}
v_a^+ &= \int_{0}^{\infty}\dot y \cdot p(a,\dot y)\,\mathrm{d}\dot y \\
&=\int_{0}^{\infty}\dot y \cdot p(a)\,p(\dot y)\,\mathrm{d}\dot y \\
&= p(a)\int_{0}^{\infty}\dot y \,p(\dot y)\,\mathrm{d}\dot y
\end{aligned}
\]
②边缘密度
\[p(a)=\frac{1}{\sqrt{2\pi}\,\sigma_y}\exp\left(-\frac{a^2}{2\sigma_y^2}\right),\qquad
p(\dot y)=\frac{1}{\sqrt{2\pi}\,\sigma_{\dot y}}\exp\left(-\frac{\dot y^2}{2\sigma_{\dot y}^2}\right)
\]
③计算速度积分 \(I=\int_{0}^{\infty}\dot y\,p(\dot y)\,\mathrm{d}\dot y\)
\[I=\int_{0}^{\infty}\dot y \cdot \frac{1}{\sqrt{2\pi}\,\sigma_{\dot y}}\,e^{-\dot y^2/(2\sigma_{\dot y}^2)}\,\mathrm{d}\dot y
\]
④变量代换：令 \(u=\frac{\dot y^2}{2\sigma_{\dot y}^2}\)
\[\mathrm{d}u=\frac{\dot y}{\sigma_{\dot y}^2}\,\mathrm{d}\dot y \implies \dot y\,\mathrm{d}\dot y=\sigma_{\dot y}^2\,\mathrm{d}u
\]
上下限：\(\dot y=0 \Rightarrow u=0;\ \dot y\to\infty \Rightarrow u\to\infty\)
⑤求解积分
\[\begin{aligned}
I&=\frac{1}{\sqrt{2\pi}\,\sigma_{\dot y}}\int_{0}^{\infty} e^{-u}\cdot \sigma_{\dot y}^2 \,\mathrm{d}u \\
&=\frac{\sigma_{\dot y}}{\sqrt{2\pi}}\int_{0}^{\infty}e^{-u}\,\mathrm{d}u \\
&=\frac{\sigma_{\dot y}}{\sqrt{2\pi}}\Big[-e^{-u}\Big]_0^\infty \\
&=\frac{\sigma_{\dot y}}{\sqrt{2\pi}}
\end{aligned}
\]
⑥合并得到Rice穿越率
\[\begin{aligned}
v_a^+ &= p(a)\cdot I \\
&= \frac{1}{\sqrt{2\pi}\,\sigma_y}\,e^{-\frac{a^2}{2\sigma_y^2}} \cdot \frac{\sigma_{\dot y}}{\sqrt{2\pi}} \\
&=\boldsymbol{\frac{\sigma_{\dot y}}{2\pi\sigma_y}\,e^{-\frac{a^2}{2\sigma_y^2}}}
\end{aligned}
\]
补充：连续分布中 \(p(a)\) 是概率密度函数值，单点概率 \(P(Y=a)=\int_{a}^{a}p(y)\mathrm{d}y=0\)；穿越率是单位时间向上穿过\(a\)的期望次数，不是单点概率。
代码
窄带过程典型样本画图
import matplotlib.pyplot as plt
import numpy as np
# 设置字体
plt.rcParams['font.sans-serif'] = ['Arial']
plt.rcParams['axes.unicode_minus'] = False
fig, ax = plt.subplots(figsize=(12, 4.5))
# ----------------- 1. 生成波形数据 -----------------
t = np.linspace(0, 11, 4000) # 提高采样率让波形更平滑
f = 1.4 # 载波频率
carrier = np.sin(2 * np.pi * f * t)
a = 1.6 # 虚线高度 y=a
# 找出载波所有的波峰时间点 (sin = 1)
peak_indices = np.where(np.diff(np.sign(np.cos(2 * np.pi * f * t))) < 0)[0]
peak_times = t[peak_indices]
# 选取右侧连续的三次波峰时间点
idx_start = int(len(peak_times) * 0.70)
t_p1, t_p2, t_p3 = peak_times[idx_start], peak_times[idx_start+1], peak_times[idx_start+2]
# 构建包络线控制点 (使用 numpy 的 interp 进行线性插值，由于采点密集，视觉上也是平滑的)
# 注意这里我们将 t_p1, t_p2, t_p3 处的包络高度设置为了 1.15 * a，确保波形"穿过去"
control_t = [0, 0.8, 1.5, 2.5, 3.5, 4.5, 5.5, t_p1, t_p2, t_p3, 10, 11]
control_env = [0, 1.3, 0.9, 1.1, 0.4, 0.3, 0.5, 1.15*a, 1.15*a, 1.15*a, 0.4, 0]
envelope = np.interp(t, control_t, control_env)
# 生成最终波形
y = envelope * carrier
# ----------------- 2. 绘图 -----------------
# 水平时间轴 (X轴)
ax.axhline(0, color='black', linewidth=1.5)
# 波形曲线 - 蓝色
ax.plot(t, y, color='#1f77b4', linewidth=1.8)
# 顶部水平虚线 (y = a)
ax.axhline(a, color='black', linewidth=1.5, linestyle='--', xmin=0.02, xmax=0.92)
ax.text(0.2, a + 0.15, r"$y=a$", fontsize=16, color='black')
# 左右两侧的垂直边界虚线
x_left = 2.0
x_right = 9.8
ax.axvline(x_left, color='black', linewidth=1.5, linestyle='--', ymin=0.1, ymax=0.9)
ax.axvline(x_right, color='black', linewidth=1.5, linestyle='--', ymin=0.1, ymax=0.9)
# 底部的双向箭头和周期 T
ax.annotate('', xy=(x_right, -2.2), xytext=(x_left, -2.2),
arrowprops=dict(arrowstyle='<|-|>', color='red', linewidth=2, mutation_scale=15))
ax.text((x_left+x_right)/2, -2.3, "$T$", fontsize=18, ha='center', va='top', color='red', fontweight='bold')
# 坐标轴箭头
# 左侧 Y 轴
ax.annotate('', xy=(0, 2.6), xytext=(0, -0.5),
arrowprops=dict(arrowstyle='-|>', color='black', linewidth=1.8, mutation_scale=20))
ax.text(0.15, 2.4, "$y(t)$", fontsize=16, color='black')
# 右侧 X 轴
ax.annotate('', xy=(11.6, 0), xytext=(10.8, 0),
arrowprops=dict(arrowstyle='-|>', color='black', linewidth=1.8, mutation_scale=20))
ax.text(11.3, -0.4, "$t$", fontsize=16, color='black')
# ----------------- 3. 细节调整 -----------------
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['bottom'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.set_xticks([])
ax.set_yticks([])
ax.set_xlim(-0.5, 12)
ax.set_ylim(-2.8, 2.8)
plt.tight_layout()
plt.show()
dt时间内发生正穿越a的概率计算
import matplotlib.pyplot as plt
import numpy as np
# 设置字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'Arial Unicode MS']
plt.rcParams['axes.unicode_minus'] = False
fig = plt.figure(figsize=(10, 8), dpi=150)
ax = fig.add_subplot(111, projection='3d')
ax.set_axis_off()
# 定义不同性质曲线的颜色
color_surface = '#1f77b4' # 曲面轮廓线：蓝色
color_axis = 'black' # 坐标轴：黑色
color_boundary = '#d62728' # 边界线 (y=a 与斜线)：红色
color_shadow = '#7f7f7f' # 阴影区域：灰色
color_intersect = '#6a0dad' # 曲面截交线：深紫色
# ================= 1. 设定视角 =================
ax.view_init(elev=25, azim=-45)
# ================= 2. 绘制挺拔的二维正态分布曲面 =================
peak_z = 4.2 # 峰值高度
sigma = 2.2 # 扩散范围
# 2.1 水平闭环圆
theta = np.linspace(0, 2*np.pi, 100)
levels = [ 0.3, 1.8, ]
for level in levels:
r = sigma * np.sqrt(-2 * np.log(level / peak_z))
x_c = r * np.cos(theta)
y_c = r * np.sin(theta)
z_c = np.full_like(x_c, level)
ax.plot(x_c, y_c, z_c, color=color_surface, lw=1)
# 2.2 交叉截线
t = np.linspace(-8.0, 8.0, 100)
z_curve = peak_z * np.exp(-t**2 / (2 * sigma**2))
ax.plot(t, np.zeros_like(t), z_curve, color=color_surface, lw=1)
ax.plot(np.zeros_like(t), t, z_curve, color=color_surface, lw=1)
# ================= 3. 绘制坐标轴 =================
# p(y, y') 垂直轴：加长到 6.5
ax.quiver(0, 0, 0, 0, 0, 6.5, color=color_axis, lw=1.5, arrow_length_ratio=0.1)
ax.text(0, 0, 6.8, r"$p(y, \dot{y})$", fontsize=16, fontfamily='serif', ha='center')
# y 轴：箭头指向负 y 方向，长度保持 8.0
ax.quiver(0, 0, 0, 0, -8.0, 0, color=color_axis, lw=1, arrow_length_ratio=0.1)
ax.plot([0, 0], [0, 8.0], [0, 0], color=color_axis, lw=1)
ax.text(0, -8.5, 0, r"$y$", fontsize=16, fontfamily='serif')
# y' 轴：长度保持 8.0
ax.quiver(0, 0, 0, 8.0, 0, 0, color=color_axis, lw=1, arrow_length_ratio=0.1)
ax.plot([-8.0, 0], [0, 0], [0, 0], color=color_axis, lw=1)
ax.text(8.5, 0, 0, r"$\dot{y}$", fontsize=16, fontfamily='serif')
# ================= 4. 绘制底面的两条边界线与阴影 =================
a_val = -4.0
# 4.1 平行于 y' 轴的直线 y = a
x_line_a = np.linspace(0, 6.0, 50)
y_line_a = np.full_like(x_line_a, a_val)
z_line_a = np.zeros_like(x_line_a)
ax.plot(x_line_a, y_line_a, z_line_a, color=color_boundary, lw=2.5)
# 4.2 斜线
slope = 0.15
y_line_slope = a_val + slope * x_line_a
ax.plot(x_line_a, y_line_slope, z_line_a, color=color_boundary, lw=2.5, ls='--')
# 4.3 阴影区域
for x_val in np.linspace(0, 4.0, 30):
y_start = a_val
y_end = a_val + slope * x_val
ax.plot([x_val, x_val], [y_start, y_end], [0.01, 0.01], color=color_shadow, lw=1)
# ================= 5. 绘制曲面上的截交线 =================
# 5.1 平面 y=a 的截交线
x_cross_a = np.linspace(0, 6.0, 100)
y_cross_a = np.full_like(x_cross_a, a_val)
z_cross_a = peak_z * np.exp(-(x_cross_a**2 + a_val**2) / (2 * sigma**2))
ax.plot(x_cross_a, y_cross_a, z_cross_a, color=color_intersect, lw=2.5)
ax.plot([0, 0], [a_val, a_val], [0, z_cross_a[0]], color=color_intersect, lw=2.5)
ax.plot([6, 6], [a_val, a_val], [z_cross_a[-1], 0], color=color_intersect, lw=2.5)
# 5.2 斜线平面的截交线
x_cross_slope = np.linspace(0, 6.0, 100)
y_cross_slope = a_val + slope * x_cross_slope
z_cross_slope = peak_z * np.exp(-(x_cross_slope**2 + y_cross_slope**2) / (2 * sigma**2))
ax.plot(x_cross_slope, y_cross_slope, z_cross_slope, color=color_intersect, lw=2.5, ls='--')
ax.plot([0, 0], [a_val, a_val], [0, z_cross_slope[0]], color=color_intersect, lw=2.5, ls='--')
ax.plot([6, 6], [y_cross_slope[-1], y_cross_slope[-1]], [z_cross_slope[-1], 0], color=color_intersect, lw=2.5, ls='--')
# ================= 6. 标注 =================
ax.scatter(0, a_val, 0, color=color_boundary, s=25)
ax.text(0, a_val - 0.5, 0, r"$a$", fontsize=16, fontfamily='serif')
ax.text(4, a_val - 1.5, 0, r"$y=a$", fontsize=14, fontfamily='serif')
ax.text(0.5, a_val + 0.4, 0, r"$\alpha$", fontsize=16, fontfamily='serif')
ax.text2D(0.05, 0.9, r"$(a)$", transform=ax.transAxes, fontsize=18, fontfamily='serif')
# 坐标轴显示范围 (z 轴拉长，需同步增大 zlim)
ax.set_xlim([-4.0, 4.0])
ax.set_ylim([-4.0, 4.0])
ax.set_zlim([0, 6.5])
plt.tight_layout()
plt.show()
dt时间内穿越概率积分区间
import matplotlib.pyplot as plt
import matplotlib.patches as patches
# 设置支持数学公式的字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['mathtext.fontset'] = 'stix'
# 创建画布
fig, ax = plt.subplots(figsize=(8, 6))
# 隐藏坐标轴边框和刻度
ax.set_axis_off()
# 设置坐标系范围
ax.set_xlim(-5, 10)
ax.set_ylim(-6, 6)
# ==========================================
# 1. 绘制坐标轴
# ==========================================
# 上方横轴 (带箭头)
ax.annotate('', xy=(9, 3), xytext=(-4, 3),
arrowprops=dict(arrowstyle='-|>', color='black', lw=2, mutation_scale=20))
# ax.plot([-4, 1], [2, 2], color='black', lw=3) # 左侧加粗
# 纵轴 (带箭头)
ax.annotate('', xy=(1, -5), xytext=(1, 5),
arrowprops=dict(arrowstyle='-|>', color='black', lw=2, mutation_scale=20))
# ax.plot([1, 1], [5, 2], color='black', lw=3) # 上方加粗
# 下方平行线 (不带箭头，在 x=8 处结束)
ax.plot([-4, 8], [-1, -1], color='black', lw=3)
# ==========================================
# 2. 绘制阴影区域
# ==========================================
# 顶点坐标：紧贴纵轴 (1, -1)，向右延伸至 (8, -1)，向上至 (8, 2)
triangle_vertices = [(1, -1), (8, -1), (8, 2)]
triangle = patches.Polygon(triangle_vertices, closed=True, facecolor='white', edgecolor='black', hatch='//', lw=1)
ax.add_patch(triangle)
# ==========================================
# 3. 添加文本标注
# ==========================================
# 原点附近的 'a'
ax.text(0.7, -1.3, 'a', fontsize=14, fontweight='bold', fontfamily='serif')
# 上方横轴最右侧的 'y'
ax.text(8.5, 2.5, 'y', fontsize=16, fontweight='bold', fontfamily='serif')
# 纵轴最上方箭头的 'y' (带点)
# ax.text(1.2, 4.8, r'$\dot{y}$', fontsize=18, fontweight='bold', fontfamily='serif')
# 纵轴最下方箭头的 'y' (带点)
ax.text(1.2, -5.2, r'$\dot{y}$', fontsize=18, fontweight='bold', fontfamily='serif')
# 【修改点】：已删除下方平行线上的标识代码
# 阴影区域内部左侧的 'α'
ax.text(4.2, -0.4, r'$\alpha$', fontsize=14, fontfamily='serif')
# 阴影区域内部右侧的 '(a-y)'
ax.text(6.8, -0.4, 'a-y', fontsize=12, fontweight='bold', fontfamily='serif')
# 在 '(a-y)' 左侧绘制一条粗短竖线，表示线段长度标记
# ax.plot([6.4, 6.4], [-0.8, -0.2], color='black', lw=3)
# 显示图像
plt.tight_layout()
plt.show()
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import make_interp_spline
# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
fig, ax = plt.subplots(figsize=(10, 5))
# ----------------- 1. 绘制坐标轴 -----------------
ax.annotate('', xy=(1.5, 9.0), xytext=(1.5, 1.0),
arrowprops=dict(arrowstyle='-|>', color='black', linewidth=2, mutation_scale=20))
ax.annotate('', xy=(9.5, 1.0), xytext=(1.0, 1.0),
arrowprops=dict(arrowstyle='-|>', color='black', linewidth=2, mutation_scale=20))
ax.text(1.7, 8.8, r'$y(t)$', fontsize=15, fontweight='bold')
ax.text(9.3, 0.6, r'$t$', fontsize=15, fontweight='bold')
# ----------------- 2. 绘制水平虚线 y=a -----------------
y_a = 6.2
ax.plot([1.5, 8.8], [y_a, y_a], color='black', linewidth=2, linestyle='--', dashes=(8, 4))
ax.text(2.5, y_a + 0.3, r'$y=a$', fontsize=14, fontweight='bold')
# ----------------- 3. 样条曲线控制点 -----------------
x_points = np.array([3.8, 4.5, 5.2, 6.0, 6.8, 7.5, 8.2])
y_points = np.array([2.0, 3.0, 4.2, 5.5, 6.5, 7.2, 7.6])
spline = make_interp_spline(x_points, y_points, k=3)
x_curve = np.linspace(3.8, 8.2, 200)
y_curve = spline(x_curve)
ax.plot(x_curve, y_curve, color='black', linewidth=2.5)
# ----------------- 4. 绘制割线及辅助线 -----------------
t_x = 5.0
t_y = spline(t_x)
x_cross = 7.3
# 计算割线斜率
slope = (y_a - t_y) / (x_cross - t_x)
def secant_line(x):
return slope * (x - t_x) + t_y
# 绘制割线
x_secant = np.linspace(4.2, 8.6, 100)
y_secant = secant_line(x_secant)
ax.plot(x_secant, y_secant, color='black', linewidth=2, linestyle='--', dashes=(6, 4))
# ----------------- 5. 绘制垂直辅助线 -----------------
# t 处的垂直实线
ax.plot([t_x, t_x], [1.0, t_y + 2.5], color='black', linewidth=1.8)
# t+dt 处的垂直虚线
ax.plot([x_cross, x_cross], [1.0, y_a + 1.5], color='black', linewidth=1.8, linestyle='-.')
# 交点处的短横线
ax.plot([x_cross - 0.15, x_cross + 0.15], [y_a, y_a], color='black', linewidth=2)
# ----------------- 6. 添加文本标注 -----------------
ax.text(t_x, 0.5, r'$t$', fontsize=15, ha='center', fontweight='bold')
ax.text(x_cross, 0.5, r'$t+dt$', fontsize=15, ha='center', fontweight='bold')
# 曲线上的 y 标注
ax.text(t_x - 0.35, t_y + 0.1, r'$y$', fontsize=15, fontweight='bold')
# 右上角的文字说明 (已向左移动，避开曲线)
ax.text(5.6, 8.2, r'$dt$内正穿越所', fontsize=13, fontweight='bold')
ax.text(5.6, 7.5, r'需的最小斜率', fontsize=13, fontweight='bold')
# ----------------- 7. 细节调整 -----------------
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['bottom'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.set_xticks([])
ax.set_yticks([])
ax.set_xlim(0, 10)
ax.set_ylim(0, 9)
plt.tight_layout()
plt.show()
随机过程穿越分析：宏观集合与微元概率
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'Arial Unicode MS', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
rng = np.random.default_rng(2023)
# ============================================================
# 生成 500 条独立的窄带高斯平稳随机过程
# ============================================================
fs = 1000.0
T_total = 10.0
t = np.arange(0, T_total, 1/fs)
f0 = 5.0
fc = 0.5
n_samples = 500
a = 0.5 # 阈值
def generate_narrowband(rng):
n = len(t)
white = rng.standard_normal(n)
F = np.fft.rfft(white)
freqs = np.fft.rfftfreq(n, d=t[1]-t[0])
F[freqs > fc] = 0
slow = np.fft.irfft(F, n)
slow = slow / np.std(slow)
env = np.abs(slow) + 0.3
phase = 2*np.pi*f0*t + 0.8*slow
return env * np.cos(phase)
Y = np.array([generate_narrowband(rng) for _ in range(n_samples)])
def find_upcrossings(t, y, a):
idx = np.where((y[:-1] < a) & (y[1:] >= a))[0]
return idx
# 微元窗口参数
t0 = 2.0
dt = 0.05
# ============================================================
# 绘图设置：2 行 1 列，取消多余的白边
# ============================================================
fig = plt.figure(figsize=(15, 10))
gs = fig.add_gridspec(2, 1, height_ratios=[1, 1], hspace=0.3)
# ---------- 1. 宏观视角（上） ----------
ax1 = fig.add_subplot(gs[0, 0])
for i in range(20):
ax1.plot(t, Y[i], color='steelblue', alpha=0.3, lw=0.8)
ax1.plot(t, Y[0], color='navy', lw=1.5, label='典型单样本')
ax1.axhline(a, color='crimson', ls='--', lw=1.5, label=f'阈值 a={a}')
ax1.axvspan(t0, t0+dt, color='darkorange', alpha=0.3, label=f'微元窗口 dt={dt}s')
ax1.set_xlim(1.5, 3.0)
ax1.set_ylim(-2.5, 2.5)
ax1.set_xlabel(r'时间 \(t\) (s)')
ax1.set_ylabel(r'幅值 \(y(t)\)')
ax1.set_title(r'【宏观视角】随机过程样本集合与微元锁定', fontsize=14)
ax1.legend(loc='upper right')
ax1.grid(alpha=0.3)
# ---------- 2. 微元切片（下） ----------
ax2 = fig.add_subplot(gs[1, 0])
window_view = (t >= t0 - 0.02) & (t <= t0 + dt + 0.02)
for i in range(n_samples):
y_win = Y[i, window_view]
t_view = t[window_view]
ups = find_upcrossings(t_view, y_win, a)
if len(ups) > 0:
ax2.plot(t_view, y_win, color='crimson', alpha=0.4, lw=1.0)
else:
ax2.plot(t_view, y_win, color='lightgray', alpha=0.3, lw=0.8)
ax2.axhline(a, color='crimson', ls='--', lw=1.5)
ax2.axvline(t0, color='k', ls=':', lw=1.2)
ax2.axvline(t0+dt, color='k', ls=':', lw=1.2)
ax2.set_title(r'【微元切片】红线代表发生正穿越 (占比 = 概率)', fontsize=14)
ax2.set_xlabel(r'时间 \(t\) (s)')
ax2.set_ylabel(r'幅值 \(y(t)\)')
ax2.grid(alpha=0.3)
fig.suptitle('随机过程穿越分析：宏观集合与微元概率', fontsize=16, y=0.95)
plt.savefig('two_panel_visual.png', dpi=180, bbox_inches='tight')
plt.show()
浙公网安备 33010602011771号