攝動(dòng)開普勒問題形式化建模與驗(yàn)證

王國(guó)輝，許京然，劉永梅，施智平，關(guān) 永

1(首都師范大學(xué) 信息工程學(xué)院,電子系統(tǒng)可靠性技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室，北京 100048)2(首都師范大學(xué) 北京成像理論與技術(shù)高精尖創(chuàng)新中心，北京 100048)3(北京城市學(xué)院 信息學(xué)部，北京 101399)4(首都師范大學(xué) 電子系統(tǒng)可靠性與數(shù)理交叉學(xué)科國(guó)家國(guó)際科技合作示范型基地，北京 100048)

1 引 言

人造地球衛(wèi)星在地球引力的作用下沿閉合軌道繞地球做周期性運(yùn)行.當(dāng)?shù)厍蚴且粋€(gè)質(zhì)量分布均勻的標(biāo)準(zhǔn)球體時(shí)，衛(wèi)星的軌道應(yīng)該為標(biāo)準(zhǔn)的橢圓形，并服從開普勒三大定律.然而，地球表面高低起伏，地球內(nèi)部結(jié)構(gòu)復(fù)雜，使得地球引力在空間分布不均.人造地球衛(wèi)星在地球引力場(chǎng)的作用下，其運(yùn)行軌道發(fā)生較為復(fù)雜的變換，產(chǎn)生攝動(dòng)問題[1].因此，攝動(dòng)開普勒問題在衛(wèi)星軌道分析研究過程中得到廣泛應(yīng)用.該領(lǐng)域?qū)Π踩耘c可靠性要求非常嚴(yán)格，微小的錯(cuò)誤將導(dǎo)致災(zāi)難性后果.

早在1920年Levi-Civita就提出二維諧振子描述平面開普勒問題[2].1964年Kustaanheimo和Stiefel提出四維諧振子描述空間開普勒問題[3].這兩種類型將開普勒問題中的微分方程轉(zhuǎn)化為線性微分方程且都是可行的.1995年Vivarelli和Vrbik提出了利用四元數(shù)解決該問題的可行性[4].直到2006年Waldvogel利用四元數(shù)合理的解決了該問題并給出了推導(dǎo)方法與過程.從二維復(fù)數(shù)到三維四元數(shù)的代數(shù)轉(zhuǎn)換，是一次突破性進(jìn)展[5].但是上述方法在解決衛(wèi)星攝動(dòng)開普勒問題的建模與分析過程中涉及到矢量代數(shù)、旋量代數(shù)、復(fù)數(shù)、四元數(shù)等多種不同的代數(shù)系統(tǒng)，在各個(gè)代數(shù)系統(tǒng)相互轉(zhuǎn)換過程中極易引入錯(cuò)誤.2008年哈爾濱工業(yè)大學(xué)方茹等[6]提出了利用幾何代數(shù)解決攝動(dòng)開普勒問題.該方法利用幾何代數(shù)將攝動(dòng)開普勒問題的建模與分析統(tǒng)一到相同代數(shù)結(jié)構(gòu)中，同時(shí)把攝動(dòng)開普勒方程轉(zhuǎn)換成旋量方程來求解，有效避免奇點(diǎn)的同時(shí)彌補(bǔ)傳統(tǒng)分析方法的不足.

由于上述建模與分析方法主要依靠紙筆演算、數(shù)值計(jì)算和計(jì)算機(jī)代數(shù)系統(tǒng).紙筆演算的方法耗時(shí)耗力，容易引入人為錯(cuò)誤；計(jì)算機(jī)數(shù)值計(jì)算方法由于計(jì)算機(jī)無法精確表示實(shí)數(shù)，數(shù)值計(jì)算不能給出精確的結(jié)果；計(jì)算機(jī)代數(shù)系統(tǒng)在邊界條件、奇異表達(dá)簡(jiǎn)化方面的處理具有缺陷，此外龐大的符號(hào)計(jì)算程序不能排除程序漏洞的存在[7].相對(duì)于傳統(tǒng)的建模與分析方法，基于定理證明的形式化方法采用嚴(yán)格的高階邏輯描述數(shù)學(xué)問題的屬性和規(guī)范，以公認(rèn)的邏輯公理和推理規(guī)則為基礎(chǔ)構(gòu)建形式化建模和驗(yàn)證過程，在安全攸關(guān)的系統(tǒng)設(shè)計(jì)中具有很大優(yōu)勢(shì).

中國(guó)科學(xué)院李洪波[8]在2002年首次提出應(yīng)用Clifford代數(shù)自動(dòng)證明共形幾何模型，創(chuàng)立了高效的幾何計(jì)算和推理新方法.2013年Harrison等曾使用HOL Light對(duì)幾何代數(shù)的基本內(nèi)容，主要包括多重矢量和基矢量做出了形式化定義，完成了幾何代數(shù)核心運(yùn)算和部分性質(zhì)的證明[9].2016年馬莎等[10]完成了共形幾何代數(shù)的高階邏輯形式化驗(yàn)證工作，為幾何代數(shù)定理證明做了鋪墊.2018年李黎明[11]在HOL Light定理證明器中構(gòu)建了幾何代數(shù)的高階邏輯定理庫(kù)，為攝動(dòng)開普勒問題的形式化建模與分析提供了必要條件.本文通過定理證明的方法對(duì)基于幾何代數(shù)的攝動(dòng)開普勒問題進(jìn)行形式化建模與分析，并給出嚴(yán)格的高階邏輯推理，從而最大程度確保數(shù)學(xué)模型的正確性和分析方法的可靠性.

2 幾何代數(shù)基礎(chǔ)及其形式化

幾何代數(shù)是由Grassmann 代數(shù)和Clifford 代數(shù)發(fā)展而來的，能夠統(tǒng)一多種代數(shù)結(jié)構(gòu)到同一框架下，從而避免多種代數(shù)系統(tǒng)之間相互轉(zhuǎn)換的問題.利用幾何代數(shù)中幾何體、幾何關(guān)系和幾何變換描述不依賴坐標(biāo)的優(yōu)點(diǎn)，可以將描述衛(wèi)星軌道運(yùn)動(dòng)的攝動(dòng)開普勒方程轉(zhuǎn)化為線性、正則旋量方程，即KS方程[12].

表1 幾何代數(shù)基本運(yùn)算在定理證明庫(kù)中的表示Table 1 Representation of the basic operations of geometric algebra in the theorem proving library

2.1 幾何代數(shù)基本運(yùn)算及其性質(zhì)的形式化

幾何代數(shù)建立在多重向量空間上，其基本元素為多重向量.幾何代數(shù)運(yùn)算主要包括內(nèi)積、外積和幾何積.在使用過程中，內(nèi)積還分為一般內(nèi)積、左縮積、右縮積和標(biāo)量積四種.在HOL Light定理證明庫(kù)中，幾何代數(shù)基本運(yùn)算如表1所示.幾何代數(shù)基本運(yùn)算不依賴研究對(duì)象的坐標(biāo)系統(tǒng)，只與坐標(biāo)基之間的相對(duì)位置有關(guān).其運(yùn)算規(guī)則符合雙線性、結(jié)合律和分配律等基本性質(zhì)，有效的統(tǒng)一了標(biāo)量運(yùn)算、矢量運(yùn)算、維度運(yùn)算和幾何運(yùn)算.

本文以一般內(nèi)積為例簡(jiǎn)要介紹其定義及性質(zhì)的高階邏輯形式化描述.其他運(yùn)算的定義與性質(zhì)可參見幾何代數(shù)的高階邏輯定理證明庫(kù)1.

定義1.幾何代數(shù)一般內(nèi)積

|- x innerga y：real^(P，Q，R)geomalg =

Productga(s t.if s={} / t={} /

～(s SUBSET t)/～(t SUBSET s) then & 0

else --(& 1)pow CARD {i，j | i IN 1..(pdimindex(：P)

+pdimindex(：Q)+pdimindex(：R))/

j IN 1..(pdimindex(：P)+pdimindex(：Q)+

pdimindex(：R))/ i IN s / j IN t / i > j} *

--(& 1)pow CARD((pdimindex(：P)+

1..pdimindex(：P)+pdimindex(：Q))

INTER s INTER t)* & 0 pow

CARD((pdimindex(：P)+pdimindex(：Q)+

1..pdimindex(：P)+pdimindex(：Q)+

pdimindex(：R))INTER s INTER t))x y

其中real^(P，Q，R)geomalg為幾何代數(shù)基本元素多重向量的類型.

性質(zhì)1.一般內(nèi)積雙線性性質(zhì)

|- bilinear f =(!x.linear(y.f x y))/

(!y.linear(y.f x y))

其中l(wèi)inear為線性性質(zhì)描述.

|-linear(f：real^M->real^N)=

(!x y.f(x+y)= f(x)+f(y))/〗(!c x.f(c*x)= c * f(x))

性質(zhì)2.一般內(nèi)積結(jié)合律性質(zhì)

|- !x y z.op x(op y z)= op(op x y)z

其中op為二元操作符.

2.2 位置矢量與旋量

由幾何代數(shù)基礎(chǔ)知識(shí)可知，歐式空間內(nèi)的旋轉(zhuǎn)和伸縮變換滿足如公式(1)所示的歐式空間內(nèi)矢量與幾何空間內(nèi)的旋量或四元數(shù)之間的關(guān)系.

(1)

2.3 多重向量微分形式化

基于幾何代數(shù)理論的衛(wèi)星攝動(dòng)KS方程涉及衛(wèi)星與地球的位置、速度和加速度的關(guān)系.因此需要構(gòu)建多重向量微分及其基本定理的形式化，為開展攝動(dòng)開普勒問題的形式化工作做準(zhǔn)備.Maggesi等[13]在HOL Light系統(tǒng)中完成了向量函數(shù)的微分，本文以此為基礎(chǔ)，實(shí)現(xiàn)多重向量微分形式化.

定義2.多重向量函數(shù)f與微分函數(shù)f′之間的關(guān)系

|-(f has_ multivector _derivative f′)net <=>

(f has_multivector_derivative(x.drop(x)% f′))net

定義3.多重向量函數(shù)f在某點(diǎn)進(jìn)行向量微分的函數(shù)值

|- multivector _derivative(f：real^1->(P，Q，R)geomalg)net=@f′.(f has_multivector _derivative f′)net

|- !x y：real^1->real^(′3，′0，′0)geomalg x′ y′ t.

(x has_ multivector _derivative x′)(at t)/

(y has_ multivector _derivative y′)(at t)

==>(( .(x t inner y t))

has_multivector _derivative

((x(t)inner y′)+x′ inner y(t)))(at t)

|-!x y：real^1->real^(′3，′0，′0)geomalg x′ y′ t.

(x has_ multivector _derivative x′)(at t)/〗(yhas_ multivector _derivative y′)(at t)

==>(( .(x t outer y t))

has_ multivector _derivative

((x(t)outer y′)+x′ outer y(t)))(at t)

|-!x y：real^1->real^(′3，′0，′0)geomalg x′ y′ t.

(x has_ multivector _derivative x′)(at t)/〗(y has_ multivector _derivative y′)(at t)

==>(( .(x t * y t))

has_ multivector _derivative

((x(t)* y′)+x′ * y(t)))(at t)

上述三個(gè)多重向量函數(shù)運(yùn)算微分的性質(zhì)形式證明，為衛(wèi)星攝動(dòng)開普勒問題形式化建模奠定了基礎(chǔ).

3 衛(wèi)星攝動(dòng)開普勒問題形式化建模

3.1 慣性坐標(biāo)系中攝動(dòng)開普勒方程形式化

人造地球衛(wèi)星在三維歐式空間的慣性坐標(biāo)系中，衛(wèi)星與地球之間有微小且不可忽略的攝動(dòng)力，在微小的攝動(dòng)力f的作用下，系統(tǒng)的攝動(dòng)開普勒方程可以表示為方程(2).

(2)

定義4.慣性坐標(biāo)系中攝動(dòng)開普勒方程形式化描述

|-multivector_derivative((multivector_derivative

(rotation_t r)(at t))(at t)=

((-k：real^1)*(r_position q t)/

norm(r_position q)pow & 3)+

(f：real^1->real^(′3，′0，′0)trip_fin_sum)

其中變量rotation_t r表示的是天體的矢量位置，它是一個(gè)與時(shí)間t有關(guān)的矢量函數(shù)；f是一個(gè)與時(shí)間t有關(guān)的多重向量函數(shù)；multivector_derivative f(at t)表示函數(shù)f對(duì)時(shí)間t的微分.k是一個(gè)實(shí)數(shù)，即公式(1)中的變量μ.

3.2 基于幾何代數(shù)的攝動(dòng)開普勒方程形式化

利用幾何代數(shù)相關(guān)理論，可以將攝動(dòng)開普勒方程轉(zhuǎn)換為運(yùn)動(dòng)的旋量方程，記為KS方程.該方程形式簡(jiǎn)單，計(jì)算方便，同時(shí)能夠有效消除引力位奇異點(diǎn)1/r，如方程(3)所示.

(3)

定義5.基于幾何代數(shù)的攝動(dòng)開普勒方程形式化描述

|-& 2 % multivector_derivative

((multivector_derivative q)(at s))(at s)-(E_DEF q t k)*(q t)=q t*r_position u t *

(f：real^1->real^(′3，′0，′0)trip_fin_sum)

其中開普勒能量E定義6所示.

定義6.開普勒能量

|-E_DEF(q：real^1->real^(′3，′0，′0)geomalg)(t：real^1)(k：real)

=(& 2 %

((conjugation(multivector_derivative q(at t))*

(multivector_derivative q(at t)))MYMMYM {})- k)

/norm(r_position q t)

4 衛(wèi)星攝動(dòng)開普勒問題形式化驗(yàn)證

本小節(jié)將利用定理證明庫(kù)中已有的和補(bǔ)充的定義、定理對(duì)基于幾何代數(shù)的攝動(dòng)開普勒方程進(jìn)行形式化推導(dǎo).根據(jù)能量守恒原則，證明上述兩種攝動(dòng)開普勒問題建模分析方法的等價(jià)性.從而驗(yàn)證了基于幾何代數(shù)的攝動(dòng)開普勒方程的正確性與完備性.

4.1 地球與衛(wèi)星相對(duì)位置矢量形式化定義

在歐式慣性坐標(biāo)系中，地球與人造地球衛(wèi)星中心相對(duì)位置矢量關(guān)系滿足公式(4)：

(4)

其中u為采用幾何代數(shù)系統(tǒng)中多重向量的四元數(shù)表示，u+為u的共軛.

定義7.地球與衛(wèi)星相對(duì)位置矢量形式化

|-r_position(q：quat)=

(geomalg_300_quat q)* mbasis{1} *

conjugation(geomalg_300_quat q)

其中，quat為四元數(shù)類型.conjugation(u)表示u的共軛.e1是一個(gè)矢量基，用mbasis{1}表示.geomalg_300_quat表示四元數(shù)表示與多重向量表示的轉(zhuǎn)換，其形式化定義如定義7所示.

定義8.四元數(shù)與多重向量轉(zhuǎn)換形式化定義

|-geomalg_300_quat(q：quat)=

(Re q)% mbasis{} +(Im1 q)% mbasis{2，3} +

(Im2 q)% mbasis{1，2} +(Im3 q)%

(--(mbasis{1，3}：real^(′3，′0，′0)geomalg))

4.2 地球與衛(wèi)星相對(duì)位置形式化相關(guān)定理

由地球與人造地球衛(wèi)星中心相對(duì)位置矢量關(guān)系公式(4)可得關(guān)系式(5)：

(5)

定理1.地球與衛(wèi)星相對(duì)位置標(biāo)量化

|-!q：quat.norm(r_position q)=

(norm(geomalg_300_quat q))pow 2

其中，norm表示矢量取模.

為了證明定理1的成立，需要引入兩個(gè)引理.即引理1位置矢量與多重矢量共軛的關(guān)系和引理2多重矢量與其共軛乘積的交互性.首先，用重寫策略(REWRITE_TAC)將定理1中的r_position和geomalg_300_quat定義展開.其次，運(yùn)用化簡(jiǎn)策略(SIMP_TAC)并結(jié)合HOL Light向量庫(kù)中的相關(guān)定理對(duì)目標(biāo)進(jìn)行化簡(jiǎn).最后，將定理目標(biāo)推導(dǎo)為實(shí)數(shù)相等，并應(yīng)用實(shí)數(shù)推導(dǎo)自動(dòng)策略(REAL_ARITH_TAC)完成證明.

引理1.位置矢量與多重矢量共軛關(guān)系

|-!q：quat.norm(r_position q)% mbasis{}=

geomalg_300_quat q *

conjugation(geomalg_300_quat q)

引理2.多重矢量與其共軛積的交互性

|-!q：quat.conjugation(geomalg_300_quat q)*

geomalg_300_quat q=

geomalg_300_quat q *

conjugation(geomalg_300_quat q)

4.3 基于幾何代數(shù)攝動(dòng)開普勒方程形式化推導(dǎo)

(6)

在形式化推導(dǎo)驗(yàn)證過程中，直接證明公式(6)成立相對(duì)復(fù)雜，所以先由多重向量的共軛性質(zhì)(u+)+=u去證明等式(7).

(7)

|-!q′ q：quat.

((geomalg_300_quat q′ * mbasis{1} *

geomalg_300_quat(cnj q))$${1，2，3}=& 0 )

==>conjugation(geomalg_300_quat q′*mbasis{1} *

geomalg_300_quat(cnj q))=

geomalg_300_quat q′*mbasis{1}*

geomalg_300_quat(cnj q)

|-!q′ q：quat.((geomalg_300_quat q′*mbasis{1}*

geomalg_300_quat(cnj q))$${1，2，3}=& 0 )

==>geomalg_300_quat q′*mbasis{1}*

geomalg_300_quat(cnj q)=

geomalg_300_quat q*mbasis{1}*

geomalg_300_quat(cnj q′)

根據(jù)引理3和引理4可以把公式(6)化簡(jiǎn)，如公式(8)所示.

(8)

其高階邏輯形式化表示如定理2所示.

|-!q′：quat q：real^1->quat t：real^1.

((geomalg_300_quat q′*mbasis{1}*

geomalg_300_quat(cnj(q t)))$${1，2，3}=& 0 /〗((q)has_multivector_derivative q′)(at t))

==>(( .(r_position(q t)))

has_multivector_derivative

(& 2 %(geomalg_300_quat q′*mbasis{1}*

geomalg_300_quat(cnj(q t)))))(at t)

公式(8)左右兩邊同乘e1u，再對(duì)時(shí)間t求導(dǎo)，可得公式(9).

(9)

其形式化描述為定理3.

|-!(q：real^1->real^(′3，′0，′0)geomalg)

(s：real^1->real^(′3，′0，′0)trip_fin_sum)(t：real^1).

& 2 % multivector_derivative((r_position q t)*

(multivectorr_derivative q(at t)))(at t)=

& 2 / norm(r_position q t)*

multivector_derivative(

(multivector_derivative q)(at s))(at s)

此時(shí)可以得到關(guān)于u對(duì)s的二次微分，如公式(10)所示.

(10)

(11)

為了證明公式(11)成立，在邏輯推導(dǎo)過程中仍需引入dt=rds表示t與s的關(guān)系，用引理5表示其形式化描述.

|-!(q：real^1->real^(′3，′0，′0)geomalg)(t：real^1)

(s：real^1->real^(′3，′0，′0)trip_fin_sum).

multivector_derivative s(at t)=

& 1 / norm(r_position q t)

通過上述推導(dǎo)，可以獲得地球與人造地球衛(wèi)星攝動(dòng)開普勒問題的幾何代數(shù)高階邏輯模型，并且形式化證明了幾何代數(shù)方法構(gòu)建衛(wèi)星攝動(dòng)開普勒問題數(shù)學(xué)模型的正確性.

4.4 攝動(dòng)開普勒問題幾何代數(shù)模型與慣性坐標(biāo)模型等價(jià)性證明

根據(jù)基于幾何代數(shù)的攝動(dòng)開普勒方程公式(2)可知，公式(11)滿足攝動(dòng)開普勒問題的旋量方程，根據(jù)能量守恒原則，可證如下公式(12)成立，

(12)

公式(12)可形式化為定理4，即開普勒能量守恒性質(zhì).

定理4.開普勒能量性質(zhì)

|-!(q：real^1->real^(′3，′0，′0)geomalg)(t：real^1)

(k：real).E_DEF q t k =

(norm(multivector_derivative(r_position q)

(at t))pow 2 / & 2)-k/norm(r_position q t)

為了證明定理4的成立，需要證明開普勒能量引理，即引理6.

引理6.開普勒能量引理

|-!(q：real^1->real^(′3，′0，′0)geomalg)(t：real^1).

norm(r_position q t)*

(norm(multivector_derivative(r_position q)

(at t))pow 2)= & 4*norm

(multivector_derivative q(at t))pow 2

通過引理6的證明我們可以得到開普勒能量在幾何代數(shù)下的形式化描述.定理4以開普勒能量守恒為形式化驗(yàn)證的最終目標(biāo).其證明過程主要借助上文提到的定義、引理和定理，運(yùn)行重寫、化簡(jiǎn)和自動(dòng)求解等策略組合實(shí)現(xiàn)邏輯推導(dǎo).從而證明攝動(dòng)開普勒問題幾何代數(shù)模型與慣性坐標(biāo)模型具有等價(jià)性.

本小節(jié)的重點(diǎn)是將攝動(dòng)開普勒旋量方程的形式化建模與驗(yàn)證.經(jīng)過過嚴(yán)密的高階邏輯形式推導(dǎo)，最大程度確保基于幾何代數(shù)的攝動(dòng)開普勒問題數(shù)學(xué)模型的正確性和分析方法的可靠性.

5 結(jié) 論

幾何代數(shù)系統(tǒng)能夠?qū)⒍喾N代數(shù)系統(tǒng)統(tǒng)一在同一個(gè)代數(shù)系統(tǒng)中，避免了代數(shù)系統(tǒng)之間轉(zhuǎn)換問題，有效的提高了衛(wèi)星攝動(dòng)開普勒問題分析的可靠性.本文以幾何代數(shù)為數(shù)學(xué)基礎(chǔ)，在高階邏輯定理證明器HOL Light中建立攝動(dòng)開普勒方程邏輯模型，對(duì)攝動(dòng)開普勒問題進(jìn)行了形式化分析與驗(yàn)證.由于解決攝動(dòng)開普勒問題的需要，本文首先補(bǔ)充了幾何代數(shù)基本運(yùn)算微分性質(zhì)定理；其次完成從歐式空間到幾何代數(shù)空間的轉(zhuǎn)換關(guān)系進(jìn)行形式化；再次形式化證明多重向量的線性性質(zhì)、共軛性質(zhì)以及共軛的線性性質(zhì)；最后采用高階邏輯對(duì)基于幾何代數(shù)的衛(wèi)星攝動(dòng)開普勒問題的數(shù)學(xué)模型進(jìn)行形式化建模與推導(dǎo)，從而最大程度確保數(shù)學(xué)模型的正確性和分析方法的可靠性.

此外，攝動(dòng)開普勒問題不僅應(yīng)用于研究衛(wèi)星軌道和姿態(tài)運(yùn)動(dòng)，也是研究?jī)蓚€(gè)帶電粒子的運(yùn)動(dòng)方程基礎(chǔ)，在測(cè)試物理理論和測(cè)量自然常數(shù)的模型系統(tǒng)中發(fā)揮了重要作用.因此下一步工作將圍繞基于幾何代數(shù)理論的微觀世界粒子運(yùn)動(dòng)方程的形式化建模與驗(yàn)證展開.