← 听觉系统模拟 · 体感系统模拟 →
简化的半规管，没有耳石膜。 让我们考虑半规管 (SCC) 的机械描述。在下面的描述中，我们将做出非常强烈的简化假设。这里的目标仅仅是理解半规管背后的基本机械原理。
我们做出的第一个强烈简化是，半规管可以被建模为一个外半径为 R，内半径为 r 的圆形管子。（对于适当的流体力学推导，参见 (Damiano 和 Rabbitt 1996) 以及 Obrist (2005)）。这个管子充满了内淋巴液。
半规管的方向可以在给定的坐标系中用一个向量 n → {\displaystyle {\vec {n}}}
θ {\displaystyle \theta } θ ˙ ≡ d θ d t {\displaystyle {\dot {\theta }}\equiv {\frac {d\theta }{dt}}} θ ¨ ≡ d 2 θ d t 2 {\displaystyle {\ddot {\theta }}\equiv {\frac {d^{2}\theta }{dt^{2}}}} ϕ {\displaystyle \phi } δ = θ − ϕ {\displaystyle \delta =\theta -\phi } 请注意，所有这些变量都是标量。我们利用了管子的角速度可以看作是头部实际角速度向量 ω → {\displaystyle {\vec {\omega }}} n → {\displaystyle {\vec {n}}}
θ ˙ = ω → ⋅ n → {\displaystyle {\dot {\theta }}={\vec {\omega }}\cdot {\vec {n}}} 其中点表示标准标量积。
为了描述内淋巴液的运动，考虑一个自由漂浮的活塞，其密度与内淋巴液相同。两个力作用在这个系统上
惯性矩 I ϕ ¨ {\displaystyle I{\ddot {\phi }}}
粘性力矩 B δ ˙ {\displaystyle B{\dot {\delta }}} 这给出了运动方程
I ϕ ¨ = B δ ˙ {\displaystyle I{\ddot {\phi }}=B{\dot {\delta }}} 代入 ϕ = θ − δ {\displaystyle \phi =\theta -\delta }
θ ˙ = δ ˙ + B I δ . {\displaystyle {\dot {\theta }}={\dot {\delta }}+{\frac {B}{I}}\delta .} 现在，让我们考虑速度阶跃 θ ˙ ( t ) {\displaystyle {\dot {\theta }}(t)} ω {\displaystyle \omega }
δ = I B ω ⋅ ( 1 − e − B I t ) {\displaystyle \delta ={\frac {I}{B}}\omega \cdot (1-e^{-{\frac {B}{I}}t})} 并且对于 t ≫ I B {\displaystyle t\gg {\frac {I}{B}}}
δ ≈ I B ω {\displaystyle \delta \approx {\frac {I}{B}}\omega } 现在，让我们推导出时间常数 T 1 ≡ I B {\displaystyle T_{1}\equiv {\frac {I}{B}}} r ≪ R {\displaystyle r\ll R}
I = m l 2 ≈ 2 ρ π 2 r 2 R 3 . {\displaystyle I=ml^{2}\approx 2\rho \pi ^{2}r^{2}R^{3}.} 根据泊肃叶-哈根方程，细管中层流速度为 v 时的力 F 为
F = 8 V ¯ η l r 2 {\displaystyle F={\frac {8{\bar {V}}\eta l}{r^{2}}}} 其中 V ¯ = r 2 π v {\displaystyle {\bar {V}}=r^{2}\pi v} η {\displaystyle \eta } l = 2 π R {\displaystyle l=2\pi R}
在扭矩 M = F ⋅ R {\displaystyle M=F\cdot R} Ω = v R {\displaystyle \Omega ={\frac {v}{R}}}
B = M Ω = 16 η π 2 R 3 {\displaystyle B={\frac {M}{\Omega }}=16\eta \pi ^{2}R^{3}} 最后，得到时间常数 T 1 {\displaystyle T_{1}}
T 1 = I B = δ r 2 8 η {\displaystyle T_{1}={\frac {I}{B}}={\frac {\delta r^{2}}{8\eta }}} 对于人类平衡系统，用实验获得的参数替换变量，得到时间常数 T 1 {\displaystyle T_{1}} ≈ {\displaystyle \approx }
G ≡ δ ω = I B = T 1 {\displaystyle G\equiv {\frac {\delta }{\omega }}={\frac {I}{B}}=T_{1}} 杯形的影响。 到目前为止，我们的讨论还没有包括杯形在 SCC 中的作用：杯形充当一个弹性膜，通过角加速度发生位移。通过其弹性，杯形将系统恢复到其静止位置。杯形的弹性在运动方程中增加了额外的弹性项。如果将其考虑在内，此方程变为
θ ¨ = δ ¨ + B I δ ˙ + K I δ {\displaystyle {\ddot {\theta }}={\ddot {\delta }}+{\frac {B}{I}}{\dot {\delta }}+{\frac {K}{I}}\delta }
求解此类微分方程的一种优雅方法是拉普拉斯变换 。拉普拉斯变换将微分方程转化为代数方程：如果信号 x(t) 的拉普拉斯变换用 X(s) 表示，则时间导数的拉普拉斯变换为
d x ( t ) d t → L a p l a c e T r a n s f o r m s ⋅ X ( s ) − x ( 0 ) {\displaystyle {\frac {dx(t)}{dt}}{\xrightarrow {LaplaceTransform}}s\cdot X(s)-x(0)} 项 x(0) 详细说明了起始条件，并且通常可以通过适当选择参考位置将其设置为零。因此，拉普拉斯变换为
s 2 θ ~ = s 2 δ ~ + B I s δ ~ + K I δ ~ {\displaystyle s^{2}{\tilde {\theta }}=s^{2}{\tilde {\delta }}+{\frac {B}{I}}s{\tilde {\delta }}+{\frac {K}{I}}{\tilde {\delta }}} 其中“~”表示拉普拉斯变换后的变量。根据上面的 T 1 {\displaystyle T_{1}} T 2 {\displaystyle T_{2}}
T 2 = B K {\displaystyle T_{2}={\frac {B}{K}}} 我们得到
δ ~ θ ~ = T 1 s 2 T 1 s 2 + s + 1 T 2 {\displaystyle {\frac {\tilde {\delta }}{\tilde {\theta }}}={\frac {T_{1}s^{2}}{T_{1}s^{2}+s+{\frac {1}{T_{2}}}}}} 对于人类， T 2 = B / K {\displaystyle T_{2}=B/K}
为了找到这个传递函数的极点，我们必须确定对于哪些 s 值，分母等于 0
s 1 , 2 = 1 T 1 ( − 1 ± 1 − 4 T 1 T 2 ) {\displaystyle s_{1,2}={\frac {1}{T_{1}}}{\Big (}-1\pm {\sqrt {1-4{\frac {T_{1}}{T_{2}}}}}{\Big )}} 由于 T 2 ≫ T 1 {\displaystyle T_{2}\gg T_{1}}
1 − x ≈ 1 − x 2 f o r x ≪ 1 {\displaystyle {\sqrt {1-x}}\approx 1-{\frac {x}{2}}forx\ll 1} 我们得到
s 1 ≈ − 1 T 1 , a n d s 2 ≈ − 1 T 2 {\displaystyle s_{1}\approx -{\frac {1}{T_{1}}},ands_{2}\approx -{\frac {1}{T_{2}}}} 通常我们感兴趣的是耳石器位移 δ {\displaystyle \delta } θ ˙ ≡ s θ ~ {\displaystyle {\dot {\theta }}\equiv s{\tilde {\theta }}}
δ ~ s θ ~ ( s ) = T 1 T 2 s ( T 1 s + 1 ) ( T 2 s + 1 ) {\displaystyle {\frac {\tilde {\delta }}{s{\tilde {\theta }}}}(s)={\frac {T_{1}T_{2}s}{(T_{1}s+1)(T_{2}s+1)}}} 对于典型的头部运动 (0.2 Hz < f < 20 Hz)，系统增益近似恒定。换句话说，对于典型的头部运动，耳石器位移与头部角速度成正比！
头部速度函数的耳石器位移的波特图，其中 T1 = 0.01 秒，T2 = 5 秒，放大系数为 (T1+ T2)/ (T1* T2)，以获得中心频率的增益约为 0。 对于线性时不变系统（LTI系统），输入和输出在频域中具有简单的关系
O u t ( s ) = G ( s ) ∗ I n ( s ) {\displaystyle Out(s)=G(s)*In(s)} 其中传递函数G(s) 可以用代数函数表示
G ( s ) = n u m ( s ) d e n ( s ) = n ( 0 ) ∗ s 0 + n ( 1 ) ∗ s 1 + n ( 2 ) ∗ s 2 + . . . d ( 0 ) ∗ s 0 + d ( 1 ) ∗ s 1 + d ( 2 ) ∗ s 2 + . . . {\displaystyle G(s)={\frac {num(s)}{den(s)}}={\frac {n(0)*{{s}^{0}}+n(1)*{{s}^{1}}+n(2)*{{s}^{2}}+...}{d(0)*{{s}^{0}}+d(1)*{{s}^{1}}+d(2)*{{s}^{2}}+...}}} 换句话说，指定分子 (n) 和分母 (d) 的系数将唯一地表征传递函数。这种表示法被一些计算工具用于模拟这种系统对给定输入的响应。
不同的工具可以用来模拟这样的系统。例如，一个时间常数为7秒的低通滤波器对1秒的输入阶跃的响应具有以下传递函数
G ( s ) = 1 7 s + 1 {\displaystyle G(s)={\frac {1}{7s+1}}} 并且可以如下模拟
使用 Simulink 进行低通滤波器的阶跃响应模拟。 如果您在命令行上工作，您可以使用 MATLAB 的控制系统工具箱 或 Python 包SciPy 的信号 模块
MATLAB 控制系统工具箱
% Define the transfer function
num = [ 1 ];
tau = 7 ;
den = [ tau , 1 ];
mySystem = tf ( num , den )
% Generate an input step
t = 0 : 0.1 : 30 ;
inSignal = zeros ( size ( t ));
inSignal ( t >= 1 ) = 1 ;
% Simulate and show the output
[ outSignal , tSim ] = lsim ( mySystem , inSignal , t );
plot ( t , inSignal , tSim , outSignal );
Python - SciPy
# Import required packages
import numpy as np
import scipy.signal as ss
import matplotlib.pylab as mp
# Define transfer function
num = [ 1 ]
tau = 7
den = [ tau , 1 ]
mySystem = ss . lti ( num , den )
# Generate inSignal
t = np . arange ( 0 , 30 , 0.1 )
inSignal = np . zeros ( t . size )
inSignal [ t >= 1 ] = 1
# Simulate and plot outSignal
tout , outSignal , xout = ss . lsim ( mySystem , inSignal , t )
mp . plot ( t , inSignal , tout , outSignal )
mp . show ()
现在考虑耳石器官的力学。由于它们由复杂的粘弹性材料组成，并具有弯曲的形状，因此它们的力学无法用分析工具来描述。但是，它们的运动可以用有限元法进行数值模拟。由此，所考虑的体积被分成许多小的体积元素，并且对于每个元素，物理方程被解析函数近似。
有限元模拟：使用小的有限元来构建一个机械模型；这里例如是椭圆囊。 这里我们将只展示粘弹性耳石材料的物理方程。每种弹性材料的运动必须服从柯西运动方程
ρ ∂ 2 u i ∂ t 2 = ρ B i + ∑ j ∂ T i j ∂ x j {\displaystyle \rho {\frac {\partial ^{2}u_{i}}{\partial t^{2}}}=\rho B_{i}+\sum _{j}{\frac {\partial T_{ij}}{\partial x_{j}}}}
其中 ρ {\displaystyle \rho } u i {\displaystyle u_{i}} B i {\displaystyle B_{i}} T i j {\displaystyle T_{ij}} x j {\displaystyle x_{j}}
对于线性弹性、各向同性材料，柯西应力张量 由下式给出
T i j = λ e δ i j + 2 μ E i j {\displaystyle T_{ij}=\lambda e\delta _{ij}+2\mu E_{ij}} 其中 λ {\displaystyle \lambda } μ {\displaystyle \mu } 拉梅常数 ； μ {\displaystyle \mu } e = d i v ( u → ) {\displaystyle e=div({\vec {u}})} E i j {\displaystyle E_{ij}}
E i j = 1 2 ( ∂ u i ∂ x j + ∂ u j ∂ x i ) . {\displaystyle E_{ij}={\frac {1}{2}}{\Big (}{\frac {\partial u_{i}}{\partial x_{j}}}+{\frac {\partial u_{j}}{\partial x_{i}}}{\Big )}.} 这导致了纳维运动方程。
ρ ∂ 2 u i ∂ t 2 = ρ B i + ( λ + μ ) ∂ e ∂ x i + μ ∑ j ∂ 2 u i ∂ x j 2 {\displaystyle \rho {\frac {\partial ^{2}u_{i}}{\partial t^{2}}}=\rho B_{i}+(\lambda +\mu ){\frac {\partial e}{\partial x_{i}}}+\mu \sum _{j}{\frac {\partial ^{2}u_{i}}{\partial x_{j}^{2}}}} 此方程适用于纯弹性、各向同性材料，并且可以用有限元技术求解。找到此方程中出现的机械参数的典型程序如下：当将材料的圆柱形样本置于应变下时，杨氏模量 E 表征长度的变化，而泊松比 ν {\displaystyle \nu } 表征直径的同时减少。拉梅常数 λ {\displaystyle \lambda } μ {\displaystyle \mu } ν {\displaystyle \nu }
E = μ ( 3 λ + 2 μ ) λ + μ {\displaystyle E={\frac {\mu (3\lambda +2\mu )}{\lambda +\mu }}} 以及
ν = λ 2 ( λ + μ ) {\displaystyle \nu ={\frac {\lambda }{2(\lambda +\mu )}}} 前庭信息的中央加工显着影响了对空间中感知的方向和运动。脑干中相应的处理过程通常可以用控制系统工具有效地建模。作为一个具体的例子，我们展示了如何模拟速度存储 的影响。
速度存储 的概念基于以下实验发现：当我们从围绕地球垂直轴的持续旋转中突然停止时，耳石器会被减速所偏转，但以约 5 秒的时间常数恢复到静止状态。然而，感知到的旋转持续的时间长得多，并且以更长的时间常数下降，通常在 15 到 20 秒之间。
前庭模型：蓝线描述了耳石器的偏转，作为对速度阶跃的响应，建模为时间常数为 5 秒的高通滤波器。绿线表示角速度的内部估计，通过一个负反馈回路中的耳石器响应内部模型和一个前馈增益因子 2 获得。 在附图中，运河对角速度刺激 ω 的响应由传递函数 C 建模，这里是一个时间常数为 5 秒的简单高通滤波器。（运河响应由耳石器的偏转决定，并且近似地与神经放电率成正比。）为了模拟时间常数的增加，我们假设中枢前庭系统有一个运河传递函数的内部模型 ， C ^ {\displaystyle {\hat {C}}} ω ^ {\displaystyle {\hat {\omega }}} k 设置为 2，模型的输出很好地再现了时间常数的增加。相应的 Python 代码可以在 [ 1]
值得注意的是，这个反馈回路可以在生理上得到解释：我们知道左右前庭核之间有很强的联系。如果这些连接被切断，感知到的旋转的时间常数会降低到半规管的外周时间常数。
中枢前庭处理通常可以用控制系统模型来描述。这里 “omega” 是头部速度， “C” 是半规管的传递函数， “k” 是一个简单的增益因子。带有 “hat” 的变量表示内部估计。 在数学上，高增益负反馈具有一个有趣的特性，它可以实际上反转负反馈回路中的传递函数：如果 k>>1 ，并且如果运河传递函数的内部模型与实际传递函数相似，则估计的角速度对应于实际的角速度。
ω ^ = ( ω C − ω ^ C ^ ) k ω ^ ( 1 + C ^ k ) = ω C k ω ^ ω = C 1 / k + C ^ → i f C ≈ C ^ k >> 1 1 {\displaystyle {\begin{aligned}&{\hat {\omega }}=(\omega C-{\hat {\omega }}{\hat {C}})k\\&{\hat {\omega }}(1+{\hat {C}}k)=\omega Ck\\&{\frac {\hat {\omega }}{\omega }}={\frac {C}{1/k+{\hat {C}}}}\,\,{\xrightarrow[{if\,C\approx {\hat {C}}}]{k>>1}}1\end{aligned}}}
![{\displaystyle {\begin{aligned}&{\hat {\omega }}=(\omega C-{\hat {\omega }}{\hat {C}})k\\&{\hat {\omega }}(1+{\hat {C}}k)=\omega Ck\\&{\frac {\hat {\omega }}{\omega }}={\frac {C}{1/k+{\hat {C}}}}\,\,{\xrightarrow[{if\,C\approx {\hat {C}}}]{k>>1}}1\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9e94ee6fa51d798b237236c8cc1dc926473505b8)