![image]()
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
from scipy import signal
from scipy.special import gamma as gamma_func
plt.rcParams['font.sans-serif'] = ['SimHei', 'Arial Unicode MS', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
def trapz_compat(y, x):
return np.trapezoid(y, x) if hasattr(np, 'trapezoid') else np.trapz(y, x)
# ============================================================
# ① 生成原始信号 (蒙特卡洛)
# ============================================================
np.random.seed(2024)
fs, T = 500, 200
N = fs*T
t = np.arange(N)/fs
f_axis = np.fft.rfftfreq(N, 1/fs)
# 双峰 PSD
S_target = (10*np.exp(-(f_axis-8)**2/(2*1.2**2)) +
4 *np.exp(-(f_axis-25)**2/(2*2.0**2)) +
0.08)
S_target[0] = 0
ph = np.random.uniform(0, 2*np.pi, len(f_axis))
x = np.fft.irfft(np.sqrt(2*S_target*fs/N)*np.exp(1j*ph), n=N)
# ============================================================
# ② PSD & ③ 谱矩
# ============================================================
f, Pxx = signal.welch(x, fs=fs, nperseg=4096, scaling='density')
m = {n: trapz_compat(f**n * Pxx, f) for n in (0,1,2,4)}
m0, m1, m2, m4 = m[0], m[1], m[2], m[4]
sigma = np.sqrt(m0)
gamma = m2/np.sqrt(m0*m4)
# ============================================================
# ④ 穿越统计
# ============================================================
N0 = (1/(2*np.pi))*np.sqrt(m2/m0) # 零穿越率
Np = (1/(2*np.pi))*np.sqrt(m4/m2) # 峰值率
def crossing_rate(a, m0, m2):
return (1/(2*np.pi))*np.sqrt(m2/m0)*np.exp(-a**2/(2*m0))
# 检测 +1σ 穿越点
thr = sigma
cross_idx = np.where((x[:-1] < thr) & (x[1:] >= thr))[0]
# ============================================================
# ⑤ Dirlik 三分量
# ============================================================
xm = (m1/m0)*np.sqrt(m2/m4)
D1 = 2*(xm - gamma**2)/(1 + gamma**2)
R = (gamma - xm - D1**2)/(1 - gamma - D1 + D1**2)
D2 = (1 - gamma - D1 + D1**2)/(1 - R)
D3 = 1 - D1 - D2
Q = 1.25*(gamma - D3 - D2*R)/D1
def c1(S): Z = S/(2*sigma); return D1/(2*sigma*Q)*np.exp(-Z/Q)
def c2(S): Z = S/(2*sigma); return D2*Z/(2*sigma*R**2)*np.exp(-Z**2/(2*R**2))
def c3(S): Z = S/(2*sigma); return D3*Z/(2*sigma)*np.exp(-Z**2/2)
def dirlik(S): return c1(S)+c2(S)+c3(S)
def rayleigh(S): return S/m0*np.exp(-S**2/(2*m0))
# ============================================================
# ⑥ 疲劳损伤
# ============================================================
b_sn, C_sn = 3.0, 1.0
ESb_ray = (2*m0)**(b_sn/2) * gamma_func(1+b_sn/2)
D_ray = N0*T/C_sn * ESb_ray
ESb_dir = trapz_compat(np.linspace(0, 8*sigma, 800)**b_sn *
dirlik(np.linspace(0, 8*sigma, 800)),
np.linspace(0, 8*sigma, 800))
D_dir = N0*T/C_sn * ESb_dir
# ============================================================
# 可视化: 6 个子图 + 参数面板
# ============================================================
fig = plt.figure(figsize=(18, 11))
gs = GridSpec(3, 3, figure=fig, hspace=0.45, wspace=0.32)
# ---- ① 原始时域 + 穿越点 ----
ax1 = fig.add_subplot(gs[0, :2])
ax1.plot(t, x, 'b-', lw=0.6, alpha=0.8, label='随机载荷 x(t)')
ax1.axhline(+thr, color='r', ls='--', lw=1.2, label=f'+1σ = {thr:.2f}')
ax1.axhline(0, color='k', ls='-', lw=0.8, alpha=0.6)
ax1.axhline(-thr, color='g', ls='--', lw=1.2, label=f'-1σ')
ax1.plot(t[cross_idx], x[cross_idx], 'ro', ms=5,
label=f'+1σ 上穿 {len(cross_idx)} 次')
ax1.set_title('① 原始时域信号与穿越点检测', fontsize=13, fontweight='bold')
ax1.set_xlabel('时间 t (s)'); ax1.set_ylabel('幅值')
ax1.legend(fontsize=9, ncol=4)
ax1.grid(alpha=0.3); ax1.set_xlim(0, 15)
# ---- ⑥ 参数面板 (右上) ----
ax_stat = fig.add_subplot(gs[0, 2]); ax_stat.axis('off')
stats_txt = (
f'谱矩:\n'
f' m0 = {m0:.4e}\n'
f' m1 = {m1:.4e}\n'
f' m2 = {m2:.4e}\n'
f' m4 = {m4:.4e}\n'
f'\n谱参数:\n'
f' σ = {sigma:.4f} γ = {gamma:.4f}\n'
f'\n穿越率:\n'
f' N0 = {N0:.3f} /s\n'
f' Np = {Np:.3f} /s\n'
f' N(+1σ) = {crossing_rate(sigma,m0,m2):.3f} /s\n'
f' 实测上穿 = {len(cross_idx)/T:.3f} /s'
)
ax_stat.text(0.02, 0.98, stats_txt, transform=ax_stat.transAxes,
fontsize=10, family='monospace', va='top',
bbox=dict(boxstyle='round,pad=0.6', facecolor='#FFF8DC',
edgecolor='gray', alpha=0.9))
ax_stat.set_title('⑥ 统计参数汇总', fontsize=12, fontweight='bold')
# ---- ② PSD ----
ax2 = fig.add_subplot(gs[1, 0])
ax2.semilogy(f, Pxx, 'b-', lw=1.3)
ax2.fill_between(f, Pxx, alpha=0.25, color='b')
ax2.set_title('② 功率谱密度 PSD', fontsize=12, fontweight='bold')
ax2.set_xlabel('频率 (Hz)'); ax2.set_ylabel('PSD')
ax2.grid(alpha=0.3, which='both'); ax2.set_xlim(0, 50)
# ---- ③ 谱矩 ----
ax3 = fig.add_subplot(gs[1, 1])
names = ['$m_0$','$m_1$','$m_2$','$m_4$']
vals = [m0, m1, m2, m4]
norm = np.array(vals)/max(vals)
bars = ax3.bar(names, norm, color=['#1f77b4','#ff7f0e','#2ca02c','#d62728'])
for bar, v in zip(bars, vals):
ax3.text(bar.get_x()+bar.get_width()/2, bar.get_height()+0.02,
f'{v:.2e}', ha='center', fontsize=8)
ax3.set_title('③ 谱矩 (归一化)', fontsize=12, fontweight='bold')
ax3.set_ylabel('归一化值'); ax3.grid(alpha=0.3, axis='y')
# ---- ④ 阈值 vs 穿越率 ----
ax4 = fig.add_subplot(gs[1, 2])
a_vals = np.linspace(0.1*sigma, 4*sigma, 100)
Na_theory = crossing_rate(a_vals, m0, m2)
Na_meas = [np.sum((x[:-1]<a)&(x[1:]>=a))/T for a in a_vals]
ax4.semilogy(a_vals/sigma, Na_theory, 'b-', lw=2, label='Rice 理论')
ax4.semilogy(a_vals/sigma, Na_meas, 'ro', ms=3, label='实测')
ax4.axvline(1.0, color='gray', ls=':', lw=1)
ax4.set_title('④ 阈值-穿越率 (Rice公式)', fontsize=12, fontweight='bold')
ax4.set_xlabel('归一化阈值 a/σ'); ax4.set_ylabel('穿越率 (次/s)')
ax4.legend(fontsize=9); ax4.grid(alpha=0.3, which='both')
# ---- ⑤ 幅值分布 (Dirlik 三分量) ----
ax5 = fig.add_subplot(gs[2, :2])
S_g = np.linspace(1e-4, 5*sigma, 600)
ax5.plot(S_g, c1(S_g), color='#1f77b4', lw=1.8, ls='--',
label=f'分量1 (指数, D1={D1:.3f})')
ax5.plot(S_g, c2(S_g), color='#ff7f0e', lw=1.8, ls='--',
label=f'分量2 (Rayleigh, D2={D2:.3f})')
ax5.plot(S_g, c3(S_g), color='#2ca02c', lw=1.8, ls='--',
label=f'分量3 (标准Rayleigh, D3={D3:.3f})')
ax5.plot(S_g, dirlik(S_g), 'k-', lw=2.5, label='Dirlik 总密度')
ax5.plot(S_g, rayleigh(S_g), 'r:', lw=2, label='标准 Rayleigh')
ax5.fill_between(S_g, c1(S_g), alpha=0.1, color='#1f77b4')
ax5.fill_between(S_g, c2(S_g), alpha=0.1, color='#ff7f0e')
ax5.fill_between(S_g, c3(S_g), alpha=0.1, color='#2ca02c')
ax5.set_title('⑤ 幅值分布: Dirlik 三分量与总密度', fontsize=13, fontweight='bold')
ax5.set_xlabel('幅值 S'); ax5.set_ylabel('概率密度 p(S)')
ax5.legend(fontsize=9, ncol=2); ax5.grid(alpha=0.3)
ax5.set_xlim(0, 4*sigma)
# ---- ⑥ 损伤累积 ----
ax6 = fig.add_subplot(gs[2, 2])
S_int = np.linspace(1e-3, 8*sigma, 800)
# 损伤密度: dD/dS = N0·T/C · S^b · p(S)
dmg_ray = N0*T/C_sn * S_int**b_sn * rayleigh(S_int)
dmg_dir = N0*T/C_sn * S_int**b_sn * dirlik(S_int)
ax6.fill_between(S_int/sigma, 0, dmg_ray, color='r', alpha=0.3,
label=f'Rayleigh 累计')
ax6.fill_between(S_int/sigma, 0, dmg_dir, color='k', alpha=0.3,
label='Dirlik 累计')
ax6.plot(S_int/sigma, dmg_ray, 'r-', lw=1.5)
ax6.plot(S_int/sigma, dmg_dir, 'k-', lw=1.5)
ax6.set_title('⑥ 损伤密度 dD/dS', fontsize=12, fontweight='bold')
ax6.set_xlabel('S/σ'); ax6.set_ylabel('dD/dS')
ax6.legend(fontsize=9); ax6.grid(alpha=0.3)
ax6.set_xlim(0, 5)
txt = f'D(Rayleigh) = {D_ray:.3e}\nD(Dirlik) = {D_dir:.3e}\n误差 = {(D_ray/D_dir-1)*100:+.1f}%'
ax6.text(0.97, 0.95, txt, transform=ax6.transAxes, fontsize=9,
family='monospace', va='top', ha='right',
bbox=dict(boxstyle='round,pad=0.4', facecolor='#FFF8DC',
edgecolor='gray', alpha=0.9))
plt.suptitle('穿越分析完整流程可视化', fontsize=16, fontweight='bold', y=0.995)
plt.tight_layout(rect=[0, 0, 1, 0.98])
plt.show()
# ============================================================
# 控制台汇总
# ============================================================
print('='*60)
print(' 穿越分析全流程结果')
print('='*60)
print(f'① 信号时长 T = {T} s, 采样 fs = {fs} Hz')
print(f'② PSD 峰: {f[np.argmax(Pxx)]:.2f} Hz')
print(f'③ 谱矩: m0={m0:.4e}, m2={m2:.4e}, m4={m4:.4e}')
print(f'④ 穿越率: N0={N0:.3f}, Np={Np:.3f}, N(+1σ)={crossing_rate(sigma,m0,m2):.3f} /s')
print(f'⑤ 分布: γ={gamma:.4f}, D1={D1:.3f}, D2={D2:.3f}, D3={D3:.3f}')
print(f'⑥ 损伤: Rayleigh={D_ray:.4e}, Dirlik={D_dir:.4e}')
print('='*60)
浙公网安备 33010602011771号