簡(jiǎn)介卡爾曼濾波是處理含噪動(dòng)態(tài)系統(tǒng)狀態(tài)估計(jì)與多傳感器數(shù)據(jù)融合的經(jīng)典遞歸算法面向信號(hào)處理、控制工程以及機(jī)器學(xué)習(xí)方向的學(xué)習(xí)者與開(kāi)發(fā)者提供一套基于MATLAB的卡爾曼濾波數(shù)據(jù)融合演示代碼。資源包共2個(gè)文件包含1個(gè)m腳本和1個(gè)txt配套說(shuō)明文件整體大小僅2KB腳本用于實(shí)現(xiàn)多源觀測(cè)信息的融合與遞推更新文本文件則可作為輸入數(shù)據(jù)或參數(shù)配置參考適合在理解算法原理后動(dòng)手運(yùn)行與調(diào)參驗(yàn)證。資源描述對(duì)基本原理亦有清晰梳理例如系統(tǒng)狀態(tài)服從高斯分布、通過(guò)預(yù)測(cè)和更新兩個(gè)步驟相接合獲得最優(yōu)估計(jì)并將應(yīng)用場(chǎng)景延伸至模式識(shí)別序列數(shù)據(jù)跟蹤、函數(shù)逼近與在線參數(shù)估計(jì)方便讀者從單一示例擴(kuò)展到實(shí)際項(xiàng)目。雖然壓縮包體積很小但演示鏈路完整能夠覆蓋從模型構(gòu)建、觀測(cè)輸入到狀態(tài)輸出的全流程。已有655人學(xué)習(xí)瀏覽通過(guò)運(yùn)行該腳本可以直觀觀察預(yù)測(cè)與更新迭代對(duì)狀態(tài)估計(jì)的影響理解噪聲協(xié)方差與權(quán)重分配的作用為后續(xù)實(shí)現(xiàn)更復(fù)雜的多傳感器融合系統(tǒng)打下扎實(shí)基礎(chǔ)。1. 兩個(gè)傳感器數(shù)據(jù)打架時(shí)卡爾曼濾波做的不是平均而是建模做 datsfusion 這類多源融合任務(wù)時(shí)我第一次拿到兩路同時(shí)測(cè)同一個(gè)位置的傳感器數(shù)據(jù)直覺(jué)反應(yīng)是加權(quán)平均噪聲大的少乘一點(diǎn)噪聲小的多乘一點(diǎn)。結(jié)果一跑就發(fā)現(xiàn)靜態(tài)場(chǎng)景勉強(qiáng)能用目標(biāo)一加速輸出立刻滯后半拍傳感器噪聲稍微變一下權(quán)重又要重新調(diào)。卡爾曼濾波把這個(gè)問(wèn)題的解法換了一個(gè)方向——它不直接決定“信誰(shuí)多少”而是給兩個(gè)傳感器各建一個(gè)噪聲模型再按模型的不確定度動(dòng)態(tài)生成融合權(quán)重。哪個(gè)傳感器此刻更可信是由協(xié)方差遞推算出來(lái)的不是預(yù)先寫死的。這篇文章就沿著這個(gè)思路把狀態(tài)方程、Q/R 協(xié)方差設(shè)定、最小可復(fù)現(xiàn)代碼、擴(kuò)展到目標(biāo)跟蹤的 EKF 路徑以及現(xiàn)場(chǎng)調(diào)參的驗(yàn)證手段一路講到底。適合做定位、導(dǎo)航、多傳感器融合以及卡爾曼濾波目標(biāo)跟蹤的工程師新手能跟著跑通代碼熟手可以重點(diǎn)看后兩章的參數(shù)邊界和坑點(diǎn)。2. 數(shù)據(jù)融合用卡爾曼濾波先把狀態(tài)方程寫對(duì)噪聲協(xié)方差才有意義2.1 狀態(tài)空間表達(dá)式把多傳感器觀測(cè)統(tǒng)一到同一個(gè)狀態(tài)向量里卡爾曼濾波處理數(shù)據(jù)融合第一步不是寫濾波公式而是把物理問(wèn)題翻譯成狀態(tài)空間表達(dá)式。狀態(tài)方程x_k A x_{k-1} B u_k w_k 測(cè)量方程z_k H x_k v_kx 是狀態(tài)向量它描述系統(tǒng)在某一時(shí)刻的完整內(nèi)部狀態(tài)。做數(shù)據(jù)融合時(shí)不同傳感器觀測(cè)的是同一個(gè)狀態(tài)的不同投影GPS 直接給位置里程計(jì)給速度積分加速度計(jì)給比力。如果把它們各自獨(dú)立地拿來(lái)做估計(jì)得到的是好幾個(gè)互相不一致的答案放進(jìn)卡爾曼濾波框架后所有傳感器都對(duì)應(yīng)到同一個(gè) x 上只是各自的 H 不同。H 的物理含義是“狀態(tài)到測(cè)量的映射”GPS 的 H 是位置那幾列取 1里程計(jì)的 H 是速度那幾列取 1。w 和 v 分別代表過(guò)程噪聲與測(cè)量噪聲它們各自的協(xié)方差矩陣 Q 和 R才是后續(xù)卡爾曼增益計(jì)算里決定融合權(quán)重的核心。常見(jiàn)做法是先在紙上把狀態(tài)向量列出來(lái)再逐個(gè)傳感器寫 H最后才碰代碼。我見(jiàn)過(guò)不少融合結(jié)果發(fā)散的項(xiàng)目根源都不是濾波公式寫錯(cuò)而是狀態(tài)向量順序和 H 矩陣列沒(méi)對(duì)齊導(dǎo)致預(yù)測(cè)和更新在融合兩個(gè)毫不相干的量。2.2 Q 和 R 的物理語(yǔ)義卡爾曼增益背后的融合權(quán)重分配卡爾曼增益 K 的表達(dá)式是K P H^T (H P H^T R)^-1K 的數(shù)值決定了最終輸出在“預(yù)測(cè)值”和“測(cè)量值”之間取哪里。P 是狀態(tài)協(xié)方差矩陣它反映當(dāng)前對(duì)狀態(tài)估計(jì)的不確定度R 是測(cè)量噪聲協(xié)方差反映傳感器讀數(shù)本身的可信度。當(dāng) R 相對(duì)于 P 很大時(shí)K 趨近于 0濾波器更相信模型預(yù)測(cè)當(dāng) R 相對(duì)于 P 很小時(shí)K 趨近于 1濾波器更相信傳感器。這就是數(shù)據(jù)融合用卡爾曼濾波的本質(zhì)它不做固定權(quán)重的平均而是讓權(quán)重隨不確定度實(shí)時(shí)變化。某一時(shí)刻 GPS 信號(hào)抖動(dòng)劇烈R 變大融合結(jié)果自然偏向里程計(jì)下一秒 GPS 恢復(fù)權(quán)重又自動(dòng)切回來(lái)。Q 描述的是模型本身的不可信程度。如果系統(tǒng)實(shí)際在加速而狀態(tài)方程用的是勻速模型那么 Q 就要留出足夠的余量來(lái)容納這種未建模的動(dòng)態(tài)Q 設(shè)得太小濾波器會(huì)過(guò)度信任預(yù)測(cè)表現(xiàn)為輸出平滑但滯后嚴(yán)重。Q 設(shè)得太大融合結(jié)果又會(huì)被過(guò)程噪聲淹沒(méi)噪聲抑制效果退化成一階低通濾波甚至更差。提示卡爾曼濾波和常見(jiàn)的一階低通濾波的差別就在這個(gè) K 是自適應(yīng)計(jì)算的。一階低通濾波的系數(shù)是常數(shù)卡爾曼的“系數(shù)”是隨 P 和 R 不斷更新的變量。2.3 定義融合系統(tǒng)的狀態(tài)向量與測(cè)量矩陣一張表說(shuō)清維度對(duì)應(yīng)以 datsfusion 里最常見(jiàn)的二維定位融合為例狀態(tài)向量一般取四個(gè)量x 方向位置、x 方向速度、y 方向位置、y 方向速度。GPS 輸出位置里程計(jì)輸出速度H 矩陣就要按這個(gè)順序去對(duì)齊。狀態(tài)向量成員物理含義對(duì)應(yīng)傳感器測(cè)量方程中的體現(xiàn)p_xx 軸位置GPSz_gps p_x v_gpsv_xx 軸速度里程計(jì)/IMUz_odom v_x v_odomp_yy 軸位置GPSz_gps p_y v_gpsv_yy 軸速度里程計(jì)/IMUz_odom v_y v_odom狀態(tài)向量一旦含速度A 矩陣自然就是勻速模型的形式位置列里帶一個(gè) dt速度列保持 1。很多初學(xué)者把狀態(tài)向量只定義為位置用兩個(gè)位置傳感器做融合也能跑通但融合結(jié)果沒(méi)有速度估計(jì)后續(xù)做目標(biāo)跟蹤時(shí)就沒(méi)有預(yù)測(cè)能力每次都要等測(cè)量到了才能動(dòng)。既然做的是數(shù)據(jù)融合能多狀態(tài)量就多狀態(tài)量這樣 Q 矩陣能把速度的變化過(guò)程也建模進(jìn)去預(yù)測(cè)出來(lái)的軌跡才連貫。H 矩陣的維度是“測(cè)量數(shù) × 狀態(tài)維數(shù)”。GPS 和里程計(jì)同時(shí)到達(dá)時(shí)測(cè)量向量是兩維H 是 2×4 的矩陣每一行對(duì)應(yīng)一個(gè)傳感器的觀測(cè)方程。寫完第一步就是檢查矩陣維度能不能相乘這一步錯(cuò)了后面代碼跑起來(lái)全是維度報(bào)錯(cuò)但報(bào)錯(cuò)信息往往會(huì)把你引向完全無(wú)關(guān)的地方。3. 用卡爾曼濾波做數(shù)據(jù)融合一個(gè)能直接跑的最小 Python 實(shí)現(xiàn)3.1 最小卡爾曼濾波類predict 與 update 兩個(gè)核心方法先寫一個(gè)一維的單變量卡爾曼濾波器類再看它怎么融合兩路傳感器。單變量場(chǎng)景里所有矩陣都退化成標(biāo)量邏輯最直觀。import numpy as np class KalmanFilter1D: 一維卡爾曼濾波用于說(shuō)明 predict / update 的完整流程 def __init__(self, q0.01, p01.0): self.q q # 過(guò)程噪聲方差模型未建模動(dòng)態(tài)的余量 self.p p0 # 初始協(xié)方差給大一點(diǎn)可以加速收斂 self.x 0.0 # 初始狀態(tài)估計(jì) def predict(self, dt1.0): # 先驗(yàn)估計(jì)狀態(tài)不變協(xié)方差累加過(guò)程噪聲 self.p self.p self.q * dt def update(self, z, r): # z 是傳感器讀數(shù)r 是該傳感器的測(cè)量噪聲方差 k self.p / (self.p r) # 卡爾曼增益 self.x self.x k * (z - self.x) # 用新息修正狀態(tài) self.p (1 - k) * self.p # 更新后協(xié)方差 return self.x這里卡爾曼增益 k 的表達(dá)式是一維化的形式K P / (P R)。當(dāng) R 遠(yuǎn)大于 Pk 趨近于 0測(cè)量修正量很小當(dāng) R 遠(yuǎn)小于 Pk 趨近于 1狀態(tài)幾乎完全跟隨測(cè)量。update 里的(z - self.x)是新息表示“測(cè)量值與當(dāng)前預(yù)測(cè)的差值”卡爾曼濾波的全部修正動(dòng)作都建立在這個(gè)差值上。r 參數(shù)按傳感器分別傳入這給了數(shù)據(jù)融合一個(gè)天然的結(jié)構(gòu)同一個(gè)濾波器不同傳感器各自帶自己的噪聲方差。3.2 雙傳感器融合的最小場(chǎng)景GPS 與里程計(jì)同時(shí)輸出位置現(xiàn)在構(gòu)造一個(gè)仿真場(chǎng)景目標(biāo)真實(shí)位置是恒定值 5.0GPS 噪聲標(biāo)準(zhǔn)差 2.0里程計(jì)噪聲標(biāo)準(zhǔn)差 0.5。注意這里用標(biāo)準(zhǔn)差構(gòu)造數(shù)據(jù)傳入濾波器時(shí)要把標(biāo)準(zhǔn)差平方成方差。np.random.seed(42) true_pos 5.0 n 100 # 兩路傳感器觀測(cè)同一位置噪聲水平差異很大 gps true_pos np.random.normal(0, 2.0, n) # GPS大噪聲 odom true_pos np.random.normal(0, 0.5, n) # 里程計(jì)小噪聲 kf KalmanFilter1D(q0.01, p01.0) states [] for t in range(n): kf.predict(dt1.0) kf.update(gps[t], r2.0**2) # GPS 的 R 為 4.0 kf.update(odom[t], r0.5**2) # 里程計(jì) R 為 0.25 states.append(kf.x) # 對(duì)比純 GPS 均值誤差 vs 融合結(jié)果誤差 gps_err np.abs(np.mean(gps) - true_pos) fusion_err np.abs(states[-1] - true_pos) print(GPS 均值誤差: %.3f, 卡爾曼融合誤差: %.3f % (gps_err, fusion_err))同一時(shí)刻兩路測(cè)量到達(dá)可以連續(xù)調(diào)用兩次 update。第一次用 GPS 更新協(xié)方差 P 會(huì)縮小第二次用里程計(jì)更新是在 GPS 更新完的基礎(chǔ)上繼續(xù)修正。順序本身不影響最終結(jié)果因?yàn)榭柭f推滿足結(jié)合律但每次 update 前不要忘記先 predict 一次。r 參數(shù)要按傳感器分別傳入GPS 的 R 是 4.0里程計(jì)是 0.25后者權(quán)重天然更大。最終融合誤差通常遠(yuǎn)小于純 GPS 的均值誤差這就是把測(cè)量噪聲模型化之后得到的收益。注意R 傳方差還是標(biāo)準(zhǔn)差是新手最容易踩的坑。傳感器手冊(cè)給的一般是“精度 ±0.5m”這是標(biāo)準(zhǔn)差傳進(jìn)公式時(shí)要平方成 0.25。不平方也能跑但濾波結(jié)果對(duì)噪聲的抑制程度會(huì)整體失真。3.3 單變量到多維融合的改動(dòng)邊界同一套公式換矩陣維度一維場(chǎng)景跑通后二維多維融合的代碼結(jié)構(gòu)不用變只把標(biāo)量換成矩陣。狀態(tài)向量 x 變成 4×1A 變成 4×4P 變成 4×4H 按傳感器測(cè)量維度取對(duì)應(yīng)的行。predict 里做F P F.T Qupdate 里先算S H P H.T R再算K P H.T np.linalg.inv(S)。建議在代碼里把矩陣維度寫清楚后先用單位矩陣做一次冒煙測(cè)試狀態(tài)向量全取 1觀測(cè)全取 1跑幾步看輸出是否發(fā)散。發(fā)散就逐項(xiàng)檢查 A、H 的維度配對(duì)而不是直接懷疑濾波公式。矩陣融合和多傳感器融合里 90% 的運(yùn)行時(shí)錯(cuò)誤都是維度錯(cuò)位一維代碼很難暴露這類問(wèn)題。4. 擴(kuò)展卡爾曼濾波與目標(biāo)跟蹤非線性融合場(chǎng)景的升級(jí)路徑4.1 為什么目標(biāo)跟蹤里要用擴(kuò)展卡爾曼濾波而不是線性卡爾曼線性卡爾曼的前提是狀態(tài)轉(zhuǎn)移和測(cè)量方程都是線性的。但目標(biāo)跟蹤場(chǎng)景里傳感器輸出的常常不是位置坐標(biāo)而是距離和角度比如雷達(dá)測(cè)出(r, θ)換算成直角坐標(biāo)是px r * cos(θ), py r * sin(θ)。這個(gè)映射關(guān)系天然帶三角函數(shù)不是線性方程。強(qiáng)行用線性卡爾曼的結(jié)果是測(cè)量更新這一步的 H 矩陣沒(méi)法表達(dá)真實(shí)的映射融合輸出會(huì)在目標(biāo)轉(zhuǎn)彎時(shí)出現(xiàn)明顯偏差。擴(kuò)展卡爾曼濾波EKF的處理方式是對(duì)非線性函數(shù)做一階泰勒展開(kāi)用雅可比矩陣代替原來(lái)的線性 H 矩陣。每一時(shí)刻的 H 都不同它依賴當(dāng)前狀態(tài)估計(jì)值所以要在每次更新時(shí)重新計(jì)算。這個(gè)“重新算 H”的動(dòng)作就是 EKF 比線性卡爾曼多出來(lái)的核心計(jì)算量。4.2 雷達(dá)測(cè)距測(cè)角場(chǎng)景的 EKF 預(yù)測(cè)與更新代碼狀態(tài)仍然取[px, py, vx, vy]過(guò)程模型繼續(xù)用勻速模型測(cè)量模型改成非線性函數(shù)并顯式寫出雅可比矩陣。import numpy as np def ekf_predict(x, P, dt, Q): # 勻速模型的狀態(tài)轉(zhuǎn)移矩陣 F np.array([ [1, 0, dt, 0], [0, 1, 0, dt], [0, 0, 1, 0], [0, 0, 0, 1] ]) x F x P F P F.T Q return x, P def ekf_update(x, P, z, R): px, py, vx, vy x r np.sqrt(px**2 py**2) # 雅可比矩陣測(cè)量方程對(duì)狀態(tài)求偏導(dǎo) H np.array([ [px / r, py / r, 0, 0], [-py / (r**2), px / (r**2), 0, 0] ]) # 預(yù)測(cè)測(cè)量值 z_pred np.array([r, np.arctan2(py, px)]) # 新息協(xié)方差、增益、狀態(tài)與協(xié)方差更新 S H P H.T R K P H.T np.linalg.inv(S) innovation z - z_pred # 角度新息要?dú)w一化到 [-pi, pi]否則角度跳變會(huì)擊穿濾波器 innovation[1] (innovation[1] np.pi) % (2 * np.pi) - np.pi x x K innovation P (np.eye(4) - K H) P return x, P, innovation雅可比矩陣的第一行是距離 r 對(duì)位置分量的偏導(dǎo)第二行是角度 θ 對(duì)位置分量的偏導(dǎo)。速度分量在測(cè)量方程里不出現(xiàn)對(duì)應(yīng)位置全為 0。角度新息的歸一化是目標(biāo)跟蹤里必須寫的細(xì)節(jié)角度從 179 度跳到 -179 度數(shù)值上變化了 358 度但真實(shí)角度變化只有 2 度。不做歸一化這個(gè)虛假的大新息會(huì)把濾波器狀態(tài)直接沖偏。4.3 融合發(fā)散前看這三個(gè)指標(biāo)新息、協(xié)方差對(duì)角線、誤差包絡(luò)EKF 跑起來(lái)之后不要只盯著最終輸出曲線看融合是否發(fā)散應(yīng)該用三個(gè)量化指標(biāo)預(yù)判。新息的均值應(yīng)該接近 0方差應(yīng)該接近 S 矩陣的對(duì)角線。如果新息長(zhǎng)期偏向一側(cè)說(shuō)明模型有系統(tǒng)偏差比如勻速模型跟不上轉(zhuǎn)彎目標(biāo)如果新息方差遠(yuǎn)大于 S 的預(yù)期說(shuō)明 R 設(shè)得太小濾波器過(guò)度相信測(cè)量。P 矩陣對(duì)角線代表每個(gè)狀態(tài)量的估計(jì)方差理想情況下應(yīng)該逐步收斂并保持穩(wěn)定。P 持續(xù)增長(zhǎng)不下降說(shuō)明更新環(huán)節(jié)沒(méi)有起到修正作用通常是 H 算錯(cuò)或者 R 設(shè)置過(guò)大P 快速收斂到極小值往往是過(guò)度自信后續(xù)一旦出現(xiàn)異常測(cè)量濾波器根本來(lái)不及反應(yīng)。誤差包絡(luò)是實(shí)際測(cè)量殘差與理論方差的比值。工程上常用一個(gè)移動(dòng)窗口統(tǒng)計(jì)實(shí)際新息的標(biāo)準(zhǔn)差與 S 的理論值做對(duì)比比值在 1.5 倍以內(nèi)算正常超過(guò) 3 倍基本可以判定濾波器模型失真需要回頭檢查 Q 和 R 的匹配度。5. 數(shù)據(jù)融合現(xiàn)場(chǎng)先調(diào)這三個(gè)參數(shù)再用新息驗(yàn)證濾波效果5.1 R 矩陣標(biāo)定用靜態(tài)數(shù)據(jù)算方差別靠手感R 是唯一能通過(guò)實(shí)驗(yàn)直接測(cè)出來(lái)的參數(shù)。做法是把傳感器放在靜止環(huán)境下采集 100 到 200 個(gè)數(shù)據(jù)點(diǎn)直接算方差。代碼就三行samples np.array([...]) # 靜止時(shí)采集的傳感器讀數(shù) r np.var(samples) # 測(cè)量噪聲方差靜態(tài)標(biāo)定得到的 R 是基線值實(shí)際運(yùn)行中還要根據(jù)場(chǎng)景微調(diào)。比如 GPS 在城市峽谷里噪聲會(huì)顯著增大這時(shí)候 R 按靜態(tài)值給濾波器會(huì)把嚴(yán)重偏離的測(cè)量當(dāng)成正常數(shù)據(jù)吸收進(jìn)去。工程上的常見(jiàn)做法是做一個(gè)簡(jiǎn)單的環(huán)境判斷當(dāng)新息連續(xù)多步超過(guò) 3 倍標(biāo)準(zhǔn)差時(shí)臨時(shí)把 R 放大 5 到 10 倍等新息回歸正常再切回來(lái)。5.2 Q 矩陣調(diào)參從滯后到跟手的現(xiàn)場(chǎng)步驟Q 沒(méi)有直接測(cè)量手段通?？啃Ч赐啤Ｕ{(diào)參步驟可以按下面這張表執(zhí)行現(xiàn)象原因判斷調(diào)整方向輸出平滑但滯后嚴(yán)重Q 偏小模型過(guò)度自信Q 增大 3~5 倍讓預(yù)測(cè)跟隨動(dòng)態(tài)輸出抖動(dòng)明顯噪聲沒(méi)壓住Q 偏大或 R 偏小先驗(yàn)證 R再減小 Q新息長(zhǎng)期不為零模型誤差或 Q 太小Q 增大或檢查狀態(tài)轉(zhuǎn)移方程P 收斂后不下降R 過(guò)大按標(biāo)定值重新給定 R調(diào)整 Q 時(shí)建議按對(duì)角元素逐項(xiàng)調(diào)不要整體縮放。位置對(duì)應(yīng)的 Q 和速度對(duì)應(yīng)的 Q 物理意義不同整體縮放一次往往會(huì)把原來(lái)合理的部分也弄壞。拿目標(biāo)跟蹤來(lái)說(shuō)位置過(guò)程噪聲來(lái)自目標(biāo)加速度的隨機(jī)變化速度過(guò)程噪聲還包含模型本身的截?cái)嗾`差兩項(xiàng)要給不同的數(shù)量級(jí)。5.3 新息序列驗(yàn)證與在線噪聲協(xié)方差估計(jì)的小技巧調(diào)完 Q 和 R 后用一個(gè)移動(dòng)窗口持續(xù)統(tǒng)計(jì)新息的均值和方差。均值應(yīng)該在 0 附近波動(dòng)窗口方差與 S 的理論值匹配。這里有一個(gè)實(shí)用的在線自適應(yīng)技巧用滑動(dòng)窗口估計(jì)實(shí)際新息方差sigma_hat如果它持續(xù)大于 S說(shuō)明 R 被低估了用兩者的比值去修正 Rwindow_size 50 innovation_history [] # 每步更新后追加新息 innovation_history.append(innovation[0]) if len(innovation_history) window_size: innovation_history.pop(0) # 窗口內(nèi)的實(shí)際新息標(biāo)準(zhǔn)差與理論 S 對(duì)比 sigma_hat np.std(innovation_history) if sigma_hat 2.0 * np.sqrt(S[0, 0]): R_adapted R_original * (sigma_hat / np.sqrt(S[0, 0]))**2這個(gè)比值修正法在噪聲緩變場(chǎng)景里非常有效能自動(dòng)適應(yīng)傳感器性能退化。但它對(duì)環(huán)境突變不敏感如果目標(biāo)突然大幅加速新息會(huì)驟增這時(shí)更合適的做法是臨時(shí)增大 Q 而不是調(diào)整 R。判斷依據(jù)很簡(jiǎn)單目標(biāo)動(dòng)態(tài)變了新息增大但協(xié)方差 P 正常傳感器壞了新息和 P 同步惡化。把這兩個(gè)信號(hào)區(qū)分開(kāi)融合系統(tǒng)才算真正具備現(xiàn)場(chǎng)可維護(hù)性。本文還有配套的精品資源點(diǎn)擊獲取