1. 一张 6×6 矩阵摆在面前，判读为什么必须靠图
第一次被问到“这批样品到底算不算强各向异性”的时候，我手里只有一张 6×6 的弹性刚度矩阵。数据完整、精度也够，但我盯着那三十多个数字看了半小时，能说出来的只有“C11 比 C33 大”——这句话对判断方向依赖性几乎没有帮助。弹性各向异性的麻烦就在这里：它描述的是波速、模量随方向变化的全景，而刚度矩阵是这套全景的压缩编码，人眼无法直接解码。我后来把整套流程写成了一个工具，也就是这篇要讲的VELAS程序，名字取的是“弹性各向异性可视化与分析”这几个词的首字母组合，定位很朴素：吃进一个刚度张量，吐出三张速度面、一组方向极值、一张各向异性强度报告。
适合读这篇的人，大概是三类。第一类是做岩石物理或岩心弹性表征的，手里有一批超声测量反演出来的 Cij，需要快速判断这块岩心是不是该按横向各向同性处理。第二类是材料方向，关注晶体、织构化金属、增材制造件里的方向依赖性，想直观看到“最软的加载方向在哪”。第三类是做地震或超声数值模拟的，需要从张量里抽出相速度、群速度、剪切波分裂这些输入参数。三类人的共同点是：不满足于知道“有各向异性”，而是需要知道“在哪个方向、强到什么程度、能不能用一个参数横向排序”。这也是我把 VELAS 做成一个带交互界面的程序，而不是一段一次性脚本的原因。
1.1 各向同性假设在哪一步开始崩掉
用各向同性近似处理实际问题是常见操作，本身没问题，问题在于近似什么时候失效。我做过一轮对比：拿同一批样品，一边用完整刚度张量算方向速度，一边用两个拉梅常数加密度算各向同性速度，然后看两者的偏差随方向怎么变。结论是——偏差不是均匀分布的小扰动，而是集中在少数几个方向上，且最大偏差往往正好落在工程上最关心的方向。比如很多沉积类样品的最软方向接近垂直于层理，而受力分析和取心设计恰恰关心这个方向。
这就带来一个实操上的判断：如果你的应用只涉及某一个固定方向，各向同性近似可能完全够用；一旦涉及方向扫描（比如井壁稳定性分析、超声波探头的多角度扫查、构件多轴加载评估），近似就会在特定角度上产生明显误差，而这种误差用全局平均值是看不出来的。VELAS 的第一个价值点就在这儿——它把“误差在哪个方向最大、多大”变成一张可以旋转查看的曲面，而不是一行统计数字。我一般会先跑一遍max|v_aniso - v_iso| / v_iso这个相对偏差场，如果峰值出现在 5% 以内，各向同性处理基本没风险；超过 15% 就得老老实实按各向异性走。
1.2 程序要回答的三个具体问题
把需求收敛之后，我发现所有使用场景都能归到三个问题上，VELAS 的功能也就按这三个问题来切分：
| 问题 | 输出物 | 典型用途 |
|---|---|---|
| 波速/模量随方向怎么变 | 三分支速度面（准P、准S1、准S2）、方向切片曲线 | 判断最软/最硬方向，确定测量角度 |
| 极值有多大、在哪个方向 | 全局极值表（方向余弦 + 数值 + 相对偏差） | 给数值模拟提供边界参数 |
| 不同样品之间怎么横向比较 | 各向异性强度指数、Thomsen 参数、分裂强度表 | 批量筛选、质量分级 |
第一类输出是给人看的，第二类是给下游程序用的，第三类是给报告用的。三者必须来自同一套计算内核，否则你会遇到“图上看着差不多、表里数值差一截”的尴尬。VELAS 的做法是：所有输出统一从同一个相速度求解器出发，极值搜索和图表渲染都只是这个求解器的不同消费者。这个设计看起来是常识，但我见过不少自研脚本，图上画的和表里报的是两套代码，最后对不上账，排查起来极其痛苦。
1.3 最小输入契约与真实数据的卡点
一篇讲工具的文章如果不讲输入格式，基本等于没讲。VELAS 要求的最小输入是：一个 6×6 的刚度矩阵，加一个标量密度。矩阵按标准 Voigt 顺序排列，即行列为 11、22、33、23、13、12 对应的下标，单位要求随密度一起给，程序内部统一换算。除此之外全部可选：样品名、测量温度、数据来源备注。
真实数据的卡点几乎全在“不完整”上。最常见的三种：只给了 21 个独立分量的一维列表，顺序还得靠猜；给了柔度矩阵而不是刚度矩阵，单位还是 1/TPa；给了 9 个分量并说明“近似取为各向同性”，但多出来的那些其实不为零。我在 VELAS 里加了一个输入检查层，专门处理这几种情况——列表按标准顺序展开并做对称性校验，柔度矩阵自动求逆并提示条件数，声明为各向同性但残留偏差超过阈值的直接给出警告而不是静默通过。这些检查写起来不到两百行，但省下来的排查时间远超预期，后面第 4 章会详细讲怎么组织这一层。
2. Christoffel 方程：把 21 个独立常数压成一个 3×3 特征值问题
弹性各向异性的数学核心其实很干净：给定传播方向 n，构造一个 3×3 的实对称矩阵，求它的特征值，开个方就得到该方向的三个相速度。这个矩阵叫 Christoffel 矩阵，表达式是Γ_ik = C_ijkl · n_j · n_l / ρ。理论上就这么一行。但真要从四阶张量落到代码里，中间有几个地方特别容易出错，而且错了之后现象很像“物理效应”，很难第一时间怀疑到代码上。
2.1 Voigt 下标映射错一位，后面全盘皆错
四阶张量到 6×6 矩阵的压缩依赖 Voigt 约定，而这个约定在不同文献里的下标顺序并不完全统一。VELAS 采用最常见的顺序：1→11、2→22、3→33、4→23、5→13、6→12。注意 4、5、6 对应的三个剪切分量，有些教材排成 23、31、12，有些排成 12、13、23，一旦你的输入数据按另一套顺序排，程序读进来的矩阵在数值上完全合法、物理上完全错位。
我最初就栽在这儿。那批样品是别人给的表格，我按自己的约定直接读，出来 P 波速度面看着挺正常，但 S 波分裂的方位角怎么都和实测对不上。排查了两天才想起来去核对下标约定，换成对方的顺序后一切吻合。从那以后 VELAS 的输入层强制带一个voigt_order字段，默认值写在配置里，不填就按默认走，填了就必须校验。这个教训值得单独强调：刚度矩阵的错位不会报错，只会给你一个看起来很像真的、但错了的结果。
2.2 六个分量的装配表达式与自检方法
Christoffel 矩阵是实对称的，六个独立分量各有一个显式表达式。下面这份是我在 VELAS 里实际用的版本，方向余弦记作n1, n2, n3，刚度矩阵用 0 起始的下标c[0][0]对应C11：
import numpy as np def christoffel(C, n): """C: 6x6 刚度矩阵(Pa), n: 单位方向向量""" c = C n1, n2, n3 = n G = np.zeros((3, 3)) G[0, 0] = (c[0,0]*n1*n1 + c[5,5]*n2*n2 + c[4,4]*n3*n3 + 2*(c[0,5]*n1*n2 + c[0,4]*n1*n3 + c[5,4]*n2*n3)) G[1, 1] = (c[5,5]*n1*n1 + c[1,1]*n2*n2 + c[3,3]*n3*n3 + 2*(c[1,5]*n1*n2 + c[3,5]*n1*n3 + c[1,3]*n2*n3)) G[2, 2] = (c[4,4]*n1*n1 + c[3,3]*n2*n2 + c[2,2]*n3*n3 + 2*(c[2,4]*n1*n2 + c[3,4]*n1*n3 + c[2,3]*n2*n3)) G[0, 1] = (c[0,5]*n1*n1 + c[1,5]*n2*n2 + c[3,4]*n3*n3 + (c[0,1] + c[5,5])*n1*n2 + (c[0,3] + c[5,4])*n1*n3 + (c[3,5] + c[1,4])*n2*n3) G[0, 2] = (c[0,4]*n1*n1 + c[3,5]*n2*n2 + c[2,4]*n3*n3 + (c[0,3] + c[5,4])*n1*n2 + (c[0,2] + c[4,4])*n1*n3 + (c[2,5] + c[3,4])*n2*n3) G[1, 2] = (c[5,4]*n1*n1 + c[1,3]*n2*n2 + c[2,3]*n3*n3 + (c[1,4] + c[3,5])*n1*n2 + (c[2,5] + c[3,4])*n1*n3 + (c[1,2] + c[3,3])*n2*n3) G[1, 0], G[2, 0], G[2, 1] = G[0, 1], G[0, 2], G[1, 2] return G这份表达式建议不要盲信，一定要做自检。最有效的自检方法是退化到各向同性：代入C11=C22=C33=λ+2μ、C12=C13=C23=λ、C44=C55=C66=μ，把所有剪切-正应变耦合项设为零，此时 Christoffel 矩阵应该严格等于μ·I + (λ+μ)·n⊗n，三个特征值分别是λ+2μ和两个重根μ。如果这条自检过不了，后面的所有图都不用看了。VELAS 的单元测试里就放着这条用例，还有一条立方对称的用例（沿 [100] 和 [111] 方向的解析解是已知的），双保险。
2.3 特征值排序带来的分支交叉问题
速度求解本身一行就够：v = sqrt(eigvalsh(Γ)/ρ)。坑在“三个分支怎么归属”。eigvalsh返回升序排列的特征值，对应下来就是“最慢、中间、最快”。按大小排序看起来天经地义，但在某些方向附近，两个剪切波速度会非常接近，一旦穿过某个方向发生顺序交换，你画出来的 S1 曲面和 S2 曲面就会在那个位置强行“拧”一下，形成一道折痕。
实测中这道折痕非常显眼，而且会被误认为是物理上的剪切波简并。VELAS 的处理办法是不按大小排序，而是按与前一采样点的匹配度排序。具体做法：球面采样时按邻接顺序遍历（比如按纬度从一极扫到另一极），对每个点算完三个特征向量后，与上一个采样点的三个特征向量做内积，选内积绝对值最大的一一配对；准P 波那一支的特征向量总是与传播方向最接近，可以先锁定，剩下两支再做匹配。这样得到的曲面是连续的，简并点附近会出现两支贴合的“接缝”，那才是真实的物理现象。
代价是遍历顺序必须严格，不能并行乱序。VELAS 里这一步是单线程串行的，牺牲一点速度换来曲面的连续性，我认为完全值得。
2.4 球面采样的三种方案与代价
速度面要在单位球上采样，采样方案直接决定图的质量和极值搜索的可靠性。我用过三种，各自的脾气不太一样：
- 经纬网格：最直观，
theta从 0 到 π、phi从 0 到 2π 均匀取点。优点是可以直接喂给plot_surface，出图快；缺点是两个极点的采样密度爆炸，而赤道附近偏疏，做面积加权统计时会严重高估极区。做展示图可以，做定量分析不行。 - Fibonacci 球：按黄金角在球面上撒点，近似等面积分布，点数可控。做极值粗筛和面积平均最合适，两万点就能覆盖得相当均匀。缺点是不能直接构成四边形网格，画曲面得走三角化。
- 细分二十面体：从正二十面体开始递归细分，三角形大小均匀，天然适合
plot_trisurf。缺点是细分层级只能按 4 倍增长，点数不够灵活。
VELAS 的实际策略是混合：定量计算（极值、平均值、面积加权）用 Fibonacci 球，展示曲面用经纬网格，两套采样共用同一个求解器。这样图上看到的东西和表里报的数字是一致的，只是采样点不同，误差在可控范围内。如果只允许选一种，我建议选 Fibonacci，因为定量结果的可信度比图的顺眼程度重要得多。
3. 速度面画出来只是开始：极值、指数与分裂的判读
三张曲面渲染出来，第一次看确实很惊艳，但很快就会发现两个问题。第一，旋转着看能感觉到“有方向性”，但说不出具体在哪个方向最强。第二，两张不同样品的曲面摆在一起，肉眼根本分不出谁的各向异性更强。可视化解决的是直觉问题，不解决量化问题，所以 VELAS 在渲染之外还做了三件事：极值搜索、强度指数、分裂分析。
3.1 从粗筛到球面局部优化
网格上取最大值这种做法，我一开始也用过，误差大到不能用。原因是速度极值出现的方向往往不在采样点上，网格越粗偏差越大；而且对于剪切波分支，极值附近曲面可能很平，网格最大值和真实极值能差出百分之几。
VELAS 的做法是两段式。第一段用 Fibonacci 球撒两万个点做粗筛，对每个分支维护一个“候选集”，保留数值排名前 5% 的点。第二段对候选点做球面局部优化：把方向参数化为球坐标(theta, phi)，用 Nelder-Mead 单纯形法最小化负速度。选 Nelder-Mead 而不是梯度法，是因为在剪切波简并点附近曲面不可导，梯度法会在那儿原地打转。
局部优化完之后还有一步对称性去重。一个立方对称的晶体，沿 [100]、[010]、[001] 三个方向的 P 波速度完全相等，优化器会给你三个几乎相同的最优解，直接列进表里会让人以为有三个不同的极值方向。VELAS 的去重逻辑是按晶体学等价方向分组，每组只保留一个代表，同时把该组的方向数列在备注里。这一步看起来是锦上添花，实际上直接决定了极值表的可读性。
提示：粗筛点数不要低于一万。我试过五千点，某些近退化样品的极值会被漏掉，因为极值所在的那一小片区域恰好落在采样点的空隙里。
3.2 Thomsen 参数只在 VTI 语境里成立
提到各向异性强度，很多人的第一反应是 Thomsen 三参数ε、γ、δ。它们确实好用，定义也简单：
ε = (C11 - C33) / (2·C33) γ = (C66 - C44) / (2·C44) δ = ((C13 + C44)² - (C33 - C44)²) / (2·C33·(C33 - C44))但这套参数的适用前提是横向各向同性且对称轴沿 3 轴（也就是 VTI）。VELAS 在报告里会把这三个值算出来，但前面加了一行状态标记：如果输入张量的对称性检测结果不是 VTI，这一行会标成“参考值，语义不成立”。这么做是因为我见过太多把 HTI 样品硬套 Thomsen 参数、然后困惑于“为什么 γ 是负的”的情况。HTI 的正确做法是先旋转到 VTI 参考系，或者直接用下一节讲的通用指数。
δ还有一个更隐蔽的坑：当C33接近C44时，分母趋于零，δ会爆掉。这种情况在近各向同性样品上真的会出现，程序必须显式拦截，不能让一个1e6量级的数字混进报告里。
3.3 一个能跨样品横向比较的通用指数
要横向排序，需要一个不依赖对称性假设的标量。VELAS 用的是 Voigt-Reuss 上下界构造的各向异性指数：分别算出 Voigt 平均和 Reuss 平均的体积模量与剪切模量，然后组合成一个偏离度。
| 量 | 表达式 |
|---|---|
| Voigt 体积模量 | K_V = [(C11+C22+C33) + 2(C12+C13+C23)] / 9 |
| Voigt 剪切模量 | G_V = [(C11+C22+C33) - (C12+C13+C23) + 3(C44+C55+C66)] / 15 |
| Reuss 体积模量 | 1/K_R = (S11+S22+S33) + 2(S12+S13+S23) |
| Reuss 剪切模量 | 15/G_R = 4(S11+S22+S33) - 4(S12+S13+S23) + 3(S44+S55+S66) |
| 各向异性指数 | A = 5·(G_V/G_R) + (K_V/K_R) - 6 |
其中S是刚度矩阵的逆（柔度矩阵）。这个指数的好处是各向同性时严格为零，且量纲无关，跨样品可比。我把一批样品按这个值排序过，结果和超声多角度扫描的直观判断基本一致，只有一两个边界情况需要人工复核。需要提醒的是，这个指数只反映“偏离各向同性的程度”，不反映“哪个方向偏离最大”，两者要配合着看。
3.4 三分支曲面的读图规则
速度面的读图有几个固定套路，我在 VELAS 的帮助文档里直接写进去了，因为新用户最容易在下面这几点上绕圈：
第一，准P 波曲面永远是最外层的那个，因为它速度最高；两个准S 波曲面在内层，且在某些方向上会贴合，形成剪切波简并。看到内层曲面“粘”在一起，不要怀疑渲染出了问题，那是真实的物理。
第二，曲面上的“尖点”和“凹坑”位置比曲面的整体形状更重要。各向异性带来的信息几乎全在极值附近的局部形态里，整体轮廓反而看不出太多。VELAS 允许在曲面上叠加极值标记点，方便定位。
第三，半径的归一化方式会改变视觉印象，这一点后面第 5 章会专门讲，是踩过的坑。
第四，相速度面不等于能量传播速度面。这一条最容易被忽略。相速度是波前法向的速度，能量实际传播的方向（群速度）是相速度面在波矢空间的梯度方向，两者只有在高对称方向上才重合。VELAS 提供群速度的数值计算（对相速度面做切向差分），但在图上默认不画，因为把相速度和群速度两条曲线画在同一张图里会让初学者彻底混乱。如果你要做射线追踪或者探头角度设计，务必用群速度那一套输出。
4. 从脚本到程序：渲染、容错与批处理上的取舍
自己用和给别人用之间隔着一道鸿沟。脚本阶段我只需要“跑出图来”，程序阶段要考虑的是：别人拿到一个格式不太标准的矩阵会不会直接崩、图在别人电脑上能不能正常转、一次处理五十个样品要花多久。VELAS 在这三件事上做的取舍，比算法部分更值得写下来。
4.1 渲染后端的三档选择
三维曲面渲染有几条路，我用过三档，适用场景差别很大：
| 方案 | 优点 | 缺点 | 适合谁 |
|---|---|---|---|
Matplotlib 的plot_surface | 零额外依赖，出静态图快，能直接导出 PDF | 交互旋转卡顿，光照简单，遮挡处理粗糙 | 只出报告插图 |
| 三角化 + 散点渲染 | 采样更均匀，能画 Fibonacci 点云 | 视觉效果一般，不好做光照 | 快速检查 |
| 专业三维库（带交互窗口） | 旋转流畅，支持光照、剖切、拾取 | 安装体积大，部分环境缺图形驱动 | 探索式分析 |
我最终的方案是双轨：批处理出报告走 Matplotlib，探索式分析开交互窗口。两者共用同一份速度计算代码，只是渲染层不同。这个分层很关键，因为一旦把计算逻辑写进渲染回调里，换后端就意味着重写。
Matplotlib 出曲面时有两个细节值得说。一是采样网格别太密，plot_surface的默认网格数是 50×50，对速度面这种光滑曲面，加密到 200×200 视觉提升有限，但渲染时间翻好几倍；我一般用 120×240，配合rstride=2, cstride=2，能兼顾平滑和速度。二是极区的退化三角形会让光照出现暗斑，把纬度采样从[0, π]稍微内缩一点（比如[1e-3, π-1e-3]）就能规避。
4.2 输入检查该放在哪一层
我最初把输入检查散落在各处：读文件时检查一次，算速度时再检查一次，画图前又检查一次。结果是错误信息五花八门，用户根本不知道自己错在哪。后来把它整个收到一个独立的校验层，所有入口都必须过这一层，规则也统一成几条：
- 对称性：
C[i][j]与C[j][i]的相对偏差超过1e-6直接报错，不静默平均。静默平均会让用户误以为数据本来就是对称的，掩盖了上游的录入错误。 - 正定性：所有特征值为正。刚度和柔度都应该是正定的，出现负特征值基本意味着单位搞错或者分量错位。
- 单位：必须显式给出，不允许“猜”。VELAS 内部统一换算到 Pa 和 kg/m³，配置里写
GPa和g/cm³的都行，但必须写。 - 对称性等级判定：自动比对各个分量，判断最接近哪种对称类（各向同性、立方、六方、正交、单斜、三斜），并把“与理想对称的最大偏差”报出来。这一步的实用价值极高——很多时候用户自己也不确定样品该按哪种对称处理，让程序给个参考再人工确认，比拍脑袋强。
第四条实现起来不复杂：把矩阵按对称类的约束条件逐项比较，算相对偏差，取最小值对应的类别。六方对称的判断要小心，因为它的约束是C11=C22、C13=C23、C44=C55、C66=(C11-C12)/2、其余为零，其中最后一条是恒等式约束，偏差大的话说明数据有问题。
4.3 批处理与报告的组织方式
单样品交互是一回事，一批样品跑一遍是另一回事。VELAS 的命令行模式接一个目录，自动扫描里面的数据文件，逐个计算、逐张出图、最后汇总成一张总表。总表按各向异性指数降序排，第一列是样品名，后面依次是三个方向的最大速度、最小速度、相对偏差、强度指数、对称类判定结果、警告信息。
这个流程里有两点经验。第一，单个样品失败不能中断整批。早期版本一个文件格式不对，整个批处理就崩了，五十个样品重跑一遍很痛苦。现在每个样品包在独立的错误捕获里，失败信息写进汇总表的备注列，其余照常。第二，图的命名要可排序。用样品名加分支名（比如S01_P.png、S01_S1.png），不要用时间戳，否则汇总的时候排序会乱。
报告我倾向输出 PDF 加 CSV 两份。PDF 给人工看，CSV 给下游脚本用。CSV 里只放数值和方向余弦，不放任何格式化字符串，避免下游解析时因为逗号或者中文标点出问题。
5. 实测中反复出现的四类问题及其排查路径
这一章是整篇里最不想写、但最该写的部分。下面四个问题，每一个都让我花过至少半天，而且每一个都有“看起来像物理、实际是程序”的迷惑性。
5.1 单位制：速度整体偏了三十倍
现象：某批样品的速度面形状完全正常，各向异性比例也对，但绝对数值比预期小了一个数量级。第一反应是密度给错了，核对之后密度没问题。
根因在单位混用。刚度矩阵用的是 GPa，密度用的是 kg/m³。v = sqrt(C/ρ)里，如果 C 按 GPa 的数值（比如 80）代入，ρ 按 kg/m³ 的数值（比如 2700）代入，得到sqrt(80/2700) ≈ 0.172，单位既不是 m/s 也不是 km/s，而正确值应该是sqrt(80e9/2700) ≈ 5443 m/s。两者差了sqrt(1e9) ≈ 31623倍……实际观测到的偏差没那么夸张，是因为中间还混了别的换算，但量级完全对不上。
这里有个极其好用的心算技巧：刚度用 GPa、密度用 g/cm³ 时，sqrt(C/ρ)的数值直接就是 km/s。因为1 GPa / (1 g/cm³) = 1e9 Pa / 1e3 kg/m³ = 1e6 (m/s)²，开方正好是 1000 m/s，也就是 1 km/s。这个换算我写在了 VELAS 的配置注释里，也写在了每一份输出的页脚，因为它是这类错误最高发的地方。
排查路径现在固定成三步：先看速度量级是否落在材料常识范围内（岩石 2–8 km/s，金属 3–7 km/s，陶瓷更高）；再看 P 波和 S 波的比值是否在 1.5–2.5 之间；最后才去看曲面形状。顺序反了就会在形状上浪费大量时间。
5.2 归一化制造出来的假各向异性
现象：两个样品，各向异性指数差了将近一倍，但把速度面并排放在一起看，形状几乎一样。
根因是渲染时的归一化方式。VELAS 早期版本的曲面半径直接取速度值，两个样品如果整体速度水平不同，画出来的曲面大小就不一样，视觉上一大一小，很容易被误读成“各向异性强度不同”。反过来，如果为了“放在一起好比较”而把每个曲面都归一化到单位球，各向异性的视觉差异反而被抹平了，两个形状差别很大的样品看起来也差不多。
最后的处理是提供三种显示模式，让用户显式选：
- 绝对模式：半径就是速度值，适合看单个样品的真实量级。
- 相对偏差模式：半径画
(v - v_mean)/v_mean，适合看方向依赖的形状，因为它把平均速度的影响剥掉了。 - 归一化模式：半径画
v / v_max，适合把多个样品叠在一张图里比形状。
默认用相对偏差模式，因为绝大多数时候用户关心的是“哪个方向软、软多少”，而不是“绝对速度是多少”。这个默认值改过之后，关于“为什么两个样品看起来一样”的疑问基本消失了。
5.3 近退化张量把参数算成 0/0
现象：某个样品的δ参数输出成了一个巨大的数，其他参数都正常。
排查下来是这样：这个样品的数据本身就接近各向同性，C33和C44的差值很小，而δ的分母里正好有(C33 - C44)这一项。分母趋近于零，分子却是有限值，结果自然爆掉。
这个问题不能靠调容差糊过去，必须显式检测。VELAS 的做法是：算δ之前先算分母的绝对值，如果小于C33的千分之一，就把δ标记为“不可靠”并输出NaN，同时在报告的警告列里写清楚原因。用户看到NaN会来问，看到1e7反而可能直接拿去用了，这是更危险的情况。
近退化还有另一个表现：两个剪切波速度几乎相同，第 2 章讲的分支匹配算法在这个位置会频繁切换配对，导致曲面上出现细小的锯齿。VELAS 在这种情况下会降低匹配的敏感度（用内积阈值而不是硬选最大），宁可让两支短暂重叠，也不制造假的分裂。
5.4 自遮挡与色标误导
现象：用户反馈“图上看到一块暗区，是不是速度特别低”。
实际上那是背面的曲面透过来被前面的曲面挡住了，加上默认色标是速度的连续映射，暗区就被误读成了低速区。三维曲面渲染在这方面的误导性很强，尤其是三张曲面叠在一起的时候。
VELAS 的应对有几条：默认只画准P 曲面，另外两张需要显式打开；色标改为按方向着色而不是按速度着色，速度信息用等值线标出来；提供一个“半透明外壳 + 内层实体”的组合模式，让内外层结构一眼可辨。另外强烈建议在报告里配一张方向切片曲线，因为二维曲线没有遮挡问题，读数也精确，是最可靠的判读依据。
提示：凡是从三维曲面图上读出“某个方向速度最低”这类结论，都建议回二维切片上复核一遍。三维图的视觉误差经常超过各向异性本身带来的差异。
6. 三类典型样品的判读记录
抽象讲完，落地还是得看具体数据长什么样。下面这三类是我处理得最多的，记录一下判读过程和结论，给遇到类似情况的同行一个参照。
6.1 横向各向同性的层状样品
这类样品最典型：沿层理方向和垂直层理方向的 P 波速度差异明显，剪切波在两个方向上的分裂行为也不同。VELAS 的对称性判定通常会给出六方对称，C11 ≈ C22、C13 ≈ C23、C44 ≈ C55，偏差都在百分之几以内。
判读的时候我先看准P 波曲面的形状：如果是标准的“橄榄形”，长轴沿对称轴，短轴在垂直于对称轴的平面内，那就是典型的 VTI。接着看 Thomsen 参数，ε一般在 0.1 到 0.3 之间，γ稍小一点，δ的变化范围最大、也最不稳定，这符合经验。最后看剪切波分裂强度，垂直方向分裂最强、沿对称轴方向分裂为零，这个规律非常稳定，可以用来交叉验证对称轴的方向有没有搞错。
有一点要提醒：这类样品的各向异性强度对不同物性参数的敏感度不一样。速度的差异可能只有百分之十几，但弹性模量的差异能到百分之几十，因为模量对速度是平方关系。做力学分析的时候要用模量那一套输出，别直接拿速度比值去代入。
6.2 立方对称的晶体类样品
立方对称的特点是三个主轴方向等价，C11=C22=C33、C12=C13=C23、C44=C55=C66，独立常数只有三个。这类样品的速度面很有意思：沿 [100] 和 [111] 两个方向的差异决定了各向异性的强弱，而且在某些材料里会出现“沿 [111] 反而更慢”的反常情况。
判读的关键是 Zener 比A = 2·C44/(C11 - C12)，它等于 1 时材料在剪切意义上各向同性。VELAS 报告里会把A和它偏离 1 的程度单独列出来。这个比值超过 1 和小于 1 对应的是两种完全不同的方向依赖模式，前者最硬方向在 [111]，后者在 [100]，画出来的曲面形态一眼就能分辨。
立方对称还有一个视觉特征：速度面有明显的四重对称瓣，沿三个主轴方向完全一样。如果画出来的图三个主轴方向不对称，那基本可以断定是数据录入出了问题，或者样品本身就不是立方对称、被误判了。
6.3 正交对称的复合材料样件
正交对称是这三个里最麻烦的，独立常数有九个，三个主轴方向互不等价。速度面既不是橄榄形也不是四重瓣，而是三轴椭球被“扭”过的样子，看第一眼很难建立直觉。
我处理这类数据的方式是先降维再看。VELAS 提供沿三个主平面的方向切片，也就是固定一个方向余弦为零、在另外两个方向构成的平面内扫角度。三张二维曲线拿到手，各向异性的主轴方向和强弱立刻清楚了，然后再回头看三维曲面，就不会迷路。
这类样品最容易出的问题是“误判成六方”。因为如果两个方向恰好接近，C11和C22差得不多，自动对称性判定可能给出六方，但残差偏大。VELAS 在报告里会同时给出“最可能的对称类”和“与理想约束的最大相对偏差”，当偏差超过 5% 时，我在报告里会加一句“建议按正交对称处理”。多花一点计算量，换来的是后续所有解释的可靠性。
