← 听觉系统模拟 · 体感系统模拟 →
简化的半规管，没有耳石膜。 让我们考虑半规管（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 .} 现在，我们考虑一个恒定幅度为 ω {\displaystyle \omega } θ ˙ ( t ) {\displaystyle {\dot {\theta }}(t)}
δ = 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}
要找到这个传递函数的极点，我们必须确定分母等于 0 的 s 值
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 中的模块signal
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”是一个简单的增益系数。“带帽”变量表示内部估计。 在数学上，具有高增益的负反馈具有有趣的性质，即它实际上可以反转负反馈回路中的传递函数：如果 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)