矩阵传递法计算层状地基瑞利波弥散曲线与位移衰减

简介：针对层状地基中瑞利波传播与衰减特性研究，这份资源将矩阵传递法的理论推导与可运行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, dtype=complex) for idx, layer in enumerate(layers[:-1]): # 这里用占位矩阵示意，工程上要替换为完整的 Thomson-Haskell 层矩阵 T = np.eye(4, dtype=complex) * 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(f"Warning: 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³)	说明
1	5	120	800	1800	填土/软黏土
2	10	180	1000	1900	淤泥质黏土
3	15	250	1200	2000	粉质黏土
4	∞	350	1500	2100	粉砂/半空间

将半空间厚度写成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_type='normal'): """绘制不同频率下归一化位移随深度的变化""" 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, label=f'{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], method='brentq') 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 米。这说明工程中进行环境振动评估时，不能只用一个“影响深度”去设计隔振沟或者桩基方案，而要先看环境振动的优势频段。如果振动能量集中在 1～5 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"

建议在每个频率点都记录条件数，并和相邻频率点的结果做平滑性检查。如果某一段频率的条件数阶跃式上升，优先怀疑是层厚与波数的乘积过大，此时把该层的厚度减半再试。这种处理往往比盲目提高搜索迭代次数更有效。

本文还有配套的精品资源，点击获取