← 聽覺系統模擬 · 體感系統模擬 →
簡化的半規管，沒有耳石膜。 讓我們考慮半規管 (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)