频域疲劳（五）穿越分析-补充完整版

频域疲劳（五）：穿越分析

穿越分析的基本结论

正穿越 \(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 穿越问题

无论是时域疲劳还是频域疲劳，都可以通过求解幅值分布来计算，不同的是计算方式。

而幅值分布可以看成如下图所示，求取峰值的分布。为了进行一般性的推导，我们先求解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}\]

上述条件可以视为y和\(\dot{y}\)的联合分布中如图楔形部分所示的体积。

积分区域：

5.2 积分计算

5.2.1 穿越概率的积分表达式

设 \(p(y,\dot{y})\) 为位移 \(y\) 与速度 \(\dot{y}\) 的联合概率密度。由 5.1 节的条件 i、ii，在 \(dt\) 时间内发生一次正穿越 \(y=a\) 的概率，等于联合密度在楔形区域上的体积：

\[P_a^+(dt)=\int_0^\infty\int_{a-\dot{y}dt}^{a}p(y,\dot{y})\,dy\,d\dot{y} \]

内层对 \(y\) 积分，外层对 \(\dot{y}\) 积分。内层积分区间宽度为 \(\dot{y}dt\)，当 \(dt\to 0\) 时趋于无穷窄。由积分中值定理（见附录推导B），对任意连续的 \(p\)，有：

\[\int_{a-\dot{y}dt}^{a}p(y,\dot{y})\,dy\;\approx\; p(a,\dot{y})\,\dot{y}dt \]

即内层积分只取决于 \(y=a\) 处的密度值，与区间内部的细节无关。代入得：

\[P_a^+(dt)=dt\int_0^\infty p(a,\dot{y})\,\dot{y}\,d\dot{y} \]

因此单位时间内的平均正穿越频率为：

\[\nu_a^+=\lim_{dt\to 0}\frac{P_a^+(dt)}{dt}=\int_0^\infty p(a,\dot{y})\,\dot{y}\,d\dot{y} \]

这就是文章开头给出的一般公式：穿越频率等于联合密度在 \(y=a\) 处沿速度方向的一阶矩。该公式对任何平稳随机过程都成立，接下来的工作只是把具体的联合密度代进去。

5.2.2 高斯过程的联合密度

对于平稳高斯过程 \(y(t)\)，位移 \(y\) 与速度 \(\dot{y}\) 不相关（证明见附录推导A），而对高斯过程而言不相关等价于独立，因此联合密度可分离：

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

代入一般公式：

\[\nu_a^+=\frac{1}{2\pi\sigma_y\sigma_{\dot{y}}}\int_0^\infty \dot{y}\, \exp\left[-\frac{a^2}{2\sigma_y^2}-\frac{\dot{y}^2}{2\sigma_{\dot{y}}^2}\right]d\dot{y} \]

其中与 \(\dot{y}\) 无关的因子 \(e^{-a^2/2\sigma_y^2}\) 可以提到积分号外：

\[\nu_a^+=\frac{1}{2\pi\sigma_y\sigma_{\dot{y}}}\,e^{-a^2/2\sigma_y^2} \underbrace{\int_0^\infty \dot{y}\,e^{-\dot{y}^2/2\sigma_{\dot{y}}^2}d\dot{y}}_{I} \]

5.2.3 关键积分 \(I\) 的计算

令 \(u=\dot{y}^2/(2\sigma_{\dot{y}}^2)\)，则 \(\dot{y}\,d\dot{y}=\sigma_{\dot{y}}^2\,du\)（推导见附录C），于是：

\[I=\int_0^\infty \dot{y}\,e^{-\dot{y}^2/2\sigma_{\dot{y}}^2}d\dot{y} =\sigma_{\dot{y}}^2\int_0^\infty e^{-u}\,du =\sigma_{\dot{y}}^2 \]

5.2.4 Rice 公式

代回，得 \(dt\) 时间内正穿越 \(y=a\) 的概率：

\[P_a^+(dt)=\frac{dt}{2\pi}\frac{\sigma_{\dot{y}}}{\sigma_y}\,e^{-a^2/2\sigma_y^2} \]

以及单位时间内的平均正穿越频率：

\[\boxed{\;\nu_a^+=\frac{1}{2\pi}\frac{\sigma_{\dot{y}}}{\sigma_y}\,e^{-a^2/2\sigma_y^2}\;} \]

这就是著名的 Rice 公式。它表明：对高斯过程，正穿越某一水平 \(a\) 的平均频率等于过程平均频率 \(\nu_0^+\) 乘以衰减因子 \(e^{-a^2/2\sigma_y^2}\)——水平 \(a\) 离均值越远，穿越越稀少。

5.2.5 特例：正穿越零电平

取 \(a=0\)，得到过程正穿越零电平的平均频率，即过程的"平均频率"：

\[\nu_0^+=\frac{1}{2\pi}\frac{\sigma_{\dot{y}}}{\sigma_y} \]

5.2.6 用功率谱密度表示

穿越频率只依赖两个参数 \(\sigma_y\) 与 \(\sigma_{\dot{y}}\)，它们都可以由功率谱密度直接算出。引入谱矩（单边谱定义，见附录推导D）：

\[m_k=\int_0^\infty \omega^k S_{yy}(\omega)\,d\omega \]

其中 \(S_{yy}(\omega)\) 为单边功率谱密度，\(m_0\) 与 \(m_2\) 分别等于位移方差与速度方差：

\[\sigma_y^2=m_0,\qquad \sigma_{\dot{y}}^2=m_2 \]

因此：

\[\nu_0^+=\frac{1}{2\pi}\sqrt{\frac{m_2}{m_0}},\qquad \nu_a^+=\nu_0^+\exp\left(-\frac{a^2}{2m_0}\right) \]

若用频率 \(f\)（Hz）表示、\(G_{yy}(f)\) 为相应的单边谱，由 \(\omega=2\pi f\) 得 \(m_k=\int_0^\infty (2\pi f)^k G_{yy}(f)\,df\)，于是：

\[\nu_0^+=\frac{1}{2\pi}\sqrt{\frac{m_2}{m_0}} =\sqrt{\frac{\int_0^\infty f^2 G_{yy}(f)\,df}{\int_0^\infty G_{yy}(f)\,df}} \]

即穿零频率等于功率谱的均方根频率（谱"重心"），这是工程上常用的形式。

5.2.7 窄带过程的平均穿越频率

对于窄带高斯过程，功率谱集中在中心频率 \(\omega_0\) 附近且带宽 \(\Delta\omega\ll\omega_0\)，此时 \(m_2\approx\omega_0^2 m_0\)（见附录推导E），故：

\[\nu_0^+=\frac{1}{2\pi}\sqrt{\frac{m_2}{m_0}}\approx\frac{\omega_0}{2\pi}=f_0 \]

即窄带过程的平均正穿越频率约等于其中心频率。这与"窄带过程可视为频率 \(f_0\) 的准正弦波，每个周期正穿越零一次"的直观完全一致，也说明宽带程度越高，穿零频率偏离中心频率越大。

得到穿越频率之后，可以进一步讨论峰值与包络的分布——在窄带近似下，包络峰值超过 \(a\) 的概率与穿越频率密切相关，这是下一节的内容。

附录：

推导

推导A：平稳高斯过程中 \(y\) 与 \(\dot{y}\) 不相关（进而独立）

对平稳过程 \(y(t)\)，自相关函数 \(R_y(\tau)=E[y(t)y(t+\tau)]\)。位移与速度的互相关为：

\[E[y(t)\dot{y}(t+\tau)]=\frac{\partial}{\partial\tau}E[y(t)y(t+\tau)]=R_y'(\tau) \]

令 \(\tau=0\)：

\[E[y\dot{y}]=R_y'(0) \]

由平稳性，\(R_y(\tau)\) 是偶函数，即 \(R_y(\tau)=R_y(-\tau)\)，其导数 \(R_y'(\tau)\) 为奇函数，故 \(R_y'(0)=0\)，于是：

\[E[y\dot{y}]=0 \]

即同一时刻的位移与速度不相关。

从频域看更直观：\(R_y(\tau)=\int_{-\infty}^{\infty}S_y(\omega)e^{i\omega\tau}d\omega\)，则 \(R_y'(0)=i\int_{-\infty}^{\infty}\omega S_y(\omega)d\omega\)；对实值平稳过程，双边谱 \(S_y(\omega)\) 是 \(\omega\) 的偶函数，被积函数 \(\omega S_y(\omega)\) 是奇函数，积分亦为零。

物理上，"平稳"意味着过程上升与下降的机会均等，位移与速度的乘积在统计意义上平均为零。对联合高斯的 \((y,\dot{y})\)，不相关即独立，故联合密度可分离为边缘密度之积。

推导B：内层积分的近似（积分中值定理）

考察内层积分

\[\int_{a-\dot{y}dt}^{a}f(y)\,dy,\qquad f(y)=p(y,\dot{y}) \]

由积分中值定理，存在 \(\xi\in(a-\dot{y}dt,\,a)\) 使得

\[\int_{a-\dot{y}dt}^{a}f(y)\,dy=f(\xi)\,\dot{y}dt \]

当 \(dt\to 0\) 时 \(\xi\to a\)，故

\[\int_{a-\dot{y}dt}^{a}p(y,\dot{y})\,dy\;\xrightarrow{\;dt\to 0\;}\;p(a,\dot{y})\,\dot{y}dt \]

等价地，也可将 \(f(y)\) 在 \(y=a\) 处一阶泰勒展开，\(f(y)=f(a)+f'(\eta)(y-a)\)，代入积分后第二项为 \(O\big((\dot{y}dt)^2\big)\) 量级，可忽略。注意 \(\dot{y}>0\) 时 \(a-\dot{y}dt<a\)，积分区间确为以 \(a\) 为右端点、长度为 \(\dot{y}dt\) 的左邻域，这与"从下方逼近 \(a\)"的穿越几何一致。

推导C：关键积分 \(I=\int_0^\infty \dot{y}\,e^{-\dot{y}^2/2\sigma^2}d\dot{y}=\sigma^2\)

作换元 \(u=\dot{y}^2/(2\sigma^2)\)，则

\[du=\frac{\dot{y}}{\sigma^2}d\dot{y}\quad\Longrightarrow\quad \dot{y}\,d\dot{y}=\sigma^2 du \]

当 \(\dot{y}:0\to\infty\) 时 \(u:0\to\infty\)，故

\[I=\sigma^2\int_0^\infty e^{-u}\,du=\sigma^2\Big[-e^{-u}\Big]_0^\infty=\sigma^2 \]

（被积函数在 \(\dot{y}\to\infty\) 时按指数衰减、在 \(\dot{y}=0\) 处取零值，广义积分收敛。）

推导D：谱矩与方差的关系

设 \(S_{yy}(\omega)\) 为单边功率谱密度（定义域 \(\omega\ge 0\)），谱矩定义为

\[m_k=\int_0^\infty \omega^k S_{yy}(\omega)\,d\omega \]

位移方差：

\[\sigma_y^2=E[y^2]=R_y(0)=\int_0^\infty S_{yy}(\omega)\,d\omega=m_0 \]

速度方差： 先求速度的自相关函数。由

\[R_{\dot{y}}(\tau)=E[\dot{y}(t)\dot{y}(t+\tau)]=-\frac{d^2}{d\tau^2}R_y(\tau) \]

（对平稳过程，速度自相关等于位移自相关的二阶负导，可由 \(E[y(t)\dot{y}(t+\tau)]=R_y'(\tau)\) 再对 \(\tau\) 求导得到），于是

\[\sigma_{\dot{y}}^2=E[\dot{y}^2]=R_{\dot{y}}(0)=-R_y''(0) =\int_0^\infty \omega^2 S_{yy}(\omega)\,d\omega=m_2 \]

综上：

\[\frac{\sigma_{\dot{y}}}{\sigma_y}=\sqrt{\frac{m_2}{m_0}} \]

推导E：窄带过程的平均频率

窄带过程的功率谱集中在中心频率 \(\omega_0\) 附近，带宽 \(\Delta\omega\ll\omega_0\)，则

\[m_2=\int_{\omega_0-\Delta\omega}^{\omega_0+\Delta\omega}\omega^2 S_{yy}(\omega)\,d\omega \approx\omega_0^2\int_{\omega_0-\Delta\omega}^{\omega_0+\Delta\omega}S_{yy}(\omega)\,d\omega =\omega_0^2 m_0 \]

故

\[\nu_0^+=\frac{1}{2\pi}\sqrt{\frac{m_2}{m_0}}\approx\frac{\omega_0}{2\pi}=f_0 \]

即窄带过程的平均正穿越频率约等于其中心频率 \(f_0\)。

代码

窄带过程典型样本画图

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()