5D Library › 嗓音生物力学、非线性源-滤波耦合与听觉-运动反馈的数学模型 · English
第四章 音区、双质量模型与非线性分岔
音区是振动机制的切换
歌者所说的「胸声」「头声」「假声」，在嗓音科学里越来越多地用喉部振动机制来描述：M0（气泡音）、M1（通常对应胸声或常规说话）、M2（头声或假声）、M3（哨音）。
证据
鲁博、昂里克与卡斯特朗戈根据电声门图与声学数据主张，音区首先应当按喉部振动机制来定义：在 M1 中声带较厚、内侧接触面较大、闭合相较长；在 M2 中声带变薄、主要是边缘振动、闭合相较短。从一个机制到另一个机制的转换往往伴随基频的突然跳变，而且上行与下行时发生跳变的音高不同（迟滞），受过训练的歌者可以在一定音域内用两种机制唱同一个音 (Roubeau et al. 2009)。从生物力学看，机制的切换主要取决于甲杓肌（声带肌，「体」）与环甲肌（拉长声带）之间的平衡，以及由此改变的覆盖层厚度与刚度 (Titze 1994)。
突然跳变、迟滞、同一参数下两种状态并存——这些都是非线性动力系统的典型特征。第二章已经看到，声道负载可以在 f0 越过 F1 时促发失稳；本章看声源本身的非线性。
双质量模型
证据
石坂与弗拉纳根 1972 年提出用两个由弹簧连接的质量块代表每侧声带的下部和上部：下部质量较大，上部较小，两者之间的耦合使上部的运动落后于下部，于是自然产生了第一章所说的相位差，模型能够自激振荡并合成出逼真的浊音 (Ishizaka and Flanagan 1972)。施泰内克与赫策尔在此基础上引入左右不对称参数 ：右侧声带的质量除以 、刚度乘以 ，使右侧的固有频率是左侧的 倍；他们画出了随不对称与声门下压变化的分岔图，报告了左右同步、次谐波与不规则振动等多种状态 (Steinecke and Herzel 1995)。
下面的程序实现这个不对称双质量模型（标准参数，单位为克、厘米、毫秒），在 800 Pa 的肺压下改变 ，统计左右两侧下部质量每秒的振动次数。
推导
# Asymmetric two-mass model (Ishizaka & Flanagan 1972; asymmetry Q of Steinecke & Herzel
# 1995), units g, cm, ms. Each fold: lower mass m1, upper m2, springs k1, k2, coupling kc,
# collision springs 3k, damping ratios 0.1 / 0.6. Areas a_i = a0 + l (x_il + x_ir); lower
# masses feel P1 = Ps [1 - (a_min / a1)^2] while open. Right fold: masses / Q, springs x
# Q.
import math
M1, M2, K1, K2, KC, D1, LEN, A0 = 0.125, 0.025, 0.08, 0.008, 0.025, 0.25, 1.4, 0.05
def make(q):
return [(M1 / s, M2 / s, K1 * s, K2 * s, KC * s, 3 * K1 * s, 3 * K2 * s,
0.2 * math.sqrt(M1 * K1), 1.2 * math.sqrt(M2 * K2)) for s in (1.0, q)]
def deriv(y, par, ps):
x1l, x2l, x1r, x2r, v1l, v2l, v1r, v2r = y
a1, a2 = A0 + LEN * (x1l + x1r), A0 + LEN * (x2l + x2r)
if a1 > 0 and a2 > 0:
p1 = ps * (1 - (min(a1, a2) / a1) ** 2)
else:
p1 = 0.0 if a1 <= 0 else ps
acc = []
for (m1, m2, k1, k2, kc, c1, c2, r1, r2), xa, xb, va, vb in zip(
par, (x1l, x1r), (x2l, x2r), (v1l, v1r), (v2l, v2r)):
col1 = c1 * a1 / (2 * LEN) if a1 < 0 else 0.0
col2 = c2 * a2 / (2 * LEN) if a2 < 0 else 0.0
acc += [(LEN * D1 * p1 - r1 * va - k1 * xa - kc * (xa - xb) - col1) / m1,
(-r2 * vb - k2 * xb - kc * (xb - xa) - col2) / m2]
return (v1l, v2l, v1r, v2r, acc[0], acc[1], acc[2], acc[3])
def run(q, ps=0.008, y0=(0.02, 0.0, -0.005, 0.0), t_end=1500.0, t_skip=500.0, dt=0.04):
par, y, t, xl, xr = make(q), tuple(y0) + (0.0,) * 4, 0.0, [], []
for _ in range(int(t_end / dt)): # classical Runge-Kutta
k1 = deriv(y, par, ps)
k2 = deriv(tuple(a + dt / 2 * b for a, b in zip(y, k1)), par, ps)
k3 = deriv(tuple(a + dt / 2 * b for a, b in zip(y, k2)), par, ps)
k4 = deriv(tuple(a + dt * b for a, b in zip(y, k3)), par, ps)
y = tuple(a + dt / 6 * (b + 2 * c + 2 * d + e)
for a, b, c, d, e in zip(y, k1, k2, k3, k4))
t += dt
if t > t_skip:
xl.append(y[0]); xr.append(y[2])
pk = lambda x: [x[i] for i in range(1, len(x) - 1) if x[i - 1] < x[i] >= x[i + 1]]
span, pl = (t_end - t_skip) / 1000.0, pk(xl)
return len(pl) / span, len(pk(xr)) / span, len({round(v, 3) for v in pl})
if __name__ == "__main__":
print("lung pressure 0.008 g/(cm ms^2) = 800 Pa; Q = right/left natural frequency")
print(" Q | left fold Hz | right fold Hz | heights* | regime")
for q in (1.0, 0.8, 0.7, 0.6, 0.55, 0.5, 0.45, 0.4, 0.35, 0.3):
fl, fr, n = run(q)
r = round(fl / fr)
reg = ("1:1, folds in step" if abs(fl / fr - 1) < 0.02 else
f"{r}:1 locked, subharmonic" if abs(fl / fr - r) < 0.05
else "not locked (irregular or quasi-periodic)")
print(f" {q:4.2f} | {fl:12.1f} | {fr:13.1f} | {n:6d} | {reg}")
print("* distinct peak heights of the left lower mass, to 0.001 cm (1 = period-1)")
# lung pressure 0.008 g/(cm ms^2) = 800 Pa; Q = right/left natural frequency
# Q | left fold Hz | right fold Hz | heights* | regime
# 1.00 | 136.0 | 136.0 | 1 | 1:1, folds in step
# 0.80 | 118.0 | 118.0 | 1 | 1:1, folds in step
# 0.70 | 101.0 | 101.0 | 1 | 1:1, folds in step
# 0.60 | 83.0 | 83.0 | 1 | 1:1, folds in step
# 0.55 | 144.0 | 72.0 | 2 | 2:1 locked, subharmonic
# 0.50 | 132.0 | 66.0 | 2 | 2:1 locked, subharmonic
# 0.45 | 120.0 | 60.0 | 2 | 2:1 locked, subharmonic
# 0.40 | 159.0 | 53.0 | 3 | 3:1 locked, subharmonic
# 0.35 | 133.0 | 44.0 | 3 | 3:1 locked, subharmonic
# 0.30 | 115.0 | 38.0 | 3 | 3:1 locked, subharmonic
# * distinct peak heights of the left lower mass, to 0.001 cm (1 = period-1)
结果呈现出清楚的阶梯。 从 1 降到 0.6 时，两侧仍然同步（1:1），只是共同的频率从 136 Hz 降到 83 Hz：较慢的右侧把较快的左侧「拖住」了。 在 0.45–0.55 之间，同步被打破，系统进入 2:1 锁频：左侧回到约 120–144 Hz，右侧只以它一半的频率振动；声门的开合于是每两个左侧周期才重复一次，频谱中出现左侧频率一半处的次谐波。 为 0.4 及以下时，变为 3:1。表中最后一栏（左侧下部质量振幅峰值的种数，精确到 0.001 cm）与此对应：1:1 时为 1，2:1 时为 2，3:1 时为 3，说明左侧的运动也随锁频比重复。这里的「左右」只是模型中的标签，不对应任何具体的病变。
在这组参数下，我们的扫描只找到了周期性的锁频状态；施泰内克与赫策尔在其他参数区域还报告了不规则振动。程序从略不对称的初始状态出发，去掉前 0.5 秒过渡后分析 1 秒；换一组参数或初始状态，阶梯的位置可能移动。
从周期到混沌：一个通用的例子
「倍周期通向混沌」是非线性动力学最著名的道路之一。为了看清它，我们暂时离开声带，用一个最常被引用的非线性振子——受迫的杜芬方程。每个驱动周期采样一次振子的位置（频闪映射）：周期解只给一个点，倍周期给两个点，以此类推；混沌时样本散成一片，而且两个初始几乎相同的轨道会以指数速度分开（最大李雅普诺夫指数为正）。
推导
# Part A. A generic nonlinear oscillator (NOT a vocal-fold model): forced Duffing equation
# x'' + 0.3 x' - x + x^3 = g cos(1.2 t), sampled once per forcing period. 1 distinct state
# = period-1, 2 = period-doubled, ...; many = chaos if the Lyapunov exponent is positive.
# Part B. A period-doubled pulse train (no noise at all) seen through Book 10's measures.
import cmath, math
DMP, W = 0.3, 1.2
def rk4(t, x, v, g, dt):
f = lambda t, x, v: (v, -DMP * v + x - x ** 3 + g * math.cos(W * t))
a1, b1 = f(t, x, v)
a2, b2 = f(t + dt / 2, x + dt / 2 * a1, v + dt / 2 * b1)
a3, b3 = f(t + dt / 2, x + dt / 2 * a2, v + dt / 2 * b2)
a4, b4 = f(t + dt, x + dt * a3, v + dt * b3)
return (x + dt / 6 * (a1 + 2 * a2 + 2 * a3 + a4),
v + dt / 6 * (b1 + 2 * b2 + 2 * b3 + b4))
def strobe(g, n_skip=200, n_keep=64, steps=200):
dt = 2 * math.pi / W / steps
x, v, y, u, t, out, lyap = 1.0, 0.0, 1.0 + 1e-8, 0.0, 0.0, [], 0.0
for k in range(n_skip + n_keep):
for _ in range(steps):
(x, v), (y, u), t = rk4(t, x, v, g, dt), rk4(t, y, u, g, dt), t + dt
# twin orbit, renormalised each period
d = math.hypot(y - x, u - v)
if k >= n_skip:
out.append(round(x, 3)); lyap += math.log(d / 1e-8)
y, u = x + (y - x) * 1e-8 / d, v + (u - v) * 1e-8 / d
return len(set(out)), lyap / n_keep
print("Part A: forced Duffing oscillator, stroboscopic samples")
print(" forcing g | distinct states | Lyapunov exponent per period | regime")
for g in (0.20, 0.27, 0.29, 0.30, 0.37, 0.50):
n, lam = strobe(g)
reg = f"period-{n}" if n <= 8 and lam < 0.01 else "chaotic" if lam > 0.01 else "?"
print(f" {g:.2f} | {n:15d} | {lam:+28.3f} | {reg}")
print("Part B: a 100 Hz pulse train whose pulses alternate in height (period doubling)")
def level(x, f): # one DFT bin, 8 kHz sampling
return abs(sum(v * cmath.exp(-2j * math.pi * f * k / 8000) for k, v in enumerate(x)))
# 8 kHz sampling, 80 samples per pulse
for ratio in (1.0, 0.9, 0.7):
amps = [1.0 if i % 2 == 0 else ratio for i in range(80)]
x = [a * math.sin(math.pi * k / 40) ** 2 if k < 40 else 0.0
for a in amps for k in range(80)]
shimmer = sum(abs(a - b) for a, b in zip(amps, amps[1:])) / 79 / (sum(amps) / 80)
sub = level(x, 50.0) / level(x, 100.0)
sub_txt = f"{20 * math.log10(sub):6.1f} dB re 100 Hz" if sub > 1e-9 else " none"
print(f" every 2nd pulse x {ratio:.1f}: shimmer {shimmer:5.1%}, jitter 0.0 %, "
f"50 Hz {sub_txt}")
# Part A: forced Duffing oscillator, stroboscopic samples
# forcing g | distinct states | Lyapunov exponent per period | regime
# 0.20 | 1 | -0.784 | period-1
# 0.27 | 2 | -0.314 | period-2
# 0.29 | 4 | -0.187 | period-4
# 0.30 | 61 | +0.398 | chaotic
# 0.37 | 5 | -0.267 | period-5
# 0.50 | 60 | +0.600 | chaotic
# Part B: a 100 Hz pulse train whose pulses alternate in height (period doubling)
# every 2nd pulse x 1.0: shimmer 0.0%, jitter 0.0 %, 50 Hz none
# every 2nd pulse x 0.9: shimmer 10.5%, jitter 0.0 %, 50 Hz -24.5 dB re 100 Hz
# every 2nd pulse x 0.7: shimmer 35.3%, jitter 0.0 %, 50 Hz -14.0 dB re 100 Hz驱动幅度 时为周期 1，0.27 时倍周期（2 个状态），0.29 时周期 4，0.30 时样本散成 61 个不同的值、李雅普诺夫指数为 +0.40，进入混沌；继续加大驱动，0.37 处出现一个周期 5 的窗口，0.50 时又是混沌（+0.60）。图三右图画出了整个过程。这说明「混沌」是一个精确的数学性质（对初始条件的敏感依赖），不是「乱」的同义词；而且混沌区里夹着周期窗口，参数的微小改变可以使振动从混沌回到周期。
程序的第二部分回到嗓音测量。一串 100 Hz 的脉冲，每隔一个脉冲高度变为 0.7 倍——这是最简单的倍周期，没有任何噪声。第 10 本第六章用的局部振幅微扰（shimmer）却高达 35.3%，基频微扰为 0，频谱中出现比 100 Hz 低约 14 dB 的 50 Hz 成分；高度变为 0.9 倍时，振幅微扰仍有 10.5%。也就是说，基频微扰与振幅微扰无法区分「有规律的倍周期」与「随机的不规则」，它们只适合描述接近周期的嗓音。
证据
赫策尔、贝里、提策与萨利赫把非线性动力学的方法（窄带语谱图中的次谐波与分岔、吸引子重构等）用于嗓音障碍患者的录音分析，认为许多病理嗓音中的不规则性来自声带振动本身的非线性现象，而不只是叠加的噪声，并主张用这些概念对嗓音分类 (Herzel et al. 1994)。贝里等人在离体喉实验中系统改变声门下压与左右不对称（内收或拉长的不对称），观察到包括单侧声带振动在内的多种振动状态；许多失稳与分岔出现在不同状态重叠的参数区域，部分录音的频谱与文献中描述的音区相符 (Berry et al. 1996)。
注
混沌吸引子不是声带息肉。 框架中「混沌吸引子 = 声带息肉」的写法不能成立。研究者在某些病理嗓音中观察到次谐波与不规则振动，但同样的现象也出现在健康嗓音的音区转换处、在离体喉与模型的不对称条件下，以及有意的歌唱技巧中；而声带息肉、结节、麻痹等不同病变都可能产生相似的声学表现，也可能没有这些表现。一个声学现象不能告诉我们原因。本书的任何程序、图表或核查表都不能用于判断某人是否有声带病变；持续的声音问题请看耳鼻喉科或嗓音专科医生。
有意的「不规则」：呼麦与摇滚失真音
证据
蒙古卡尔格拉（kargyraa）唱法。 林德斯塔德等人用高速影像、声学频谱与逆滤波研究这种低沉的「喉音」，发现室带（假声带）与声带同时振动，室带大约每两个声带周期闭合一次，形成倍周期，听到的音高约为声带频率的一半 (Lindestad et al. 2001)。贝利、昂里克与佩洛松对倍周期发声中的声带与室带振动做了生理描述，并建立了室带参与的空气动力学模型 (Bailly et al. 2010)。摇滚「失真」音（dist）。 赞格·博尔奇等人发现，失真音中的非周期性来自气流脉冲幅度的低频调制，这种调制由声门上方黏膜的振动造成——它部分阻挡气道，降低部分脉冲的幅度——而声带本身仍大致周期地振动；他们推测声门上黏膜的振动由较高声门下压带来的大气流驱动 (Zangger Borch et al. 2004)。
这些研究说明，次谐波与倍周期也可以是受控的艺术手段。它们的机制各不相同：卡尔格拉主要是室带参与的倍周期，摇滚失真音是声门上结构对声源的调制。另一类常被称为「呼麦」的泛音唱法，通过精细调节共振峰突出单个泛音，机制又不相同，本书不展开。框架中把呼麦与撕裂音一概写成「室带参与、次谐波」，只对其中一部分成立。
假说 · 尚未验证
混声迟滞假说（尚未确立）。经过系统的「混声」训练，歌者在缓慢上下滑音中从 M1 跳到 M2 与从 M2 跳回 M1 的音高之差（迟滞宽度）会缩小，跳变处的基频突变幅度也会减小。可以检验的预测：在一项对未受训练成人的随机对照训练研究中，训练组迟滞宽度的缩小显著大于对照组（例如只做音准练习）。如果两组的迟滞宽度变化没有差别，假说就被推翻。
本章小结
音区可以理解为喉部振动机制（M0–M3）的切换，跳变与迟滞是非线性系统的典型特征。双质量模型用两块耦合的质量自然地产生上下缘相位差；加入左右不对称后，示意模型从 1:1 同步走到 2:1、3:1 锁频，出现次谐波。用通用的杜芬振子可以看到经典的倍周期通向混沌，混沌是对初始条件敏感依赖的精确性质。有规律的倍周期会使振幅微扰高达 35%，说明这类指标只适合接近周期的嗓音。某些病理嗓音中观察到不规则振动，但它不能用于诊断；卡尔格拉与摇滚失真音则是对非线性的有意利用。下一章从声源转向耳朵：歌者怎样通过一个有延迟的回路听自己。
参考文献
Fan, Y. (2026). Vocal Biomechanics, Non-Linear Source-Filter Coupling and Auditory-Motor Feedback: From Titze's Self-Sustained Vocal Fold Oscillation and Vocal-Tract Inertance to the Lombard Effect, Read Beside Classical Chinese Vocal Aesthetics [嗓音生物力学、非线性源-滤波耦合与听觉-运动反馈的数学模型] (Version 1.0.0). 5D Library. https://doi.org/10.5281/zenodo.23770084 · Licensed under Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International (CC BY-NC-ND 4.0).