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