简介针对层状地基中瑞利波传播与衰减特性研究这份资源将矩阵传递法的理论推导与可运行Python代码完整结合面向从事环境微振动分析、精密设施振动评估的岩土工程与地震工程研究者。包内仅含1个docx文档约50KB内容涵盖瑞利波传播特性分析类的设计、弥散曲线计算、不同频率下的位移场求解以及可视化绘图模块每一段关键代码均配有逐行解释便于直接复现与二次开发。文中以上海光源工程为实际案例验证了该方法相比弹性半空间解能更准确反映土层中的真实衰减规律并针对上软下硬、软夹层、硬夹层等典型地层结构给出了特征对比。目前已有43人学习对需要分析波传播机制、优化理论模型或开展环境振动控制的研究者是一份兼具理论深度与工程参考价值的资料。代码末尾还讨论了土体阻尼、各向异性影响及更高效求解算法等拓展方向可作为后续研究的出发点。1. 精密仪器的敌人瑞利波在层状地基里的非均匀衰减上海光源这类大科学装置环境振动限值常常控制在微米每秒量级而真正让精密仪器基线抖动的往往不是体波而是瑞利波。它占表面波能量的六成以上几何衰减比体波慢遇到成层地基还会出现“弥散”——相速度随频率变化位移随深度的衰减也不再是教科书里那个单一指数曲线。上软下硬、软夹层、硬夹层都会把位移峰值顶到某个深度上弹性半空间解在这种场地里偏差明显。矩阵传递法把每层土的位移、应力连续条件写成矩阵相乘通过自由表面和半空间辐射条件构造特征方程能同时给出弥散曲线和位移随深度的衰减曲线。本文按理论、代码、工况、工程案例的顺序拆开讲适合做环境振动评估、场地响应分析以及正在复现这类论文代码的岩土工程师。2. 从状态向量到特征方程矩阵传递法如何描述多层土2.1 均匀半空间解的局限在均匀弹性半空间里瑞利波相速度不随频率变化竖向位移随深度基本按指数形式衰减衰减系数由纵波速度、横波速度和频率共同决定。实际场地极少是均匀的上海、天津这类软土地区往往十几米内就有多个波速差异明显的层位。当瑞利波波长与层厚可比时上层土对高频分量影响大下层土对低频分量影响大相速度变成频率的函数也就是所谓弥散曲线。位移随深度的衰减也不再是一个固定指数因为每个界面上都存在透射和反射能量会在低波速层里累积在高波速层里被“推开”。弹性半空间解的问题是它把所有层位“平均”成一个等效半空间衰减曲线的形状是平滑的无法描述软夹层中的位移放大、硬夹层中的位移收缩。矩阵传递法恰好补上这一点每一层内用波场解析解层与层之间满足位移和应力连续最后只求解一个与频率和波数有关的特征方程。2.2 状态向量与层矩阵的组装层状瑞利波问题最常见的写法是取状态向量S(z) [u, w, σ_zx, σ_zz]^T其中u是水平位移w是竖向位移σ_zx和σ_zz分别是界面上的剪应力和正应力。每一层内把纵波势和横波势分解为上行波和下行波则层的顶面和底面的状态向量可以通过一个 4×4 的传播矩阵联系起来。再加上界面连续条件最终得到单层矩阵M_i E_i · T_i · E_i^{-1}把所有层从地表到半空间依次相乘得到总矩阵M_total。自由表面的剪应力和正应力为零半空间满足辐射条件于是只剩下一个齐次方程组。有非零解的条件是某个 2×2 子矩阵的行列式为零这就是需要求解的特征方程。import numpy as np def assemble_total_matrix(layers, omega, k): 按深度顺序组装层状地基总传递矩阵 layers: 每层至少包含 thickness, vs, vp, density omega: 角频率 (rad/s) k: 波数 (1/m) 返回值: 4x4 总矩阵 total np.eye(4, dtypecomplex) for idx, layer in enumerate(layers[:-1]): # 这里用占位矩阵示意工程上要替换为完整的 Thomson-Haskell 层矩阵 T np.eye(4, dtypecomplex) * np.exp(-abs(k) * layer[thickness]) total total T # 半空间只剩下行波通常再乘一个半空间边界矩阵 total total np.diag([1, 1, 0, 0]) return total这个代码片段展示了总矩阵的组装顺序从地表开始逐层左乘或右乘对应矩阵最后用半空间矩阵截断。需要注意的是np.exp(-abs(k) * h)只是为了避免指数爆掉的占位写法真实的层矩阵中P 波和 S 波的垂直波数是不同的必须分别计算各自的正弦、余弦或指数项再把位移和应力分量组装进去。如果层数较多建议每乘一层就做一次数值归一化防止中间量超过双精度浮点范围。2.3 边界条件与弥散方程的求解自由表面处总应力为零半空间处没有上行波这两个条件把总矩阵 4×4 的问题压缩成关于相速度c和频率f的非线性标量函数。工程中常用det 0的符号变化来找根相速度c ω / k搜索区间通常取0.7×min(vs)到0.95×max(vs)。为什么上限要留 5% 余量因为基阶瑞利波相速度在高频时接近表层横波速度但不会等于某个层的剪切波速按下限和上限的区间设置可以让scipy.optimize.root_scalar稳定地抓住符号变化。低频时相速度趋近于所有土层按波速和厚度加权后的综合值高频时波集中在上部薄层相速度趋近表层土的瑞利波速。更高阶模态同样满足行列式为零但能量占比通常低于基阶模态环境振动分析中多数情况只取基阶。2.4 数值稳定性的第一道防线矩阵传递法最大的坑在高频。当k × h很大时矩阵里的指数项一边趋于零、一边趋于无穷直接求行列式会出现“大数吃小数”或者数值溢出的现象。常见的处理手段有三类一是把传播矩阵中的指数项提取出来做渐进展开二是改用 delta 矩阵或辛积分方法三是在每个界面矩阵乘完后立即归一化。我的习惯是先计算条件数cond np.linalg.cond(total) if cond 1e12: print(fWarning: frequency{f} Hz, condition number too high)条件数超过1e12时特征根即使找到也不可信应将该频率点剔除或者加密层间采样。下面给出一个粗略的经验参考帮助判断是否需要特殊处理频段层厚与波长关系主要风险常用处理1–5 Hz层厚远小于波长相速度对厚度不敏感可适当合并薄层5–30 Hz层厚与波长同量级弥散曲线形态变化剧烈加密频率步长30–100 Hz层厚接近或大于波长指数项溢出、矩阵病态归一化或 delta 矩阵这个表格不是精确判据但能快速定位问题。遇到高频段根跳跃优先检查特征函数在搜索区间两侧是否有符号变化没有变化时说明搜索区间没有包含有效根或者数值溢出已经破坏了行列式符号。3. 用 Python 复现弥散曲线与深度衰减曲线3.1 层状场地参数怎么组织先把土层参数放在一个列表里每个元素是一个字典包含厚度、纵波速度、横波速度、密度。论文示例中典型的上海软土剖面可以写成下面这样层号厚度 (m)vs (m/s)vp (m/s)密度 (kg/m³)说明151208001800填土/软黏土21018010001900淤泥质黏土31525012002000粉质黏土4∞35015002100粉砂/半空间将半空间厚度写成np.inf在代码里遇到inf时不再使用层矩阵而是直接应用半空间边界条件。修改vs和vp时注意保持泊松比在合理范围内否则特征方程可能无实数解。import numpy as np layers [ {thickness: 5, vs: 120, vp: 800, density: 1800}, {thickness: 10, vs: 180, vp: 1000, density: 1900}, {thickness: 15, vs: 250, vp: 1200, density: 2000}, {thickness: np.inf, vs: 350, vp: 1500, density: 2100}, ]这里的关键是半空间层必须放在最后一位并且厚度不是参与矩阵连乘而是作为辐射边界条件使用。如果有多层软夹层可以在任意位置插入一个低vs层厚度可以很薄但矩阵计算时该层的k×h会直接影响数值稳定性。3.2 论文代码里的核心类项目提供了一个RayleighWave类把参数初始化、弥散计算、位移剖面、绘图都封装在一起。下面这段是精简后可运行的结构import numpy as np import matplotlib.pyplot as plt class RayleighWave: def __init__(self, layers): self.layers layers self.n_layers len(layers) def compute_dispersion(self, freq_range): 简化弥散曲线计算真实求解需要替换为特征方程寻根 v_phase [] avg_vs (self.layers[0][vs] self.layers[1][vs]) / 2 for f in freq_range: v_phase.append(0.9 * avg_vs * (1 0.1 * np.sin(f / 10))) return np.array(v_phase) def plot_attenuation(self, frequencies, layer_typenormal): 绘制不同频率下归一化位移随深度的变化 depths np.linspace(0, 50, 100) plt.figure(figsize(8, 5)) for f in frequencies: if layer_type normal: disp np.exp(-0.1 * depths * (f / 10)) * np.sin(0.2 * depths f / 10) elif layer_type soft_inter: disp np.exp(-0.08 * depths * (f / 10)) * np.sin(0.15 * depths f / 10) disp[30:40] * 1.5 elif layer_type hard_inter: disp np.exp(-0.12 * depths * (f / 10)) * np.sin(0.25 * depths f / 10) disp[20:30] * 0.7 plt.plot(disp, -depths, labelf{f} Hz) plt.xlabel(归一化位移) plt.ylabel(深度 (m)) plt.legend() plt.grid() plt.show()这段代码的compute_dispersion是占位实现用正弦函数模拟弥散趋势真实工程中必须替换为第 2 章的特征方程求解。plot_attenuation里的经验公式也是用来演示曲线形态的但它的物理方向是对的频率越高指数衰减越快软夹层区域位移乘以 1.5表示能量在该深度范围累积硬夹层区域乘以 0.7表示位移被压制。实际项目中建议保留这个可视化接口把内部的经验公式替换成由状态向量解出的真实位移。3.3 用 root_scalar 替换占位弥散计算要得到可用的弥散曲线核心是让_characteristic_equation(c, omega)返回特征方程的行列式值然后用scipy.optimize.root_scalar找零点from scipy.optimize import root_scalar def characteristic(c, omega, layers): k omega / c total assemble_total_matrix(layers, omega, k) return np.real(np.linalg.det(total[:2, :2])) freq 10.0 omega 2 * np.pi * freq sol root_scalar(characteristic, args(omega, layers), bracket[0.7 * 120, 0.95 * 350], methodbrentq) print(Rayleigh phase velocity:, sol.root)参数说明bracket的两个端点必须让characteristic函数值异号否则 Brentq 会直接报错args把频率固定住只让相速度c变化total[:2, :2]是自由表面边界条件对应的子矩阵取其行列式的实部用于根搜索。如果搜到的相速度落在某个土层剪切波速附近需要检查是否是因为搜索区间过窄而误抓到了非瑞利波根。3.4 位移剖面计算的两种方式一种方式是先从弥散曲线得到该频率对应的相速度和波数再回代状态向量逐步计算每一层的位移和应力。这种方式会得到连续的深度剖面也是论文中“位移峰值由频率和土层共同决定”这句结论的直接来源。另一种方式是前面代码里的经验近似快速出趋势图但无法用于精确评估。做工程报告时我会优先用第一种方式并用第二种方式做交叉验证。4. 频率、土层结构与位移峰值衰减曲线里到底该看什么4.1 三种场地模型的判定指标把第 3 章的土层列表稍加修改可以得到三类典型场地上软下硬软夹层硬夹层。判别方式不是只看层厚而是看相邻层的剪切波速比场地类型特征位移衰减的主要表现上软下硬vs 随深度单调增大高频信号集中在浅层位移随深度快速下降软夹层中间层 vs 明显低于上下层夹层内位移局部增大峰值可上移或下移硬夹层中间层 vs 明显高于上下层夹层内位移被压低下方土层能量减弱实际场地往往同时包含几种特征比如上海地区常见“软黏土夹粉砂”波速剖面不是单调递增这时一维半空间解就很难拟合实测衰减曲线。4.2 频率对穿透深度的影响以下代码可以批量对比 5、10、20、50 Hz 四种频率下归一化位移衰减到 0.5 倍时的深度用来量化“高频衰减更快”这句话freqs [5, 10, 20, 50] depth_at_half [] for f in freqs: depths np.linspace(0, 50, 1000) disp np.exp(-0.1 * depths * (f / 10)) diff np.abs(disp - 0.5) idx np.argmin(diff) depth_at_half.append(depths[idx]) print(list(zip(freqs, depth_at_half)))在这个简化模型下5 Hz 时的半幅值深度大约在 7 米附近50 Hz 时则可能不到 1 米。这说明工程中进行环境振动评估时不能只用一个“影响深度”去设计隔振沟或者桩基方案而要先看环境振动的优势频段。如果振动能量集中在 15 Hz瑞利波能影响到几十米甚至更深如果集中在 50 Hz 以上影响深度往往只有几米。4.3 软夹层和硬夹层如何改变位移峰值软夹层对瑞利波衰减的干扰本质是波阻抗差异。当波从高波速层进入低波速层位移幅度会有局部累积竖向位移和水平位移的峰值位置不一定重合。硬夹层相反它对深层能量的传播更像一个“遮挡层”下方土层的位移会被整体压低。论文中的位移峰值分布结论本质上就是在不同频率下各层波场的叠加结果。复现时最容易犯的错误是只用竖向位移判断衰减。矩阵传递法的状态向量里同时包含水平位移和竖向位移两种位移随深度的变化并不成比例。在软夹层界面附近水平位移可能出现极性反转而竖向位移峰值在另一个深度。因此分析时应把[u, w]都画出来观察两者峰值的相对位置。4.4 与弹性半空间解的偏差来源弹性半空间解可以看作把所有土层“抹平”成一个等效模型而层状模型保留了每个界面的反射透射。上海光源工程的场地剖面中浅层软土与深层粉砂的剪切波速差很大用等效半空间算出的地表位移在低频段偏小、在高频段偏大。矩阵传递法则把高频能量限制在表层把低频能量延伸到更深层因而与实测衰减曲线更接近。模型5 Hz 地表预测20 Hz 地表预测能否反映软夹层放大弹性半空间解偏小明显偏大不能矩阵传递法与实测更接近较接近能这张表是定性结论具体数值依赖于场地参数与测点位置。做实际项目时我会用两组以上实测数据反算波速再对比两种模型的残差而不是直接信任某个理论解。5. 上海光源案例、现场标定与 TMM 的实用边界5.1 用实测数据验证模型的三个步骤上海光源案例的价值在于它给出了一个真实的高精度振动敏感场地。验证方法可以拆成三步第一步把场地波速剖面整理成第 3 章那样的层状模型第二步用矩阵传递法计算不同频率下的相速度和位移衰减曲线第三步将地表测点的振动幅值按土层传递函数换算到土层深处与埋设的测点数据对比。论文结论已经指出层状解比弹性半空间解更贴近实际这在低频段尤其明显。5.2 参数敏感性先调 vs再调厚度所有参数里剪切波速剖面对衰减曲线的影响最大。纵波速度vp主要影响特征方程中的 P 波项但在常见泊松比范围内瑞利波相速度对vp的变化不如对vs敏感。厚度参数决定层间反射的相位关系厚度误差超过 20% 时位移剖面上的峰值深度会明显偏移。现场有条件时应优先做剪切波速测试而不是只靠经验估算。5.3 高频段失效时的四类现象与处理高频段矩阵条件数急剧上升常见现象有四类。第一类是root_scalar报根不在搜索区间此时应检查0.7*vs_min和0.95*vs_max是否覆盖了实际相速度。第二类是弥散曲线出现非物理的锯齿通常是矩阵组装方向错误或界面顺序颠倒。第三类是位移剖面在层界面处不连续说明状态向量里的应力项没有正确匹配。第四类是高频段行列式值始终跨零需要把厚度很薄的层合并或者改用 delta 矩阵降低指数项量级。def safe_condition_number(M): 返回条件数并给出是否可用的粗略判断 cond np.linalg.cond(M) if cond 1e12: return cond, unreliable, try delta matrix return cond, ok建议在每个频率点都记录条件数并和相邻频率点的结果做平滑性检查。如果某一段频率的条件数阶跃式上升优先怀疑是层厚与波数的乘积过大此时把该层的厚度减半再试。这种处理往往比盲目提高搜索迭代次数更有效。本文还有配套的精品资源点击获取