状态空间模型
狀態(tài)空間模型是動(dòng)態(tài)時(shí)域模型，以隱含著的時(shí)間為自變量。狀態(tài)空間模型包括兩個(gè)模型：
一是狀態(tài)方程模型，反映動(dòng)態(tài)系統(tǒng)在輸入變量作用下在某時(shí)刻所轉(zhuǎn)移到的狀態(tài)；
二是輸出或量測方程模型，它將系統(tǒng)在某時(shí)刻的輸出和系統(tǒng)的狀態(tài)及輸入變量聯(lián)系起來。
狀態(tài)空間模型分類
??狀態(tài)空間模型按所受影響因素的不同分為：
（1）確定性狀態(tài)空間模型
（2）隨機(jī)性狀態(tài)空間模型
??狀態(tài)空間模型按數(shù)值形式分為：
（1）離散空間狀態(tài)模型
（2）連續(xù)空間狀態(tài)模型
狀態(tài)空間模型按所描述的動(dòng)態(tài)系統(tǒng)可以分為：
（1）線性的與非線性的
（2）時(shí)變的與時(shí)不變的
5.jpg (105.47 KB, 下載次數(shù): 15)
下載附件
二、系統(tǒng)的狀態(tài)空間
離散事件隨機(jī)性系統(tǒng)的概念是系統(tǒng)理論中最基本的概念。
離散事件隨機(jī)性系統(tǒng)的狀態(tài)，是指系統(tǒng)內(nèi)部的可能運(yùn)動(dòng)狀態(tài)和可能儲(chǔ)能狀態(tài)。系統(tǒng)在k=k0時(shí)刻的狀態(tài)，是在k<k0時(shí)以系統(tǒng)內(nèi)部儲(chǔ)能的積累結(jié)果，并在k=k0時(shí)以系統(tǒng)要素儲(chǔ)能的方式表現(xiàn)出來，還將影響系統(tǒng)在k>k0時(shí)的外部行為。
用隨機(jī)向量序列來描述系統(tǒng)在任一時(shí)刻的狀態(tài)向量，稱為狀態(tài)向量法，也稱為狀態(tài)空間法。狀態(tài)向量表示為：
6.jpg (10.91 KB, 下載次數(shù): 15)
下載附件
其中
7.jpg (6.88 KB, 下載次數(shù): 14)
下載附件
（k=1，2，1…，n）為第i個(gè)狀態(tài)向量。三、系統(tǒng)的輸入輸出
輸入輸出狀態(tài)概念
引入狀態(tài)向量是為了對系統(tǒng)內(nèi)部結(jié)構(gòu)進(jìn)行數(shù)學(xué)描述，但在許多情況下，系統(tǒng)的狀態(tài)是無法直接量測到的，有時(shí)甚至全部不能量測到。在實(shí)際工作中，能量測到的量之勢系統(tǒng)的輸入與輸出。
輸入輸出變量表示
系統(tǒng)的輸入也是隨時(shí)間而變的一組變量，表示為：
8.jpg (9.92 KB, 下載次數(shù): 13)
下載附件
稱為輸入向量，其分量
9.jpg (8.07 KB, 下載次數(shù): 12)
下載附件
（i=1，2，…，r)稱為輸入變量。系統(tǒng)所受的隨機(jī)干擾也是隨時(shí)間而變的一組變量，表示為
10.jpg (11.01 KB, 下載次數(shù): 15)
下載附件
稱為系統(tǒng)的動(dòng)態(tài)模型噪聲，它是系統(tǒng)的一種特殊輸入向量。系統(tǒng)的輸出也是隨時(shí)間而變的一組變量，表示為：
11.jpg (9.29 KB, 下載次數(shù): 11)
下載附件
稱為輸出向量，其分量? ???（i=1，2，…，m)稱為輸入變量。量測系統(tǒng)也會(huì)受到隨機(jī)噪聲的污染，表示為：
12.jpg (8.93 KB, 下載次數(shù): 15)
下載附件
稱為系統(tǒng)的量測噪聲。
四、狀態(tài)空間模型
狀態(tài)空間模型定義
狀態(tài)空間模型是描述動(dòng)態(tài)系統(tǒng)的完整模型，它表達(dá)了由于輸入引起系統(tǒng)內(nèi)部狀態(tài)的變化，并由此使輸出發(fā)生的變化。
狀態(tài)空間模型的不同形式
如，線性時(shí)不變模型的狀態(tài)方程可表示為：
13.jpg (5.84 KB, 下載次數(shù): 10)
下載附件
輸出方程為：
14.jpg (3.36 KB, 下載次數(shù): 12)
下載附件
五、狀態(tài)空間模型的建立
? 例 1? ? 某養(yǎng)魚場為了反映池塘魚種的變化，請你幫助建立狀態(tài)空間模型。
解答： （1）取狀態(tài)向量X（k）為k時(shí)刻3個(gè)魚種的數(shù)量：
15.jpg (19.19 KB, 下載次數(shù): 13)
下載附件
輸入向量為：16.jpg (25.87 KB, 下載次數(shù): 14)
下載附件
（2）狀態(tài)轉(zhuǎn)移矩陣
17.jpg (8.44 KB, 下載次數(shù): 15)
下載附件
式中： p1，p2，p3為鯽魚、青魚和鯉魚的生長率，這里為p1=0.1，p2=0.13，p3=0.08。（3）輸入矩陣仍定為常陣
18.jpg (8.26 KB, 下載次數(shù): 11)
下載附件
（4）輸出矩陣或預(yù)測矩陣C為3×3維單位陣，這樣輸出向量或量測向量就等同于狀態(tài)向量，狀態(tài)空間模型：
19.jpg (7.91 KB, 下載次數(shù): 13)
下載附件
即：
20.jpg (33.42 KB, 下載次數(shù): 12)
下載附件
21.jpg (4.77 KB, 下載次數(shù): 14)
下載附件
22.jpg (13.19 KB, 下載次數(shù): 11)
下載附件
Kalman濾波簡介
Kalman濾波是一種線性濾波與預(yù)測方法，原文為：A New Approach to Linear Filtering and Prediction Problems。文章推導(dǎo)很復(fù)雜，看了一半就看不下去了，既然不能透徹理解其原理，但總可以通過實(shí)驗(yàn)來理解其具體的使用方法。
Kalman濾波分為2個(gè)步驟，預(yù)測(predict)和校正(correct)。預(yù)測是基于上一時(shí)刻狀態(tài)估計(jì)當(dāng)前時(shí)刻狀態(tài)，而校正則是綜合當(dāng)前時(shí)刻的估計(jì)狀態(tài)與觀測狀態(tài)，估計(jì)出最優(yōu)的狀態(tài)。預(yù)測與校正的過程如下：
預(yù)測：
1.png (6.58 KB, 下載次數(shù): 12)
下載附件
校正：
2.png (11.75 KB, 下載次數(shù): 15)
下載附件
公式1是狀態(tài)預(yù)測，公式2是誤差矩陣預(yù)測，公式3是kalman增益計(jì)算，公式4是狀態(tài)校正，其輸出即是最終的kalman濾波結(jié)果，公式5是誤差矩陣更新。各變量說明如下表：
4.jpg (120.25 KB, 下載次數(shù): 12)
下載附件
算法實(shí)現(xiàn)與分析
Kalman濾波最復(fù)雜的計(jì)算應(yīng)該就是公式3中的矩陣求逆，考慮到實(shí)現(xiàn)的方便性，采用matlab來簡單實(shí)現(xiàn)，本文主要是分析kalman濾波中各個(gè)變量的作用和對濾波結(jié)果的影響。具體代碼如下：
| 1 2 3 4 5 6 7 8 9 | function filter = Kalman(filter) ??%predict ??predict_x = filter.A * filter.x + filter.B * filter.u; ??filter.P = filter.A * filter.P * filter.A' + filter.Q; ??%correct ??filter.K = filter.P * filter.H' / (filter.H * filter.P * filter.H' + filter.R); ??filter.x = predict_x + filter.K * (filter.z - filter.H * predict_x); ??filter.P = filter.P - filter.K * filter.H * filter.P; end |
在matlab中，kalman濾波實(shí)際上就是上面那5個(gè)公式，而難點(diǎn)卻是在測試代碼中針對不同問題各個(gè)變量的初始化上，下面來逐個(gè)分析。
1.建立模型，明確觀測量，系統(tǒng)狀態(tài)以及其轉(zhuǎn)移方程(下面這段公式太多，通過word寫好后截圖)
3.png (121.08 KB, 下載次數(shù): 13)
下載附件
2.初始化噪聲協(xié)方差矩陣
經(jīng)過上面一步，只有PQRK四個(gè)矩陣還未確定了。顯然增益矩陣K是不需要初始化的，P是誤差矩陣，初始化可以是一個(gè)隨機(jī)的矩陣或者0，只要經(jīng)過幾次的處理基本上就能調(diào)整到正常的水平，因此也就只會(huì)影響前面幾次的濾波結(jié)果。
Q和R分別是預(yù)測和觀測狀態(tài)協(xié)方差矩陣，一般可以簡單認(rèn)為系統(tǒng)狀態(tài)各維之間(即上面的a和b)相互獨(dú)立，那么Q和R就可以設(shè)置為對角陣。而這兩個(gè)對角線元素的大小將直接影響著濾波結(jié)果，若Q的元素遠(yuǎn)大于R的元素，則預(yù)測噪聲大，從而更相信觀測值，這樣可能使得kalman濾波結(jié)果與觀測值基本一致；反之，則更相信預(yù)測，kalman濾波結(jié)果會(huì)表現(xiàn)得比較規(guī)整和平滑；若二者接近，則濾波結(jié)果介于前面兩者之間，根據(jù)實(shí)驗(yàn)效果看也缺乏實(shí)際使用價(jià)值。
以上幾個(gè)矩陣確定后，對于狀態(tài)x，由于0時(shí)刻我們沒有任何關(guān)于該系統(tǒng)的知識(shí)，可以使用0時(shí)刻的測量值z0來初始x0，預(yù)測從k=1開始；也可以初始化-1時(shí)刻的狀態(tài)，當(dāng)然這個(gè)狀態(tài)實(shí)際是未知的，也就可隨機(jī)取。2種方式都可以，但使用0時(shí)刻測量值來初始化狀態(tài)，可以使得前面幾次預(yù)測更準(zhǔn)確。
3.實(shí)驗(yàn)分析
首先使用下面代碼生成一組數(shù)據(jù)存在z.mat中：
[C] 純文本查看 復(fù)制代碼 ?| 1 2 3 4 5 6 7 8 | interval = pi/18; t = 1:interval:100*pi; len = size(t, 2); a = t + 4 * (rand(1,len)-0.5); b = t .* sin(t/10) +? 10 * (rand(1,len)-0.5); z = [a; b]; save('z.mat','z'); plot(z(1,:),z(2,:),'o') |
可以看出其近似為一條振幅不斷增大的正弦曲線疊加一個(gè)隨機(jī)噪聲。繪制出來如下：
4.png (23.9 KB, 下載次數(shù): 10)
下載附件
如果使用上面推導(dǎo)的恒定狀態(tài)系統(tǒng)模型，代碼與實(shí)驗(yàn)結(jié)果如下：
[C] 純文本查看 復(fù)制代碼 ?| 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 | clear close all clc dim_observe = 2;????????? %觀測值維數(shù) n = dim_observe;? %狀態(tài)維數(shù)，觀測狀態(tài)每個(gè)維度都有1個(gè)速度，故需乘2 filter.A = eye(n);%[1,0,1,0;0,1,0,1;0,0,1,0;0,0,0,1]; filter.B = 0; filter.u = 0; filter.P = eye(n); filter.K = zeros(n); filter.H = eye(n);%[1,0,0,0;0,1,0,0]; cQ = 1e-8; cR = 1e-2; filter.Q = eye(n) * cQ;??????? %這里簡單設(shè)置Q和R對角線元素都相等，設(shè)為不等亦可 filter.R = eye(dim_observe) * cR; filter.x = zeros(n,1); %初始狀態(tài)x0 load('z.mat'); figure(1),subplot(2,2,1), t = 1; out = []; for i=1:size(z,2) ??filter.z = z(:,i); ??filter = Kalman(filter); ??plot(filter.x(1),filter.x(2),? 'r*');hold on??????? ??plot(filter.z(1),filter.z(2),? 'bo');??????? hold on ??out=[out filter.x]; %???????? pause(.5) end figure(1), str = sprintf('cQ = %e, cR = %e', cQ, cR); title(str) %畫局部放大 subplot(2,2,2), plot(out(1,:),out(2,:),? 'r*');hold on??????? plot(z(1,:),z(2,:),? 'bo');??????? hold on axis([120 170 80 200]) |
5.png (32.82 KB, 下載次數(shù): 13)
下載附件
可以看出濾波結(jié)果完全滯后于測量數(shù)據(jù)，其根本原因在于建立的模型存在問題。
如果采用上面推導(dǎo)的物體運(yùn)動(dòng)模型則只需要修改部分代碼，主要是矩陣A和H，以及其他矩陣對應(yīng)的維數(shù)，具體如下：
[C] 純文本查看 復(fù)制代碼 ?| 1 2 3 4 5 6 7 8 | dim_observe = 2;????? %觀測值維數(shù) n = 2 * dim_observe;? %狀態(tài)維數(shù)，觀測狀態(tài)每個(gè)維度都有1個(gè)速度，故需乘2 filter.A = [1,0,1,0;0,1,0,1;0,0,1,0;0,0,0,1]; filter.B = 0; filter.u = 0; filter.P = eye(n); filter.K = zeros(n); filter.H = [1,0,0,0;0,1,0,0]; |
運(yùn)行結(jié)果如下圖，藍(lán)色為觀測數(shù)據(jù)，紅色為kalman濾波數(shù)據(jù)，右側(cè)為局部放大圖。可以看出經(jīng)過濾波后的數(shù)據(jù)相當(dāng)平滑，這里Q和R中元素的量級分別為cQ和cR，下圖結(jié)果可以看到cR比cQ多了6個(gè)數(shù)量級。
6.png (39.06 KB, 下載次數(shù): 14)
下載附件
（1）
增加幾組結(jié)果用于對比分析，對于的cQ和cR見圖的標(biāo)題。
7.png (39.18 KB, 下載次數(shù): 14)
下載附件
（2）
8.png (40.65 KB, 下載次數(shù): 13)
下載附件
（3）
9.png (41.63 KB, 下載次數(shù): 14)
下載附件
（4）
10.png (41.2 KB, 下載次數(shù): 13)
下載附件
（5）
11.png (41.14 KB, 下載次數(shù): 11)
下載附件
（6）
首先看圖1和2，cR與cQ大小均相差了3個(gè)數(shù)量級，而二者的比值相同，則kalman濾波結(jié)果相同。
再看圖2~圖6，cR/cQ在不斷減小，kalman濾波結(jié)果的平滑性也在不斷降低，到圖5和6中，濾波結(jié)果完全和觀測值相同，說明此時(shí)kalman濾波已經(jīng)完全相信觀測值了。原因在于cR/cQ過小，系統(tǒng)認(rèn)為預(yù)測噪聲的方差很大，不值得信賴，而觀測值的噪聲方差小，可信度高。
根據(jù)上面的實(shí)驗(yàn)結(jié)果，可以看出Kalman濾波應(yīng)用中的幾個(gè)問題：
1.模型建立的正確性從根本上決定了濾波效果的正確性。
上面使用物體靜止模型進(jìn)行濾波，結(jié)果完全不對，而使用勻速運(yùn)動(dòng)模型則能達(dá)到較好的效果。從根本上講，上面的數(shù)據(jù)也不是勻速運(yùn)動(dòng)的，為何結(jié)果會(huì)基本正確？看看第一個(gè)使用靜止模型的濾波結(jié)果，雖然我們假定了物體是靜止的，但由于觀測數(shù)據(jù)的作用，kalman濾波結(jié)果也會(huì)有相應(yīng)的運(yùn)動(dòng)而不是完全靜止，也就是說濾波器在不停地修正這個(gè)狀態(tài)，而在勻速運(yùn)動(dòng)模型中，物體的速度我們認(rèn)為不變，但同樣地kalman濾波器也會(huì)不停地修正這個(gè)速度，濾波器中計(jì)算的速度實(shí)質(zhì)的偏離了真實(shí)速度的，因此最終也會(huì)有相應(yīng)的偏差，不過這個(gè)偏差在我們?nèi)菰S范圍內(nèi)，也就可以大膽使用了。
如果能確定物體是勻變速直線運(yùn)動(dòng)，使用相應(yīng)帶加速度的模型會(huì)得到更準(zhǔn)確的效果。但是越嚴(yán)格的模型其適用范圍也相應(yīng)越小。
2.影響濾波結(jié)果平滑性的因素是cR/cQ，這個(gè)值反映了我們對于預(yù)測和觀測值的信任程度；其值越大則越相信預(yù)測結(jié)果，濾波結(jié)果平滑性好；反之則越相信觀測結(jié)果，濾波結(jié)果越偏向于觀測值。一般我們使用kalman濾波器是為了能平滑數(shù)據(jù)的波動(dòng)，因此應(yīng)盡量保證cR/cQ稍大，上面的測試結(jié)果該值在1e4以上數(shù)據(jù)較為平滑。
總結(jié)
- 上一篇： 如何用PADS打开AD的PCB文件？
- 下一篇： iOS录音功能