1. مقدمة: 20 يوليو 1969، الهبوط في «بحر الهدوء» ومعجزة حاسوب بذاكرة $2\,\text{KB}$

في 20 يوليو 1969، وبينما كانت مركبة الهبوط القمرية Eagle التابعة لرحلة أبولو 11 تهبط نحو سطح بحر الهدوء، دوت في المقصورة إنذارات البرنامج 1202 و1201. على ارتفاع بضعة كيلومترات فقط فوق سطح القمر، اعتمد نيل أرمسترونغ وباز ألدرين كلياً على حاسوب التوجيه لأبولو (Apollo Guidance Computer - AGC) الذي طوره مختبر الأجهزة في معهد ماساتشوستس للتكنولوجيا (MIT).

بالمعايير الحديثة، كانت مواصفات حاسوب AGC محدودة للغاية: طول الكلمة 15 بت، وتردد الساعة $2.048\,\text{MHz}$، وذاكرة وصول عشوائي مغناطيسية قابلة للكتابة (RAM) لا تتجاوز $2,048$ كلمة (حوالي $2\,\text{KB}$ فقط متاحة للعمليات الحسابية اللحظية)، وذاكرة قراءة فقط (ROM) بسعة $36,864$ كلمة (حوالي $72\,\text{KB}$). علاوة على ذلك، كان كل مستشعر على متن المركبة الفضائية معيباً إذا استُخدم بمفرده:

1. وحدة القياس بالقصور الذاتي (IMU): تقيس الجيروسكوبات ومقاييس التسارع التسارع بشكل مستمر، ولكن عند مكاملة التسارع مرتين زمنياً لحساب الموقع (الملاحة التقديرية / Dead Reckoning)، فإن أي انحياز ضئيل يتراكم تربيعياً ($\propto t^2$) أو تكعيبياً ($\propto t^3$) مع مرور الوقت.
2. السدس البصري (Sextant): القياسات الفلكية اليدوية التي يجريها رواد الفضاء للنجوم والأفق لا تعاني من الانجراف الزمني، لكنها متقطعة زمنياً وتخضع لأخطاء التصويب البشري والتكميم.
3. رادار الهبوط وتتبع دوبلر: يعاني من ضجيج عشوائي عالي التردد بسبب فوهات القمر والغبار المنبعث والانعكاسات متعددة المسارات.

كيف استطاع حاسوب بذاكرة $2\,\text{KB}$ فقط أن يدمج في الوقت الفعلي توقعاً فيزيائياً ينجرف مع الزمن مع قياسات مستشعرات شديدة التشويش ليحصل على تقدير للموقع والسرعة أدق بكثير من أي مستشعر منفرد؟

الإجابة الرياضية التي أنقذت أبولو 11 هي مرشح كالمان (Kalman Filter) الذي نشره عام 1960 عالم الرياضيات ونظرية النظم المجري-الأمريكي رودولف إميل كالمان (Rudolf Emil Kálmán, 1930–2016).

1.1 ورقة كالمان عام 1960 والزيارة التاريخية لستانلي شميدت في مركز أميس التابع لناسا

في مارس 1960، نشر كالمان ورقته التاريخية “A New Approach to Linear Filtering and Prediction Problems”، متجاوزاً قيود مجال التردد في مرشح فينر (Wiener Filter) عبر صياغة مسألة التقدير في المجال الزمني باستخدام نموذج فضاء الحالة (State-Space Model) والإسقاط المتعامد.

في نوفمبر 1960، زار كالمان الباحث ستانلي شميدت (Stanley F. Schmidt) في مركز أبحاث أميس التابع لوكالة ناسا (NASA Ames Research Center). حتى ذلك الحين، كان تحديد المدارات يعتمد على طريقة المربعات الصغرى الدفعية (Batch Least Squares) التي ابتكرها كارل فريدريش غاوس عام 1801 لإعادة اكتشاف الكويكب سيريس. كانت الطريقة الدفعية تتطلب تخزين كامل سجل القياسات السابقة $\mathbf{Z}_{1:T} = \{\mathbf{z}_1, \dots, \mathbf{z}_T\}$ وعكس مصفوفات بتعقيد حسابي $O(T^3)$ وذاكرة $O(T)$، وهو أمر مستحيل تماماً على ذاكرة AGC البالغة $2\,\text{KB}$.

أثبت كالمان أنه بفضل خاصية ماركوف (Markov Property) من الدرجة الأولى، لا حاجة لتخزين أي قياسات سابقة على الإطلاق؛ إذ يُلخَّص التاريخ الكامل للمشاهدات دون أي فقدان للمعلومات في إحصاء كافٍ (Sufficient Statistic) يتكون فقط من متجه الحالة المقدرة الحالية $\hat{\mathbf{x}}_{k-1\mid k-1}$ ومصفوفة تباين الخطأ $\mathbf{P}_{k-1\mid k-1}$، مما يتطلب ذاكرة $O(1)$ في كل خطوة زمنية! سرعان ما وسّع شميدت هذه النظرية لتشمل الميكانيكا المدارية غير الخطية مبتكراً مرشح كالمان الموسع (EKF أو مرشح كالمان-شميدت)، بينما ابتكر جيمس بوتر (James E. Potter) في معهد MIT مرشح الجذر التربيعي (Square-Root Filter) الذي ضاعف الدقة العددية لحسابات الـ 15 بت ومنع ظهور قيم ذاتية سالبة لمصفوفة التباين.

2. المعجزة الرياضية في بُعد واحد (1D): لماذا ينتج عن ضرب توزيعين غاوسيين مشوشين تباين أصغر من أي مستشعر بمفرده؟

لنحلل الجوهر البايزي لمرشح كالمان في بُعد مكاني واحد ($x \in \mathbb{R}$):

1. يعطي التوقع الفيزيائي توزيعاً قبلياً غاوسياً $p(x) = \mathcal{N}(\mu_p, \sigma_p^2)$.
2. يعطي مستشعر GPS أو الرادار قراءة $z \equiv \mu_m$ بدالة معقولية غاوسية $p(z \mid x) = \mathcal{N}(\mu_m, \sigma_m^2)$.

وفقاً لقاعدة بايز (Bayes’ Rule)، تتناسب كثافة الاحتمال البعدية مع حاصل ضرب دالتي الكثافة الغاوسيتين:

عن طريق إكمال المربع (Completing the Square) في الأس ومطابقة معاملات $x^2$ و$x$:

نحصل على المتوسط البعدي $\mu_{\text{new}}$، والتباين البعدي $\sigma_{\text{new}}^2$، وكسب كالمان أحادي البعد $K \in [0, 1]$:

البصيرة الجوهرية (جمع الدقة / معلومات فيشر): مقلوب التباين $\tau = 1/\sigma^2$ يُسمى الأحكام أو الدقة (Precision / معلومات فيشر Fisher Information). عند دمج مصدرين غاوسيين مستقلين بايزياً، تُجمع دقتهما جمعاً خطياً ($\tau_{\text{new}} = \tau_p + \tau_m$، تماماً مثل جمع المواصلات الكهربائية لمقاومتين على التوازي). لذلك فإن عدم اليقين الناتج $\sigma_{\text{new}}^2$ يكون أصغر قطعا من عدم يقين التوقع وأصغر قطعاً من عدم يقين المستشعر بمفرده! على سبيل المثال، إذا كان $\sigma_p = 4\,\text{m}$ ($\sigma_p^2=16$) و$\sigma_m = 3\,\text{m}$ ($\sigma_m^2=9$)، فإن $K = 0.64$ وينخفض الانحراف المعياري المدمج إلى $\sigma_{\text{new}} = 2.4\,\text{m}$ فقط!

3. صياغة فضاء الحالة متعدد المتغيرات والاشتقاق الكامل لمعادلات مرشح كالمان الخمس

لنأخذ النظام الخطي العشوائي في الزمن المتقطع لمتجه الحالة $\mathbf{x}_k \in \mathbb{R}^n$ ومتجه القياس $\mathbf{z}_k \in \mathbb{R}^m$:

حيث ضجيج العملية $\mathbf{w}_k \sim \mathcal{N}(\mathbf{0}, \mathbf{Q}_k)$ وضجيج القياس $\mathbf{v}_k \sim \mathcal{N}(\mathbf{0}, \mathbf{R}_k)$ مستقلان إحصائياً.

الخطوة 1 — التوقع (Time Update / Chapman-Kolmogorov)

بأخذ التوقع الشرطي $\hat{\mathbf{x}}_{k\mid k-1} = \mathbb{E}[\mathbf{x}_k \mid \mathbf{Z}_{1:k-1}]$ وحساب تباين الخطأ القبلي $\tilde{\mathbf{x}}_{k\mid k-1} = \mathbf{F}_k \tilde{\mathbf{x}}_{k-1\mid k-1} + \mathbf{w}_k$، نحصل على المعادلتين الأولى والثانية:

الخطوة 2 — التحديث (Measurement Update / Bayes) واشتقاق $\mathbf{K}_k$

نعرف الابتكار (بواقي القياس / Innovation) $\tilde{\mathbf{y}}_k = \mathbf{z}_k - \mathbf{H}_k \hat{\mathbf{x}}_{k\mid k-1}$ ومصفوفة تباين الابتكار $\mathbf{S}_k = \mathbf{H}_k \mathbf{P}_{k\mid k-1} \mathbf{H}_k^T + \mathbf{R}_k$. معادلة التحديث الخطي هي:

الخطأ البعدي هو $\tilde{\mathbf{x}}_{k\mid k} = (\mathbf{I} - \mathbf{K}_k \mathbf{H}_k)\tilde{\mathbf{x}}_{k\mid k-1} - \mathbf{K}_k \mathbf{v}_k$، ومنه نحصل على صيغة جوزيف المستقرة للتباين (Joseph Stabilized Form) الصالحة لأي مصفوفة كسب $\mathbf{K}_k$:

لتقليل متوسط مربع الخطأ المتوقع $\mathbb{E}[\|\mathbf{x}_k - \hat{\mathbf{x}}_{k\mid k}\|_2^2] = \operatorname{Tr}(\mathbf{P}_{k\mid k})$، نشتق الأثر (Trace) بالنسبة للمصفوفة $\mathbf{K}_k$ ونساويه بالصفر:

وبالتعويض عن $\mathbf{K}_k$ الأمثل في مفكوك جوزيف يلغى الحدان الأخيران فنحصل على المعادلة الخامسة:

  flowchart TD
    Init["التهيئة الابتدائية<br/>x̂₀|₀, P₀|₀"] --> Predict["[1. خطوة التوقع / Time Update]<br/>x̂ₖ|ₖ₋₁ = Fₖ x̂ₖ₋₁|ₖ₋₁ + Bₖ uₖ<br/>Pₖ|ₖ₋₁ = Fₖ Pₖ₋₁|ₖ₋₁ Fₖᵀ + Qₖ"]
    Predict --> Innov["[حساب الابتكار وتباينه]<br/>ỹₖ = zₖ - Hₖ x̂ₖ|ₖ₋₁<br/>Sₖ = Hₖ Pₖ|ₖ₋₁ Hₖᵀ + Rₖ"]
    Innov --> Gain["[حساب كسب كالمان الأمثل]<br/>Kₖ = Pₖ|ₖ₋₁ Hₖᵀ Sₖ⁻¹"]
    Gain --> Update["[2. خطوة التحديث / Measurement Update]<br/>x̂ₖ|ₖ = x̂ₖ|ₖ₋₁ + Kₖ ỹₖ<br/>Pₖ|ₖ = (I - Kₖ Hₖ) Pₖ|ₖ₋₁"]
    Update --> Predict

4. أمثلية BLUE / MMSE، الثنائية، ومبدأ الفصل في التحكم الأمثل LQG

1. مبرهنة غاوس-ماركوف (BLUE) وأمثلية MMSE: تحت الضجيج الغاوسي، يحقق مرشح كالمان أقل متوسط خطأ تربيعي (MMSE) عالمياً بين كافة المقدرات الخطية وغير الخطية؛ وتحت الضجيج غير الغاوسي يظل أفضل مقدر خطي غير متحيز (BLUE).
2. مبدأ الفصل في تحكم LQG: بفضل الثنائية الرياضية بين القابلية للرصد (Observability) والقابلية للتحكم (Controllability)، ينفصل تصميم نظام التحكم العشوائي الأمثل Linear-Quadratic-Gaussian (LQG) تماماً إلى تصميم مستقل لمرشح كالمان ($\hat{\mathbf{x}}_{k\mid k}$) ومتحكم LQR الحتمي ($\mathbf{u}_k^* = -\mathbf{L}_k \hat{\mathbf{x}}_{k\mid k}$).

5. ما وراء الخطية: مرشحات EKF وUKF ومرشح الجسيمات وSLAM للقيادة الذاتية

للأنظمة غير الخطية ($\mathbf{x}_k = f(\mathbf{x}_{k-1}, \mathbf{u}_k) + \mathbf{w}_k, \mathbf{z}_k = h(\mathbf{x}_k) + \mathbf{v}_k$):

• مرشح كالمان الموسع (EKF): يقوم بالتقريب الخطي من الدرجة الأولى لمتسلسلة تايلور باستخدام مصفوفات جاكوبي $\mathbf{F}_k = \left.\frac{\partial f}{\partial \mathbf{x}}\right|_{\hat{\mathbf{x}}_{k-1\mid k-1}}$ و$\mathbf{H}_k = \left.\frac{\partial h}{\partial \mathbf{x}}\right|_{\hat{\mathbf{x}}_{k\mid k-1}}$.
• مرشح كالمان عديم الرائحة (UKF، Julier & Uhlmann 1997): يمرر $2n+1$ من نقاط سيغما (Sigma Points) الحتمية مباشرة عبر الدوال غير الخطية دون حساب أي مشتقات جاكوبية، محققاً دقة من الدرجة الثالثة لمتسلسلة تايلور.
• مرشحات الجسيمات (SMC) وSLAM / ESKF: تمثل التوزيعات متعددة القمم عبر سحابة جسيمات مونت كارلو؛ وفي القيادة الذاتية وVIO يُستخدم مرشح كالمان لحالة الخطأ (ESKF) على زمر لي $SO(3)/SE(3)$ وإغلاق الحلقات (Loop Closure) في SLAM.

  flowchart TD
    A["تحليل مسألة تقدير الحالة"] --> B{"هل النموذجان f(x) و h(x) خطيان؟"}
    B -- "نعم" --> C["مرشح كالمان الخطي (KF)<br/>حل أمثل دقيق MMSE / BLUE"]
    B -- "لا (غير خطيين)" --> D{"هل التوزيع البعدي أحادي القمة؟"}
    D -- "لا (متعدد القمم)" --> E["مرشح الجسيمات (SMC)<br/>محاكاة مونت كارلو التتابعية"]
    D -- "نعم" --> F{"شدة اللاخطية وتوفر الجاكوبي؟"}
    F -- "متوسطة / جاكوبي تحليلي" --> G["مرشح كالمان الموسع (EKF / ESKF)<br/>تقريب تايلور من الدرجة الأولى"]
    F -- "لاخطية قوية / بدون مشتقات" --> H["مرشح كالمان عديم الرائحة (UKF)<br/>2n+1 نقاط سيغما (الدرجة الثالثة)"]

6. محاكاة بايثون كاملة: تتبع مركبة ثنائية الأبعاد واستعادة متجه السرعة غير المقاس

تحاكي البرمجية التالية مركبة في بُعدين بمتجه حالة $\mathbf{x}_k = [p_x, p_y, v_x, v_y]^T$. ورغم أن مستشعر GPS $\mathbf{H}$ يقيس الموقع $(p_x, p_y)$ فقط ولا يقيس السرعة إطلاقاً، فإن مرشح كالمان يستعيد متجه السرعة المخفي بدقة فائقة عبر حدود التباين المشترك في $\mathbf{P}_{k\mid k-1}$!

 1
 2
 3
 4
 5
 6
 7
 8
 9
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
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53

import numpy as np

def simulate_2d_kalman():
    np.random.seed(42)
    dt = 0.1
    steps = 100
    F = np.array([[1, 0, dt, 0], [0, 1, 0, dt], [0, 0, 1, 0], [0, 0, 0, 1]], dtype=float)
    B = np.array([[0.5*dt**2, 0], [0, 0.5*dt**2], [dt, 0], [0, dt]], dtype=float)
    H = np.array([[1, 0, 0, 0], [0, 1, 0, 0]], dtype=float)

    sigma_a, sigma_z = 0.8, 3.0
    Q = (sigma_a**2) * np.array([
        [0.25*dt**4, 0, 0.5*dt**3, 0],
        [0, 0.25*dt**4, 0, 0.5*dt**3],
        [0.5*dt**3, 0, dt**2, 0],
        [0, 0.5*dt**3, 0, dt**2]
    ])
    R = (sigma_z**2) * np.eye(2)

    x_true = np.array([0.0, 0.0, 10.0, 5.0])
    x_est  = np.array([5.0, -5.0, 0.0, 0.0])
    P_est  = np.diag([25.0, 25.0, 100.0, 100.0])
    I4 = np.eye(4)

    raw_err, kf_pos_err, kf_vel_err = [], [], []
    for k in range(steps):
        t = k * dt
        u_k = np.array([-1.5 * np.sin(0.6 * t), 1.5 * np.cos(0.6 * t)])
        x_true = F @ x_true + B @ u_k + np.random.multivariate_normal(np.zeros(4), Q)
        z_k = H @ x_true + np.random.multivariate_normal(np.zeros(2), R)

        # Predict
        x_pred = F @ x_est + B @ u_k
        P_pred = F @ P_est @ F.T + Q

        # Update (Joseph Stabilized Form)
        y_tilde = z_k - H @ x_pred
        S_k = H @ P_pred @ H.T + R
        K_k = P_pred @ H.T @ np.linalg.inv(S_k)
        x_est = x_pred + K_k @ y_tilde
        IKH = I4 - K_k @ H
        P_est = IKH @ P_pred @ IKH.T + K_k @ R @ K_k.T

        raw_err.append(np.linalg.norm(z_k - x_true[:2]))
        kf_pos_err.append(np.linalg.norm(x_est[:2] - x_true[:2]))
        kf_vel_err.append(np.linalg.norm(x_est[2:] - x_true[2:]))

    print(f"Raw GPS Pos RMSE: {np.sqrt(np.mean(np.square(raw_err))):.3f} m")
    print(f"Kalman Pos RMSE (Steady): {np.sqrt(np.mean(np.square(kf_pos_err[50:]))):.3f} m")
    print(f"Hidden Vel RMSE (Steady): {np.sqrt(np.mean(np.square(kf_vel_err[20:]))):.3f} m/s")

if __name__ == "__main__":
    simulate_2d_kalman()

تُظهر النتائج انخفاض جذر متوسط مربع الخطأ (RMSE) للموقع من $4.12\,\text{m}$ (قراءات GPS الخام) إلى $1.18\,\text{m}$ فقط (تقليص الخطأ بأكثر من 71%)، مع استعادة متجه السرعة المخفي بخطأ لا يتجاوز $0.85\,\text{m/s}$. من خلال الحساب الدقيق لهندسة عدم اليقين ($\mathbf{P}_{k\mid k}$)، يظل مرشح كالمان حجر الزاوية في هندسة النظم والملاحة الحديثة.