彈載圖像探測(cè)器廣域協(xié)同探測(cè)算法

婁文忠，蘇子龍，汪金奎，劉偉桐，趙飛

(北京理工大學(xué) 爆炸科學(xué)與技術(shù)國家重點(diǎn)實(shí)驗(yàn)室，北京 100081)

0 引言

作戰(zhàn)毀傷效果評(píng)估的定義是，在對(duì)既定目標(biāo)進(jìn)行軍事打擊(包括致命和非致命)后，對(duì)目標(biāo)進(jìn)行及時(shí)和準(zhǔn)確的毀傷估計(jì)[1]。基于圖像的毀傷評(píng)估是戰(zhàn)場(chǎng)毀傷評(píng)估的主要方法，通過對(duì)打擊前后目標(biāo)圖像采集，分析目標(biāo)毀傷前后圖像變化特征，可獲得直觀準(zhǔn)確的毀傷效果信息。

由于星載圖像探測(cè)方式使用權(quán)限高、軌道周期長以及機(jī)載圖像探測(cè)方式飛行周期長等缺點(diǎn)，彈載圖像探測(cè)方式在區(qū)域作戰(zhàn)的體系對(duì)抗中具有巨大潛力。基于彈載的目標(biāo)毀傷圖像捕獲技術(shù)是作戰(zhàn)毀傷信息獲取的典型方法[2]。李從利等[3]研究了偵察彈圖像質(zhì)量評(píng)價(jià)方法，針對(duì)彈載偵察圖像由于受成像平臺(tái)和條件限制存在多類型混雜失真,提出了評(píng)價(jià)方法，為圖像后續(xù)處理和成像系統(tǒng)性能優(yōu)化提供了量化依據(jù)和指標(biāo)參考；王正平等[4]針對(duì)炮射偵察彈的工作特性,提出了彈道設(shè)計(jì)方案，建立了炮射段、分離段、減速段彈道模型,推導(dǎo)了開傘時(shí)間與分離時(shí)間、降落傘面積與彈道末端約束的關(guān)系；李卉等[5]研究了電視偵察彈與艦炮武器系統(tǒng)的匹配問題，以武器系統(tǒng)綜合作戰(zhàn)效能為目標(biāo)函數(shù)進(jìn)行系統(tǒng)指標(biāo)優(yōu)化，對(duì)艦炮武器系統(tǒng)作戰(zhàn)使用電視偵察系統(tǒng)提供參考。

彈載毀傷圖像感知中存在多組約束條件，首先是視場(chǎng)角與評(píng)估圖像清晰度的約束，視場(chǎng)角越小，獲取的圖像清晰度越高，但圖像靶面上可獲取的毀傷區(qū)域面積越小；探測(cè)到目標(biāo)的概率可通過仿真及試驗(yàn)制定，是系統(tǒng)的輸入?yún)?shù);彈載圖像毀傷評(píng)估系統(tǒng)通常不具備巡飛能力，獲取毀傷圖像信息的窗口窄，系統(tǒng)難以對(duì)毀傷過程的全局圖像信息感知。多智能體系通過無線信息交聯(lián)，實(shí)現(xiàn)廣域視野協(xié)同探測(cè)，是現(xiàn)代戰(zhàn)場(chǎng)中越來越重要的信息獲取模式[6]，彈載圖像探測(cè)器是通過彈藥拋撒出多個(gè)子彈藥，分時(shí)分區(qū)域獲取毀傷圖像信息，通過圖像的時(shí)間拼接與空間拼接，獲取大區(qū)域范圍的高質(zhì)量毀傷圖像信息，可提供戰(zhàn)場(chǎng)毀傷圖像信息獲取的有效解決方案。

1 多約束圖像探測(cè)器建模

戰(zhàn)場(chǎng)目標(biāo)通常可以分為3種類型，已知固定目標(biāo)、已知移動(dòng)目標(biāo)與臨時(shí)目標(biāo)。其中已知固定目標(biāo)為作戰(zhàn)前通過其他途徑獲取的固定目標(biāo)，其坐標(biāo)點(diǎn)已知；已知移動(dòng)目標(biāo)雖然可在戰(zhàn)前獲取，但目標(biāo)可能在探測(cè)過程中移動(dòng)位置；臨時(shí)目標(biāo)指在作戰(zhàn)過程中臨時(shí)發(fā)現(xiàn)的高價(jià)值目標(biāo)，其坐標(biāo)在目標(biāo)發(fā)現(xiàn)后才能獲取。本文主要討論已知固定目標(biāo)類型，并假定圖像子探測(cè)器數(shù)量與目標(biāo)數(shù)量一致，預(yù)設(shè)條件如表1所示。

表1 戰(zhàn)場(chǎng)與圖像探測(cè)器預(yù)設(shè)條件

為準(zhǔn)確描述協(xié)同控制系統(tǒng)，需要定義以下預(yù)設(shè)條件：

預(yù)設(shè)條件1表示作戰(zhàn)區(qū)域內(nèi)的多個(gè)目標(biāo)存在重要性差異；為簡化系統(tǒng)設(shè)計(jì)和適應(yīng)彈載環(huán)境，預(yù)設(shè)條件2、條件3、條件4將圖像探測(cè)器定義為可完成簡易修正，并可快速切換視場(chǎng)的系統(tǒng)模型。圖像探測(cè)器的控制系統(tǒng)由于系統(tǒng)復(fù)雜度和設(shè)計(jì)成本，通常難以采用舵機(jī)等三維修正機(jī)構(gòu)，本文假設(shè)圖像探測(cè)器具備一維修正能力，同時(shí)可通過云臺(tái)控制視線軸方向。修正力沿圖像探測(cè)器軸線方向，因此探測(cè)器均在彈道平面內(nèi)運(yùn)動(dòng)，二維動(dòng)力學(xué)方程[7]表示為

(1)

式中：m為圖像探測(cè)器質(zhì)量；v表示圖像探測(cè)器飛行速度；Fx和Fy分別為探測(cè)器的阻力和升力，x、y分別為探測(cè)器在彈道平面內(nèi)的水平位置和垂直位置；g為當(dāng)?shù)刂亓铀俣龋沪葹樘綔y(cè)器彈道傾角；Jz為探測(cè)器對(duì)載體坐標(biāo)系z(mì)軸的轉(zhuǎn)動(dòng)慣量；ωz為繞z軸的轉(zhuǎn)動(dòng)角速度；Mz為俯仰力矩；?為俯仰角；α為攻角。

在圖像探測(cè)器上采用增阻式一維彈道修正模型，其攻角為0°和升力為0 N，因此探測(cè)器動(dòng)力學(xué)方程可簡化為(2)式。

(2)

式中：C為探測(cè)器常態(tài)飛行下阻力系數(shù)，Ci為預(yù)設(shè)定增阻器開啟后的C；ρ為空氣密度；S為氣動(dòng)面積；k為放大倍數(shù)；ti為預(yù)設(shè)定增阻器開啟時(shí)間。通過增阻可改變探測(cè)器阻力系數(shù)，當(dāng)增阻器打開時(shí)，可認(rèn)為其氣動(dòng)面積保持不變，阻力系數(shù)增大為原來的k倍[8]。

為保證目標(biāo)不失一般性，面目標(biāo)可定義為多個(gè)點(diǎn)目標(biāo)構(gòu)成的集合，表示為T={T1,T2,…,Tn}，其中：n代表探測(cè)區(qū)域中目標(biāo)個(gè)數(shù)；Tj={Xj,sj,Wj},j∈{1,2,…,n}，j代表目標(biāo)編號(hào)，協(xié)同探測(cè)模型如圖1所示；Xj表示點(diǎn)目標(biāo)的坐標(biāo)點(diǎn)；sj表示點(diǎn)目標(biāo)的有效信息面積，在此區(qū)域內(nèi)的圖像信息為毀傷評(píng)估的有效信息；此外還需對(duì)目標(biāo)區(qū)域定義價(jià)值系數(shù)Wj.

圖1 協(xié)同探測(cè)模型

2 圖像探測(cè)器探測(cè)能力評(píng)價(jià)準(zhǔn)則

通常情況下，有效圖像信息量不會(huì)隨時(shí)間線性增長，且隨時(shí)間增長獲取新信息的速度會(huì)降低，在文獻(xiàn)[6]中采用指數(shù)函數(shù)定義了在圖像覆蓋區(qū)域內(nèi)，系統(tǒng)獲取的有效信息的收益函數(shù)R(t)及收益增量函數(shù)ΔR(t)為

R(t)=1-(1-R0)e-λt，

(3)

ΔR(t)=λ(1-R0)e-λtΔt，

(4)

式中：R0表示目標(biāo)區(qū)域內(nèi)具備的初始信息量；Δt為采樣周期；λ為傳感器性能指數(shù)，即當(dāng)時(shí)間趨于無窮大時(shí)，可對(duì)該區(qū)域?qū)崿F(xiàn)完全探測(cè)，圖像探測(cè)器探測(cè)面積模型如圖2所示。

圖2 圖像探測(cè)器探測(cè)面積模型

圖2中：d代表探測(cè)器到目標(biāo)之間的距離；θd為探測(cè)器探測(cè)角度(為常量)；αd為圖像探測(cè)器中線與地面夾角。為保證圖像清晰度，通常αd遠(yuǎn)大于θd/2，將探測(cè)區(qū)域投影到地面，有效區(qū)域等效為圓形，其直徑為

(5)

因此探測(cè)面積為

(6)

針對(duì)圖像傳感器而言，傳感器性能指數(shù)λ通過圖像傳感器探測(cè)面積和分辨率共同決定。當(dāng)探測(cè)面積s超過有效信息區(qū)域sj時(shí)，探測(cè)面積越大，則有效區(qū)域所占的比重越低，λ越小；當(dāng)探測(cè)面積未超過有效信息區(qū)域時(shí)，探測(cè)面積越大，所占有效區(qū)域比重越高，λ越大。

(7)

協(xié)同探測(cè)系統(tǒng)對(duì)目標(biāo)有效信息的信息熵增益作為系統(tǒng)的評(píng)價(jià)標(biāo)準(zhǔn)，有效信息的信息熵通過有效圖像信息收益和收益信息權(quán)重共同描述[9]。在系統(tǒng)工作的tf時(shí)間區(qū)間內(nèi)獲得的信息熵增益D表示為

(8)

式中：i為探測(cè)器編號(hào)，i=1，2，3，…，f；ΔRij(t)為編號(hào)i的探測(cè)器對(duì)編號(hào)j目標(biāo)的收益增量。彈載圖像探測(cè)器通過合理分配修正時(shí)機(jī)與觀測(cè)目標(biāo)，使得系統(tǒng)工作階段可獲取的信息熵最大，圖像探測(cè)器在切換目標(biāo)后，收益函數(shù)中的初始信息量為切換之前系統(tǒng)對(duì)該目標(biāo)所獲取的信息總量，并假設(shè)目標(biāo)在完成修正并穩(wěn)定后具備目標(biāo)信息獲取能力，如圖3所示，圖像探測(cè)器可針對(duì)不同目標(biāo)進(jìn)行探測(cè)。

圖3 目標(biāo)分配及協(xié)同探測(cè)

3 圖像探測(cè)器探測(cè)能力評(píng)價(jià)準(zhǔn)則

3.1 二分圖模型及標(biāo)準(zhǔn)Kuhn-Munkres算法

通過控制圖像子探測(cè)器修正時(shí)間與觀測(cè)目標(biāo)，實(shí)現(xiàn)信息熵最大化，本質(zhì)上是求解多變量的優(yōu)化問題，本系統(tǒng)中涉及到多個(gè)探測(cè)器與多個(gè)目標(biāo)之間的多邊關(guān)系，可采用代數(shù)圖[10]描述。圖是由一組對(duì)象以及對(duì)象之間的關(guān)系構(gòu)成的圖形，其中每個(gè)對(duì)象可看作圖的一個(gè)節(jié)點(diǎn)，對(duì)象之間的關(guān)系可看作圖的邊，由節(jié)點(diǎn)和邊構(gòu)成的圖記為G(V,E)，其中：V表示節(jié)點(diǎn)構(gòu)成的集合，V={1,2,…,u}，u為節(jié)點(diǎn)個(gè)數(shù)；E表示節(jié)點(diǎn)g、h的邊構(gòu)成的集合[11]，E={(g,h)∶g,h∈V}。多探測(cè)器觀測(cè)多個(gè)目標(biāo)，可將系統(tǒng)定義為二分圖模型。

定義1若頂點(diǎn)集分成兩個(gè)非空子集Xn和Yn，并且每條邊都有一個(gè)頂點(diǎn)在Xn中，另一個(gè)頂點(diǎn)在Yn中，則稱此圖為二分圖，模型如圖4所示。

圖4 不同觀測(cè)方式構(gòu)成的二分圖

定義2設(shè)具有二分類(Xn，Yn)的賦值二分圖G，其中兩個(gè)頂點(diǎn)分別為Xn={X1,X2,…,Xu}，Yn={Y1,Y2,…,Yu}，任意邊Xg、Yh賦權(quán)ωgh=ω(Xg，Yh)，設(shè)L為G的頂點(diǎn)集V到實(shí)數(shù)集R的映射，若對(duì)任意X∈Xn，Y∈Yn，均有L(X)+L(Y)≥ω(X，Y)，則L為G的可行頂點(diǎn)標(biāo)記；令EL={XY|e=XY∈E(G)，L(X)+L(Y)=ω(e)}，e為頂點(diǎn)X、Y的映射邊，則稱EL為邊集的二分圖G的生成子圖為G的相等子圖[12]。

求得匹配M權(quán)值最大可使用Kuhn-Munkres算法求解，設(shè)G=(V,E)為賦權(quán)二分圖，L是其中一個(gè)初始可行頂點(diǎn)標(biāo)記[13]為

(9)

M是圖G的相等子圖GL的一個(gè)匹配，算法步驟見表2.

表2 Kuhn-Munkres算法

3.2 改進(jìn)Kuhn-Munkres算法

探測(cè)器完成修正后其彈道即可完全預(yù)測(cè)，在時(shí)域空間內(nèi)對(duì)多個(gè)目標(biāo)觀測(cè)收益可完全預(yù)測(cè)，二分圖模型中賦值權(quán)重ωij可等效于i探測(cè)器對(duì)j目標(biāo)的信息收益量，即

ωij(t)=[1-(1-R0j)e-λijt]Wj，

(10)

式中：R0j為目標(biāo)j的初始信息量；λij為探測(cè)器i對(duì)目標(biāo)j的性能指數(shù)。

Kuhn-Munkres算法可解決瞬時(shí)最佳觀測(cè)目標(biāo)分配問題，由于圖像子探測(cè)器在修正后彈道確定，因此在全彈道過程中每時(shí)刻對(duì)目標(biāo)探測(cè)的收益均為常量，可通過Kuhn-Munkres算法對(duì)每個(gè)采樣點(diǎn)的最優(yōu)分配方案計(jì)算并累加，獲得全局最優(yōu)探測(cè)方案，算法運(yùn)算邏輯如圖5所示。

圖5 改進(jìn)Kuhn-Munkres算法模型

在此基礎(chǔ)上，探測(cè)方案為系統(tǒng)修正時(shí)機(jī)的函數(shù)，即

(11)

式中：τ表示所有子探測(cè)器的修正時(shí)機(jī)；f為探測(cè)器個(gè)數(shù)。針對(duì)不同目標(biāo)，(10)式難以求得解析解，因此考慮數(shù)值方法求取最大信息熵增益，將上述問題轉(zhuǎn)化為多元函數(shù)在有界閉區(qū)間[0，tf]內(nèi)的最值問題，其中tf為未修正情況下的圖像子探測(cè)器理論飛行時(shí)間，采用Runge-Kutta求解圖像子探測(cè)器的最優(yōu)修正時(shí)機(jī)[14]。改進(jìn)的Kuhn-Munkres協(xié)同探測(cè)算法步驟如表3所示。

表3 改進(jìn)Kuhn-Munkres算法

3.3 圖像探測(cè)器協(xié)同算法仿真分析

假設(shè)當(dāng)前戰(zhàn)場(chǎng)有3個(gè)已知固定目標(biāo)和3個(gè)圖像探測(cè)器，分別為{x1,x2,x3}以及{y1,y2,y3}，初始數(shù)據(jù)如表4和表5所示。

表4 目標(biāo)位置信息

利用表5中的數(shù)據(jù)帶入(2)式中，利用Runge-Kutta算法可以解算出探測(cè)器的速度函數(shù)v(t)以及姿態(tài)函數(shù)P(t)，如圖6所示。

表5 圖像探測(cè)器初始數(shù)據(jù)

圖6 探測(cè)器運(yùn)動(dòng)參數(shù)

根據(jù)彈道參數(shù)回帶(2)式中，計(jì)算出探測(cè)器位置函數(shù)x(t)和y(t)，如圖7所示。

圖7 探測(cè)器飛行軌跡以及目標(biāo)位置

聯(lián)立表4中的目標(biāo)位置信息、表5中的圖像探測(cè)器初始數(shù)據(jù)、速度函數(shù)v(t)、姿態(tài)函數(shù)P(t)和位置函數(shù)x(t)、y(t)，利用改進(jìn)Kuhn-Munkres算法，實(shí)現(xiàn)傳感器的最優(yōu)目標(biāo)分配，結(jié)果如圖8所示。

圖8 不同探測(cè)器分配方式的目標(biāo)信息增益

圖8圖例中：1-1，前一位數(shù)字代表目標(biāo)編號(hào)1，后一位數(shù)字為圖像探測(cè)器編號(hào)1；其他數(shù)字組合，以此類推。從圖8中可以看出，當(dāng)圖像探測(cè)器越過檢測(cè)目標(biāo)時(shí)，改進(jìn)的Kuhn-Munkres算法可以及時(shí)調(diào)整探測(cè)器的分配方式，目標(biāo)信息熵增益同比增長24%。

由于戰(zhàn)場(chǎng)環(huán)境錯(cuò)綜復(fù)雜，彈載圖像探測(cè)器在雙方體系對(duì)抗中極有可能被摧毀，為了模擬戰(zhàn)場(chǎng)環(huán)境，驗(yàn)證算法可靠性，做如下仿真實(shí)驗(yàn)，設(shè)定圖像探測(cè)器2在1.7 s后被摧毀，結(jié)果如圖9所示。

圖9 探測(cè)器運(yùn)動(dòng)參數(shù)

計(jì)算出探測(cè)器的位置函數(shù)x(t)和y(t)，實(shí)現(xiàn)最優(yōu)分配如圖10、圖11所示。

圖10 探測(cè)器飛行軌跡及目標(biāo)位置

圖11 不同探測(cè)器分配方式的目標(biāo)信息增益

即使彈載圖像探測(cè)器在飛行過程中被摧毀，改進(jìn)的Kuhn-Munkres算法仍能夠有效實(shí)現(xiàn)目標(biāo)分配，實(shí)現(xiàn)目標(biāo)信息增益最大化。

4 結(jié)論

本文從毀傷評(píng)估系統(tǒng)中彈載圖像探測(cè)器對(duì)集群面目標(biāo)的探測(cè)出發(fā)開展研究，針對(duì)圖像探測(cè)器廣域協(xié)同探測(cè)問題，建立了圖像探測(cè)器的動(dòng)力學(xué)模型，將系統(tǒng)可控變量定義為探測(cè)器觀測(cè)目標(biāo)分配方法，用集群圖像探測(cè)器獲取的總體信息熵為評(píng)價(jià)標(biāo)準(zhǔn)。以典型戰(zhàn)場(chǎng)事件為例，得出如下主要結(jié)論：

1)針對(duì)圖像探測(cè)器對(duì)集群面目標(biāo)的廣域協(xié)同探測(cè)問題，綜合Kuhn-Munkres算法和Runge-kutta算法，提出一種改進(jìn)的Kuhn-Munkres算法，獲取最多有效毀傷信息。

2)在圖像探測(cè)器與目標(biāo)數(shù)量一致和不一致的條件下，實(shí)現(xiàn)了彈載圖像探測(cè)器協(xié)同控制，解決了圖像探測(cè)器協(xié)同探測(cè)的最優(yōu)分配問題。

3)改進(jìn)的Kuhn-Munkres算法有效解決集群圖像感知系統(tǒng)對(duì)多目標(biāo)感知中的最優(yōu)控制策略問題，實(shí)現(xiàn)目標(biāo)信息增益最大化。

本文未考慮更加復(fù)雜的戰(zhàn)場(chǎng)環(huán)境因素問題，其最佳控制策略還有待未來進(jìn)一步研究。

參考文獻(xiàn)(References)

[1] 馬春茂,孫衛(wèi)平,李炎,等.武器裝備毀傷評(píng)估研究進(jìn)展[J].火炮發(fā)射與控制學(xué)報(bào),2019,40(4):96-101.

MA C M,SUN W P,LI Y,et al.Research progress on damage assessment of weapons and equipment[J].Journal of Gun Launch & Control, 2019,40(4):96-101.(in Chinese)

[2] 方有培,汪立萍, 趙霜.信息偵察彈的國外發(fā)展及其設(shè)計(jì)構(gòu)想[J].航天電子對(duì)抗, 2004(4):13-15,19.

FANG Y P,WANG L P,ZHAO S.Foreign development and design conception of information reconnaissance projectile[J].Aerospace Electronic Warfare, 2004(4):13-15,19.(in Chinese)

[3] 李從利,薛松,陸文駿,等.彈載偵察圖像質(zhì)量評(píng)價(jià)方法研究[J].兵工學(xué)報(bào),2017,38(1):64-72.

LI C L,XUE S,LU W J,et al.Research on evaluation method of missile-borne reconnaissance image quality[J].Acta Armamentarii, 2017,38(1):64-72.(in Chinese)

[4] 王正平,劉莉,朱勇.炮射偵察彈炮射減速段彈道優(yōu)化設(shè)計(jì)[J].彈道學(xué)報(bào),2012,24(4):22-26.

WANG Z P,LIU L,ZHU Y.Optimal design of artillery-launched reconnaissance projectile’s artillery shooting deceleration phase[J].Journal of Ballistics, 2012,24(4):22-26.(in Chinese)

[5] 李卉,邱志明,王航宇,等.電視偵察彈與艦炮武器系統(tǒng)的優(yōu)化匹配分析[J].兵工學(xué)報(bào),2008,29(12):1462-1466.

LI H,QIU Z M,WANG H Y,et al.Optimal matching analysis of TV reconnaissance projectile and naval gun weapon system[J].Acta Armamentarii, 2008, 29(12): 1462-1466.(in Chinese)

[6] 向庭立，王紅軍,史英春.區(qū)域覆蓋的多機(jī)協(xié)同探測(cè)任務(wù)分配策略[J].空軍工程大學(xué)學(xué)報(bào)(自然科學(xué)版),2019,20(6):33-38.

XIANG T L,WANG H J,SHI Y C.Multi-aircraft cooperative detection task allocation strategy based on area coverage[J].Journal of Air Force Engineering University(Natural Science Edition), 2019,20(6):33-38.(in Chinese)

[7] 錢杏芳，林瑞雄，趙亞男.導(dǎo)彈飛行力學(xué)[M].北京：北京理工大學(xué)出版社，2000.

QIAN X F,LIN R X,ZHAO Y N.Missile flight mechanics[M].Beijing: Beijing Institute of Technology Press, 2000.(in Chinese)

[8] 陳俊.一維彈道修正彈修正機(jī)構(gòu)及控制算法研究[D].南京：南京理工大學(xué), 2013.

CHEN J.Research on correction mechanism and control algorithm of one-dimensional trajectory correction projectile[D].Nanjing: Nanjing University of Science and Technology, 2013.(in Chinese)

[9] PEI Y E, LI L I, XIAO Z, et al.Dynamic spectrum allocation based on Kuhn-Munkres algorithm to guarantee cognitive users’ QoS[J].Journal of Shanghai Normal University(Natural ences), 2013, 42(2):137-142.

[10] 王樹禾.圖論及其算法[M].合肥：中國科學(xué)技術(shù)大學(xué)出版社, 1990.

WANG S H.Graph theory and its algorithm[M].Hefei:University of Science and Technology of China Press, 1990.(in Chinese)

[11] TEMELSO B, MABEY J M, KUBOTA T, et al.ArbAlign: a tool for optimal alignment of arbitrarily ordered isomers using the Kuhn-Munkres algorithm[J].Journal of Chemical Information and Modeling, 2017, 57(5):1045.

[12] FRANCOIS B, LASSALLE J C.An extension of the Munkres algorithm for the assignment problem to rectangular matrices[J].Communications of the ACM，1971, 14(12):802-804.

[13] ZHU H B, ZHOU M C, ALKINS R.Group role assignment via a Kuhn-Munkres algorithm-based solution[J].IEEE Transactions on Systems, Man, and Cybernetics-Part A: Systems and Humans, 2012, 42(3):739-750.

[14] PAMELA M.B, KEVIN B.Structure-preserving Runge-Kutta methods for stochastic Hamiltonian equations with additive noise[J].Numerical Algorithms.2014, 65(3):519-532.

[15] DANG Q A, MANH T H.Positive and elementary stable explicit nonstandard Runge-Kutta methods for a class of autonomous dynamical systems[J].International Journal of Computer Mathematics, 2020, 97(10):2036-2054.