• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    基于滑??刂频母呱璞蕊w行器協(xié)同制導(dǎo)方法

    2025-03-20 00:00:00郭博鐵鳴范文慧李傳旭
    關(guān)鍵詞:滑模控制參數(shù)優(yōu)化穩(wěn)定性

    摘 要:針對(duì)多個(gè)高升阻比飛行器在飛行中間段的時(shí)間協(xié)同問(wèn)題,提出一種基于滑模控制的高升阻比飛行器協(xié)同制導(dǎo)方法,建立多個(gè)高升阻比飛行器協(xié)同飛行力學(xué)模型,設(shè)計(jì)規(guī)劃飛行器在中間段制導(dǎo)過(guò)程中的飛行程序與針對(duì)飛行器側(cè)向機(jī)動(dòng)的滑??刂撇呗浴Mㄟ^(guò)設(shè)計(jì)標(biāo)稱軌跡,對(duì)控制參數(shù)進(jìn)行優(yōu)化,進(jìn)而生成各個(gè)飛行器不同初始條件所需的側(cè)向過(guò)載并通過(guò)側(cè)向過(guò)載得到所需的傾側(cè)角指令,以實(shí)現(xiàn)對(duì)飛行器飛行時(shí)間的控制,從而使多飛行器能夠同時(shí)到達(dá)目標(biāo)點(diǎn)上方設(shè)定范圍;考慮飛行器始末條件和狀態(tài)約束,使飛行器能夠滿足協(xié)同任務(wù)需求。利用李雅普諾夫穩(wěn)定性判據(jù)證明系統(tǒng)的穩(wěn)定性以及滑模面非奇異性,仿真結(jié)果表明,該協(xié)同制導(dǎo)策略具備一定抗干擾性,能夠滿足異地非同步發(fā)射的多個(gè)飛行器協(xié)同制導(dǎo)需求。

    關(guān)鍵詞: 高升阻比飛行器; 協(xié)同制導(dǎo); 滑??刂? 參數(shù)優(yōu)化; 穩(wěn)定性

    中圖分類號(hào): V 448.232

    文獻(xiàn)標(biāo)志碼: ADOI:10.12305/j.issn.1001 506X.2025.02.24

    Cooperative guidance method of high lift to drag ratio aircraft

    based on sliding mode control

    GUO Bo1, TIE Ming1,*, FAN Wenhui2, LI Chuanxu1

    (1. Science and Technology on Space Physics Laboratory, Beijing 100076, China;

    2. Department of Automation, Tsinghua University, Beijing 100084, China)

    Abstract:Aiming at the problem of time coordination of multiple high lift to drag ratio aircrafts in mid flight, a cooperative guidance method of high lift to drag ratio aircraft based on sliding mode control is proposed. A cooperative flight mechanics model of multiple lift to drag ratio aircrafts is established, and the flight program of the aircraft in mid flight guidance and the sliding mode control strategy for lateral maneuvering of the aircraft are designed and planned. By designing nominal trajectories, the control parameters are optimized to generate the lateral overload required by different initial conditions of each aircraft and obtain the required tilt angle command through the lateral overload, so as to control the flight time of the aircraft, so that multiple aircrafts can reach the set range above the target point at the same time, and consider the starting and ending conditions and state constraints of the aircraft, enabling the aircraft to meet collaborative mission requirements. Lyapunov stability criterion is used to prove the stability of the system and the non singularity of the sliding mode surface. Simulation results show that the cooperative guidance strategy has a certain anti interference property, and can meet the requirements of the cooperative guidance of multiple aircrafts launched asynchronously.

    Keywords:high lift to drag ratio aircraft; cooperative guidance; sliding mode control; parameter optimization; stability

    0 引 言

    高升阻比飛行器能夠在20~100 km的高層大氣飛行,相比一般的飛行器,該類飛行器飛行速度更快、飛行航程更遠(yuǎn)、飛行空域更廣,在一定時(shí)間內(nèi)改變飛行速度和飛行方向的能力更強(qiáng),因此研究該類飛行器具有重大意義1-2。然而,隨著工程需要,在實(shí)際場(chǎng)景中單個(gè)飛行器需逐漸勝任難度更高的任務(wù),多飛行器協(xié)同工作能力更強(qiáng),能夠勝任更高難度的任務(wù)。因此,多飛行器協(xié)同任務(wù)規(guī)劃在近些年來(lái)愈發(fā)受到關(guān)注3,需要考慮研究多飛行器協(xié)同的相關(guān)技術(shù)。

    對(duì)比單飛行器的任務(wù)規(guī)劃,多飛行器的協(xié)同規(guī)劃難度更高,需要考慮的約束更復(fù)雜,除了需要同時(shí)考慮多個(gè)飛行器的過(guò)載、動(dòng)壓、熱約束外4,還需要同時(shí)考慮在不同始末條件下多飛行器飛行時(shí)間的控制5

    實(shí)現(xiàn)多個(gè)飛行器針對(duì)同一目標(biāo)點(diǎn)的協(xié)同制導(dǎo)方式主要分為基于到達(dá)時(shí)間約束的飛行器獨(dú)立控制制導(dǎo)和基于信息交互的飛行器在線協(xié)同制導(dǎo)6-7。飛行器獨(dú)立時(shí)間控制制導(dǎo),即預(yù)先設(shè)定每個(gè)飛行器的飛行時(shí)間,在飛行過(guò)程中不進(jìn)行多個(gè)飛行器間的信息交互;基于信息交互的在線協(xié)同制導(dǎo)則利用多飛行器在飛行過(guò)程中的互相通信,實(shí)現(xiàn)時(shí)間協(xié)同8-9。而當(dāng)飛行器在大氣層高速飛行時(shí),飛行器周?chē)鷷?huì)形成等離子體鞘層,干擾通信10。因此,本文研究多飛行器的協(xié)同制導(dǎo),并采用更符合實(shí)際的獨(dú)立時(shí)間控制。

    針對(duì)多飛行器協(xié)同制導(dǎo)方面的研究,早期有些學(xué)者專門(mén)研究針對(duì)單個(gè)飛行器到達(dá)靜止目標(biāo)點(diǎn)的時(shí)間控制11-12。在此基礎(chǔ)上,又有更多協(xié)同制導(dǎo)方案被提出。文獻(xiàn)[13]首次將攻擊時(shí)間誤差反饋引入比例導(dǎo)引律,提出撞擊時(shí)間控制制導(dǎo)(impact time control guidance,ITCG),實(shí)現(xiàn)飛行器針對(duì)靜止目標(biāo)的協(xié)同。文獻(xiàn)[14]提出一種關(guān)于剩余時(shí)間和剩余距離的計(jì)算公式,通過(guò)預(yù)測(cè)剩余時(shí)間、剩余距離和末速度之間的關(guān)系建立能滿足到達(dá)時(shí)間約束的再入制導(dǎo)律。文獻(xiàn)[15]通過(guò)在橫向剖面使用預(yù)測(cè)校正制導(dǎo)調(diào)節(jié)時(shí)間以通過(guò)牛頓迭代法調(diào)節(jié)該剖面以滿足剩余距離和剩余時(shí)間的約束。王肖等人研究通過(guò)在高度-速度剖面內(nèi)設(shè)計(jì)參考軌跡,結(jié)合側(cè)向航向角走廊,通過(guò)在線數(shù)值預(yù)測(cè)剩余飛行航程和時(shí)間校正實(shí)現(xiàn)飛行器的時(shí)間協(xié)同制導(dǎo)16。周宏宇等人提出一種基于改進(jìn)粒子群優(yōu)化算法的多飛行器協(xié)同規(guī)劃算法,將強(qiáng)化學(xué)習(xí)方法同粒子群優(yōu)化算法結(jié)合,通過(guò)多次重構(gòu)飛行剖面最終完成協(xié)同制導(dǎo)17。You等人提出一種可以針對(duì)三維空間中多個(gè)飛行器的領(lǐng)航跟隨固定時(shí)間協(xié)同算法,使跟隨飛行器能夠通過(guò)算法,利用視線角度由領(lǐng)航飛行器帶領(lǐng)實(shí)現(xiàn)協(xié)同制導(dǎo)18。

    針對(duì)多飛行器協(xié)同控制,本文選取控制方式為滑??刂?。滑??刂苿?dòng)態(tài)響應(yīng)快,對(duì)參數(shù)攝動(dòng)不敏感,對(duì)外部的干擾具有較好的魯棒性,常被用來(lái)解決帶有碰撞角度約束的飛行器制導(dǎo)問(wèn)題19-20。文獻(xiàn)[21]提出一種考慮落角約束的滑模制導(dǎo)律,該方法通過(guò)改變分?jǐn)?shù)階參數(shù)實(shí)現(xiàn)對(duì)飛行器軌跡以及飛行時(shí)間的調(diào)節(jié)。文獻(xiàn)[22]和文獻(xiàn)[23]均為通過(guò)對(duì)法向過(guò)載使用滑??刂品椒ú⒃O(shè)定期望時(shí)間和期望角度,實(shí)現(xiàn)針對(duì)靜止目標(biāo)的協(xié)同制導(dǎo)。文獻(xiàn)[24]提出一種協(xié)同多個(gè)飛行器飛行時(shí)間和角度的三維自適應(yīng)滑模制導(dǎo)律,實(shí)現(xiàn)多飛行器在期望時(shí)刻的分布式協(xié)同制導(dǎo),并保證閉環(huán)系統(tǒng)的全局固定時(shí)間收斂特性。

    對(duì)于高升阻比飛行器來(lái)說(shuō),飛行器的動(dòng)力學(xué)模型都具有強(qiáng)非線性、強(qiáng)不確定性,其速度和加速度都是實(shí)時(shí)變化的。而且,相比末制導(dǎo)段,飛行器在中間段的飛行時(shí)間更長(zhǎng),距離更遠(yuǎn),有更足夠的空間來(lái)調(diào)節(jié)飛行時(shí)間25。因此,本文研究?jī)?nèi)容將進(jìn)一步考慮在氣動(dòng)力加速度隨大氣密度變化的情況下,多個(gè)飛行器在中間段的協(xié)同制導(dǎo),即多個(gè)飛行器從不同初始點(diǎn)開(kāi)始,依靠氣動(dòng)力通過(guò)滑模控制使飛行器在飛行時(shí)間相同的情況下到達(dá)同一目標(biāo)點(diǎn)上空,為末端制導(dǎo)做準(zhǔn)備。

    最后,通過(guò)仿真實(shí)驗(yàn)驗(yàn)證該方法的正確性和有效性。

    1 飛行器動(dòng)力學(xué)模型

    1.1 動(dòng)力學(xué)方程

    本文主要研究多飛行器在經(jīng)歷火箭助推段、自由段以及拉起段后在大氣層中無(wú)動(dòng)力飛行的中間段制導(dǎo)過(guò)程,假設(shè)第i枚飛行器的的動(dòng)力學(xué)模型可以表示如下:

    dvidt=-Dimi-gisin γi(1)

    dridt=visin γi(2)

    dNidt=vicos γisin iricos i(3)

    didt=vicos γicos iri(4)

    didt=Lisin σimivicos γi+Viricos γisinitan i(5)

    dγidt=Licos σimivi-givi-viricos γi(6)

    式中:下標(biāo)i代表第i枚飛行器;v為飛行器飛行速度;r為飛行器的地心距;N為飛行器所處經(jīng)度;為飛行器所處緯度;為飛行器航向角;γ為飛行器彈道傾角;D為飛行器受到的氣動(dòng)阻力;L為飛行器受到的氣動(dòng)升力;m為飛行器質(zhì)量;g為重力加速度大小;σ為飛行器傾側(cè)角。

    氣動(dòng)升力和氣動(dòng)阻力可表示如下:

    L=12ρv2SrefCL

    D=12ρv2SrefCD(7)

    式中:Sref為特征面積;CL和CD分別為升力系數(shù)和阻力系數(shù),可以分別由飛行器的攻角和馬赫數(shù)通過(guò)擬合計(jì)算得到26。

    h和ρ分別為飛行器所在高度和所在高度處的大氣密度:

    ρ=ρ0e-βh(8)

    h=r-R0(9)

    式中:ρ0為標(biāo)準(zhǔn)大氣密度;β為一常數(shù);R0為地球的平均半徑。

    取ρ0=1.225 kg/m3;β=1.406 4×10-4;R0=6 374.4 km。

    1.2 飛行約束條件

    每一枚飛行器在大氣層高速飛行的過(guò)程中,首先需要滿足氣動(dòng)熱約束、過(guò)載約束和動(dòng)壓約束:

    Qi=KQρiv3.15i≤Qmax(10)

    qi=0.5ρiv2i≤qmax(11)

    ni=L2i+D2imigi≤nmax(12)

    KQ=C1Rdρ0(g0R0)3.15(13)

    式中:Q為熱流密度;q為動(dòng)壓;n為總過(guò)載;KQ為熱流常數(shù);C1為形狀因子;Rd為前緣駐點(diǎn)處半徑;g0為地球表面加速度;Qmax為最大熱流密度;qmax為最大動(dòng)壓;nmax為最大過(guò)載。

    除此之外,還要求飛行器在滿足到達(dá)終端位置約束、終端高度約束、總航程約束的同時(shí),還需要保證多個(gè)飛行器到達(dá)的時(shí)間相同。

    設(shè)Srest,i為第i枚飛行器的剩余待飛航程,有以下關(guān)系式27

    Srest,i=R0·arccos[cos icos ′cos(θi-θ′)+

    sin isin ′](14)

    式中:i,θi分別表示第i枚飛行器當(dāng)前位置的緯度和經(jīng)度;′,θ′分別表示終點(diǎn)位置的緯度和經(jīng)度。

    設(shè)tf,i表示第i枚飛行器的到達(dá)時(shí)間,tdf為各個(gè)飛行器的期望到達(dá)時(shí)間,需要滿足關(guān)系:

    tf,1=tf,2=…=tf,i=tdf(15)

    2 多飛行器協(xié)同制導(dǎo)設(shè)計(jì)

    2.1 滑模面設(shè)計(jì)

    針對(duì)靜止目標(biāo)的制導(dǎo)幾何可以繪制在如圖1所示的參考系中。λ為飛行器相對(duì)目標(biāo)的視線角;γ為飛行器的航向角;θ表示航向角和視線角的差值,即飛行器前置角,前置角可以描述飛行器飛行方位角與目標(biāo)點(diǎn)的偏差,由此可以寫(xiě)成如下關(guān)系式:

    θ=λ-γ(16)

    vr=vcos θ(17)

    vλ=vsin θ(18)

    還可以由此得到視線角和航向角和速度分量vλ以及側(cè)向過(guò)載aM之間的關(guān)系:

    λ·=vλr(19)

    γ·=aMv(20)

    對(duì)于飛行器無(wú)動(dòng)力飛行的中段制導(dǎo)過(guò)程,飛行器的目標(biāo)是到達(dá)預(yù)定位置并與末制導(dǎo)段進(jìn)行交班。交班區(qū)域位置設(shè)計(jì)為在以目標(biāo)點(diǎn)經(jīng)緯度為圓心的圓形范圍內(nèi),當(dāng)飛行器到達(dá)該圓形范圍內(nèi)即完成任務(wù),飛行器還需在飛行末時(shí)刻滿足高度約束。

    針對(duì)多飛行器協(xié)同制導(dǎo),還需滿足多飛行器到達(dá)交班點(diǎn)的時(shí)間相同。因此,本文采用的是預(yù)測(cè)-校正制導(dǎo)方法,通過(guò)剩余飛行距離預(yù)估時(shí)間對(duì)飛行控制指令進(jìn)行實(shí)時(shí)校正以實(shí)現(xiàn)對(duì)飛行時(shí)間的控制28。

    假設(shè)交班點(diǎn)與終點(diǎn)的所剩距離為S;飛行器當(dāng)前位置距終點(diǎn)位置為S;飛行器到交班點(diǎn)所需的預(yù)估剩余時(shí)間為tgo;預(yù)估到達(dá)時(shí)刻為tf;當(dāng)前時(shí)刻為t;則可以得到如下關(guān)系式29

    tf=t+tgo(21)

    Srest=S-S(22)

    tgo=Srestv1+0.52N-1θ2(23)

    式中:N為比例導(dǎo)引系數(shù),一般設(shè)N=3??梢杂纱嗽O(shè)計(jì)滑模面為

    s=tf-tdf=t+tgo-tdf=tgo-tdgo(24)

    式中:tdf表示期望到達(dá)時(shí)間;tdgo表示期望剩余時(shí)間?;C鎠最終趨于穩(wěn)定。

    2.2 滑模制導(dǎo)律設(shè)計(jì)

    令r=Srest,通過(guò)對(duì)滑模面求導(dǎo),設(shè)計(jì)滑模制導(dǎo)律函數(shù)30

    s·=tgo·-tdgo·=r·v1+0.52N-1θ2-rθv(2N-1)λ·+

    1+rθv2(2N-1)am(25)

    基于式(25),可以利用李雅普諾夫非線性控制理論設(shè)計(jì)滑模控制的趨近律。

    aM=aeqM+aMconM+aswM(26)

    式中:aM為控制飛行器所需要的側(cè)向過(guò)載;aeqM為利用s·=0設(shè)計(jì)的等效控制分量,aeqM可以表示為

    aeqM=v2(2N-1)r·cos θ-1θ+0.5v2rθcos θ+VMλ·(27)

    同時(shí),為了使滑模面最終趨于零并保證系統(tǒng)滿足李雅普諾夫穩(wěn)定性條件,設(shè)計(jì)aswM和aMconM,可以寫(xiě)成如下關(guān)系式:

    aswM=-M(psgn(θ)+1)sgn(s)(28)

    aMconM=-h(huán)(θ)θks(29)

    在這里,可以引入一連續(xù)正函數(shù):

    h(θ)=sgn(θ)θ,|θ|lt;ε1

    1-ε1ε2-ε1|θ|+ε1(ε2-1)ε2-ε1, ε1≤|θ|≤ε2

    1, 其他(30)

    式中:k,s,M,p,ε1,ε2均為正常數(shù)且需要滿足pgt;1。

    接下來(lái),對(duì)系統(tǒng)的穩(wěn)定性進(jìn)行證明,證明過(guò)程可參考文獻(xiàn)[30],考慮李雅普諾夫函數(shù)為

    V=12s2(31)

    上述滑模面和趨近律的設(shè)計(jì),可以保證李雅普諾夫函數(shù)是正定的,接下來(lái)只需要證明其導(dǎo)數(shù)是負(fù)定,即可保證系統(tǒng)是穩(wěn)定的。

    對(duì)李雅普諾夫函數(shù)求導(dǎo)可得

    V·=ss·=-rh(θ)v2(2N-1)ks2-Mrv2(2N-1)(p|θ|+θ)|s|(32)

    分別取式(32)中兩項(xiàng)單獨(dú)討論,令

    A=rh(θ)v2(2N-1)(33)

    B=Mrv2(2N-1)(p|θ|+θ)|s|(34)

    由于h(θ)的函數(shù)值域均為正,k,M,p為正常數(shù)且pgt;1,所以當(dāng)θ≠0時(shí),A和B都是正定的。這意味著當(dāng)θ=0時(shí)不能滿足趨近律s=0,因此需要證明θ=0不是一個(gè)穩(wěn)定吸引子。

    下面考慮θ=0的情況,由式(16)求導(dǎo)得到:

    θ·=λ·-γ·(35)

    由式(18)~式(20)和式(26)可知,當(dāng)θ=0時(shí),能夠滿足λ·=0,aeqM=0,aMconM=0,由此可以得到:

    γ·=aswMv=-Mvsgn(s)(36)

    將式(36)和式(35)聯(lián)立得到:

    θ·=aswMv=Mvsgn(s)(37)

    因此,當(dāng)θ=0且s≠0時(shí),θ·≠0,系統(tǒng)并沒(méi)有趨于穩(wěn)定,所以θ=0并非一個(gè)穩(wěn)定性因子。僅當(dāng)s=0時(shí),θ=0才是一個(gè)吸引子。因此,由式(32)可以證明系統(tǒng)滿足李雅普諾夫穩(wěn)定性條件,接下來(lái)討論制導(dǎo)指令產(chǎn)生奇異性的問(wèn)題。理論上,當(dāng)r為零或無(wú)限接近于零時(shí)制導(dǎo)指令會(huì)產(chǎn)生奇異性現(xiàn)象,然而在實(shí)際飛行過(guò)程中,一般會(huì)將飛行結(jié)束條件設(shè)置為進(jìn)入S*范圍內(nèi)。此時(shí),Srest負(fù)定,Srest為零的狀態(tài)只是到達(dá)結(jié)束條件下的某一個(gè)短暫時(shí)刻,即不會(huì)在滑模面趨于穩(wěn)態(tài)時(shí)連續(xù)為零。因此,在式(26)中的滑模制導(dǎo)律不會(huì)導(dǎo)致奇異產(chǎn)生。同時(shí),另一個(gè)奇異條件是θ等于零或無(wú)限接近于零,而在式(27)中,令:

    T=cos θ-1θ(38)

    依據(jù)洛必達(dá)法則,當(dāng)θ無(wú)限趨于零的時(shí)候,T也是趨于零的,而不會(huì)變成無(wú)窮大項(xiàng),所以θ不會(huì)導(dǎo)致奇異情況產(chǎn)生。

    同樣,可以證明所需過(guò)載aM的另外兩個(gè)分量aMconM和aswM不會(huì)出現(xiàn)奇異,本文設(shè)計(jì)的滑模制導(dǎo)律不會(huì)出現(xiàn)奇異問(wèn)題。

    另外,為了抑制控制系統(tǒng)在運(yùn)行過(guò)程中產(chǎn)生的抖振現(xiàn)象,現(xiàn)將不連續(xù)函數(shù)sgn(x)改寫(xiě)為一連續(xù)函數(shù)sgmf(x):

    sgmf(x)=211+e-αz-12, αgt;0(39)

    2.3 參數(shù)優(yōu)化求解

    由于飛行器飛行速度、加速度都是實(shí)時(shí)變化的,所以為了防止飛行過(guò)程中傾側(cè)角變化出現(xiàn)抖振、能量不足等原因?qū)е聼o(wú)法在限定時(shí)間到達(dá)終端位置,可以對(duì)參數(shù)M進(jìn)行尋優(yōu)計(jì)算?,F(xiàn)采用提前規(guī)劃標(biāo)稱軌跡的方式,利用遞推最小二乘法對(duì)合適的M值進(jìn)行計(jì)算,其最小二乘指標(biāo)為

    J=min(tdf-tf,i)(40)

    本文采用遞推最小二乘法迭代求解合適的M值,利用飛行時(shí)間反向計(jì)算所需要的M值,以下給出使用該方法的計(jì)算過(guò)程。

    設(shè)計(jì)變量:

    ξ=M(41)

    建立如下殘差方程:

    g(ξ)=tdf-tf(42)

    由最小二乘法的極值原理,可以得到如下的迭代算法:

    ξk+1=ξk-A-1kBk(43)

    式中:

    Ak=∑Ni=1gTiξkgiξk(44)

    Bk=∑Ni=1gTiξkgi(ξk)(45)

    當(dāng)估計(jì)值與上一步的差的范數(shù)滿足一定精度時(shí)就停止迭代。

    2.4 縱向與側(cè)向軌跡規(guī)劃

    得到側(cè)向過(guò)載指令aM后,再由aM求得所需的程序角指令,由于飛行器采用傾側(cè)轉(zhuǎn)彎(bank to turn,BTT)控制方法,側(cè)向控制指令均由傾側(cè)角σ產(chǎn)生,根據(jù)側(cè)向過(guò)載aM,升力L和飛行器質(zhì)量m得到傾側(cè)角指令為

    σ=arcsin(aMmL)(46)

    由傾側(cè)角指令產(chǎn)生的側(cè)向過(guò)載是有限的,所以側(cè)向過(guò)載aM不能無(wú)限制地施加。另外,飛行器需要維持穩(wěn)定,因此需要對(duì)傾側(cè)角限幅,故限定傾側(cè)角范圍為σ∈[-60°,60°]。

    除傾側(cè)角外,飛行器的程序角還包括攻角,為滿足約束條件,還需要對(duì)飛行器的攻角α進(jìn)行規(guī)劃:

    α=amax, v≥v1

    aL/D max+amax-aL/D maxv1-v2(v-v2), v2lt;vlt;v1

    aL/D max, v≤v2(47)

    式中:amax為飛行器攻角的最大值;αL/D max為飛行器最大升阻比情況下的攻角,可由升阻比函數(shù)CL/CD求出;v1為一給定的較大速度值;v2為一給定的較小速度值。

    因此,所研究的問(wèn)題可以進(jìn)一步轉(zhuǎn)化為通過(guò)利用航程-傾側(cè)角剖面以及速度-攻角剖面設(shè)計(jì)的飛行走廊。

    飛行器的軌跡通過(guò)四階五步龍格庫(kù)塔法積分得到,每0.1 s進(jìn)行一次積分,當(dāng)滿足終端位置約束時(shí)結(jié)束積分,假設(shè)終端位置為進(jìn)入距離目標(biāo)點(diǎn)半徑為S的圓區(qū)域內(nèi)時(shí)即為飛行器完成任務(wù),最優(yōu)指標(biāo)為每一枚飛行器的到達(dá)時(shí)間一致且均為期望到達(dá)時(shí)間tfd,即

    J=min(tfd-tf,i), i=1,2,…,N(48)

    式中:k為飛行器的數(shù)量。

    3 仿真分析

    3.1 仿真背景及參數(shù)設(shè)定

    假設(shè)共有3枚飛行器且該3枚飛行器均為同種類的飛行器,各個(gè)飛行器的質(zhì)量均為1 000 kg,特征面積均設(shè)計(jì)為0.51 m2,氣動(dòng)參數(shù)參考使用CAV H(common aero vehicle high)飛行器的氣動(dòng)參數(shù),該類飛行器的過(guò)程約束可以參考文獻(xiàn)[31]并設(shè)置熱流上限為Qmax=1 200 kw/m2,過(guò)載上限為Nmax=4 g,動(dòng)壓上限為q=400 kPa。

    假設(shè)3枚飛行器(i=1,2,3)在經(jīng)過(guò)前面階段的飛行后于同一時(shí)刻開(kāi)始中間段飛行,分別由不同初始位置出發(fā),其初始條件分別如表1所示,而飛行器的目標(biāo)點(diǎn)與終端約束情況如表2所示。

    由表1可知,設(shè)定的3枚飛行器初始條件不同,因此其剩余航程均不相同。3枚飛行器剩余航程分別為4 939.4 km,5 062.0 km,5 005.2 km。飛行器初始前置角均為零,因此在初速度、初始高度、初始航跡傾角都相近的情況下,如果不通過(guò)合適側(cè)向機(jī)動(dòng)改變飛行時(shí)間,3枚飛行器最終到達(dá)終點(diǎn)的時(shí)間必不相同。

    3.2 多飛行器協(xié)同制導(dǎo)仿真分析

    對(duì)飛行器制導(dǎo)律的參數(shù)指標(biāo)進(jìn)行設(shè)計(jì),設(shè)k=0.001,p=1.3,ε1=0.001,ε2=0.015。

    對(duì)于期望到達(dá)時(shí)間tdf的設(shè)定有一定局限范圍,既不能過(guò)大也不能過(guò)小。因?yàn)榧纫WC其大于剩余航程最長(zhǎng)的飛行器在無(wú)側(cè)向機(jī)動(dòng)情況下的飛行時(shí)間,還要保證不能因飛行時(shí)間過(guò)長(zhǎng)導(dǎo)致能量提前消耗殆盡以致其無(wú)法到達(dá)終點(diǎn)。經(jīng)計(jì)算,能夠滿足3枚飛行器均可到達(dá)的期望時(shí)間可取范圍為tdf∈[1 366,1 404]s,因此設(shè)置tdf=1 386 s。通過(guò)第2節(jié)的最小二乘法計(jì)算出3枚飛行器的控制參數(shù)M值分別如表3所示。

    利用所設(shè)計(jì)的指標(biāo)進(jìn)行仿真,得到3枚飛行器在不同初始條件下的協(xié)同制導(dǎo)仿真結(jié)果。圖2展示了規(guī)劃出3枚飛行器的航程-高度飛行剖面,3枚飛行器在經(jīng)過(guò)不斷的跳躍飛行最終能夠完成所需的全部航程。圖3展示了3枚飛行器經(jīng)緯度變化的軌跡,可以看出飛行器由不同初始位置出發(fā),在經(jīng)過(guò)不斷側(cè)向機(jī)動(dòng)后最終能夠到達(dá)終點(diǎn)區(qū)域。

    圖4和圖5展示的是3枚飛行器飛行過(guò)程中高度和速度的變化。飛行器高度在不斷的來(lái)回波動(dòng)中緩慢下降至終端約束的高度范圍內(nèi),而飛行器速度隨著能量的消耗不斷減小。

    圖6展示了3枚飛行器飛行過(guò)程中動(dòng)壓、過(guò)載、熱流的變化情況,這些變化量在飛行過(guò)程中始終符合過(guò)程約束,不會(huì)超過(guò)限定值。

    圖7展示了3枚飛行器由制導(dǎo)指令產(chǎn)生的傾側(cè)角變化過(guò)程,為了改變飛行時(shí)間而調(diào)整飛行方向,傾側(cè)角需要根據(jù)制導(dǎo)指令在限定范圍內(nèi)實(shí)時(shí)調(diào)整。

    圖8展示了3枚飛行器在飛行過(guò)程中攻角的變化情況,為維持飛行高度以及實(shí)現(xiàn)更大航程,飛行器在飛行的過(guò)程中需要調(diào)整攻角,使攻角隨速度的衰減不斷增大。

    圖9和圖10分別展示了3枚飛行器的航向角變化和前置角變化,航向角和前置角在不斷地側(cè)向調(diào)整中,逐漸恢復(fù)到初始時(shí)的角度。

    3枚飛行器初始的航向角就是初始點(diǎn)與目標(biāo)點(diǎn)之間的方位角,因此在飛行器不進(jìn)行任何側(cè)向飛行的情況下,飛行器剛好能以最短時(shí)間準(zhǔn)確飛達(dá)目標(biāo)點(diǎn)上空。同理,初始前置角均接近于零,這意味著飛行器初始飛行方向是準(zhǔn)確朝著目標(biāo)點(diǎn)的,但是由于飛行器需要依靠側(cè)向機(jī)動(dòng)調(diào)整飛行時(shí)間,所以要進(jìn)行飛行方向的調(diào)整,而側(cè)向機(jī)動(dòng)的指令則由計(jì)算得到的傾側(cè)角實(shí)現(xiàn)。在完成對(duì)時(shí)間的調(diào)整后,前置角最終會(huì)逐漸恢復(fù)到零。航向角同樣如此,在完成對(duì)飛行時(shí)間的調(diào)整后也會(huì)逐漸恢復(fù)到初始航向角大小,進(jìn)而繼續(xù)準(zhǔn)確朝著目標(biāo)點(diǎn)進(jìn)發(fā)。

    最終,3枚飛行器均能在設(shè)定時(shí)間到達(dá)指定目標(biāo)區(qū)域內(nèi),其終端狀態(tài)如表4所示??梢?jiàn),飛行器均到達(dá)了距離目標(biāo)點(diǎn)140 km范圍處,3枚飛行器與期望到達(dá)時(shí)間的偏差分別為0.5 s、0.5 s、1 s??梢哉J(rèn)為在誤差范圍內(nèi),飛行器協(xié)同制導(dǎo)的任務(wù)是成功的。

    3.3 對(duì)比分析

    本文研究方法的思路主要是通過(guò)側(cè)向機(jī)動(dòng)以達(dá)到改變飛行器飛行時(shí)間的目的。而文獻(xiàn)[21]提出的方法主要針對(duì)基于縱向落角約束的預(yù)測(cè)-校正制導(dǎo)研究。因此,對(duì)比該方法的仿真如圖11~圖16所示。

    從結(jié)果看,文獻(xiàn)[21]方法通過(guò)預(yù)估的末時(shí)刻落角反饋制導(dǎo)指令,改變縱向飛行剖面。使用該方法同樣能夠達(dá)到協(xié)同制導(dǎo)的目的,并使飛行器在約束范圍內(nèi)到達(dá)距離目標(biāo)點(diǎn)140 km范圍內(nèi),完成交班任務(wù),最終的到達(dá)時(shí)間與期望到達(dá)時(shí)間的偏差分別為0.7 s、0.1 s、0.5 s,可以認(rèn)為在誤差范圍內(nèi)。

    但是,該方法通過(guò)改變制導(dǎo)律中的分?jǐn)?shù)階參數(shù),僅僅只改變飛行器的縱向軌跡,沒(méi)有通過(guò)橫側(cè)向制導(dǎo)改變飛行器的橫側(cè)向軌跡,因此飛行器只可始終沿初始航向角飛行。而飛行器在實(shí)際飛行中難免受到干擾,不能在無(wú)側(cè)向機(jī)動(dòng)控制的情況下始終保持初始方位角,最終可能造成飛行方向的偏差。

    相比文獻(xiàn)[21]中的方法,本文研究方法更優(yōu)異的方面在于能夠通過(guò)直接對(duì)飛行器的橫側(cè)向軌跡進(jìn)行控制約束,使飛行器在飛行過(guò)程中糾正到達(dá)時(shí)間偏差的同時(shí),還可以對(duì)飛行方向進(jìn)行修正,最終令飛行器始終能夠在制導(dǎo)指令下向目標(biāo)位置飛行。

    除此之外,利用與本文研究方法同樣的思路,使用傳統(tǒng)的比例-積分-微分(proportion integration differentiation,PID)控制器法設(shè)置一組對(duì)照實(shí)驗(yàn)。該方法與本文研究思路類似,同樣利用預(yù)估剩余時(shí)間tgo和期望剩余時(shí)間tdgo設(shè)計(jì)關(guān)于側(cè)向過(guò)載的PID制導(dǎo)律:

    aM=KP(tgo-tdgo)+KI∫tf0(tgo-tdgo)dt+

    KD(tgo·-tdgo·)(49)

    式中:KP、KI、KD分別為PID控制參數(shù)。利用與前面所述處理滑??刂茀?shù)相同的方法,固定KI、KD,通過(guò)最小二乘迭代計(jì)算調(diào)節(jié)參數(shù)KP。

    圖17~圖19分別展示了在傳統(tǒng)PID控制方法下的3枚飛行器航程-高度飛行剖面、過(guò)程約束以及傾側(cè)角指令的變化。由圖17可以看出,3枚飛行器最終在飛行終點(diǎn)處均能夠到達(dá)指定范圍內(nèi)的高度。由圖18可以看出,3枚飛行器均能在過(guò)程約束范圍內(nèi)飛行。但是,通過(guò)圖19得到,3枚飛行器的傾側(cè)角指令在飛行最后時(shí)段發(fā)生了在上下幅值之間的高頻率來(lái)回抖振,這意味著需要控制傾側(cè)角變化的飛行器舵發(fā)生高頻率來(lái)回?cái)[動(dòng),顯然是不符合實(shí)際的。另外,在使用該傳統(tǒng)PID控制方法得到的結(jié)果中,3枚飛行器到達(dá)時(shí)刻分別為1 388.2 s、1 387.4 s、1 386.0 s,與期望到達(dá)時(shí)間的最大偏差達(dá)到了2.2 s,結(jié)果指標(biāo)略遜于本文所述方法。因此,通過(guò)仿真實(shí)驗(yàn)可知,本文使用方法相比于傳統(tǒng)PID控制方法要更優(yōu)。

    3.4 蒙特卡羅仿真分析

    在實(shí)際飛行過(guò)程中,飛行器會(huì)受到各種不確定性因素的干擾,所以需要對(duì)飛行器在受到外界擾動(dòng)情況下的飛行條件進(jìn)行仿真驗(yàn)證。因此,為了進(jìn)一步驗(yàn)證該方法的抗干擾性,選取一組數(shù)據(jù)利用蒙特卡羅打靶法進(jìn)行仿真。該組數(shù)據(jù)為滿足正態(tài)分布條件的初始偏差和外界擾動(dòng),如表5所示,偏差限為3σ。

    接下來(lái)以飛行器3為例,驗(yàn)證本文研究方法中飛行器在微小擾動(dòng)情況下的實(shí)際效果。通過(guò)100次蒙特卡羅打靶仿真,得到如下結(jié)果。

    圖20為飛行器在擾動(dòng)條件下高度隨時(shí)間變化情況,飛行器飛行過(guò)程中高度受到影響,產(chǎn)生微小變化,但是在終點(diǎn)仍能符合預(yù)定高度約束范圍。圖21為飛行器到達(dá)時(shí)間與目標(biāo)到達(dá)時(shí)間tdf偏差的密度分布,在滿足正態(tài)分布的偏差條件下,飛行器到達(dá)時(shí)間偏差同樣符合正態(tài)分布,與預(yù)期情況一致。但是正態(tài)分布的均值不為零,原因是最初由最小二乘法得到的最優(yōu)控制參數(shù)M計(jì)算出的飛行器到達(dá)時(shí)間本就有微小偏差,而此正態(tài)分布的均值即為最初到達(dá)時(shí)間與目標(biāo)到達(dá)時(shí)間tdf的差。但是,該正態(tài)分布的時(shí)間偏差的最大值為1.7 s,仍然在可以接受的誤差范圍內(nèi)。因此,可以認(rèn)為該多飛行器協(xié)同制導(dǎo)方法能克服干擾因素的影響。

    圖22展示了飛行器到達(dá)時(shí)間偏差與其對(duì)應(yīng)的末速度,直觀展示了飛行器末速度隨時(shí)間偏差的變化,這些偏差均在可接受范圍內(nèi)。

    4 結(jié) 論

    針對(duì)多高升阻比飛行器的協(xié)同制導(dǎo)技術(shù)是未來(lái)復(fù)雜環(huán)境與任務(wù)亟需發(fā)展的重點(diǎn)研究方向。本文提出一種通過(guò)預(yù)估剩余時(shí)間和預(yù)估剩余航程的多飛行器協(xié)同制導(dǎo)方法,該方法可以在保證制導(dǎo)準(zhǔn)確的同時(shí)調(diào)節(jié)飛行時(shí)間,實(shí)現(xiàn)多飛行器協(xié)同制導(dǎo)。

    通過(guò)設(shè)計(jì)基于滑模制導(dǎo)律并利用最小二乘迭代法計(jì)算控制參數(shù),規(guī)劃了飛行器飛行所需攻角和傾側(cè)角指令。通過(guò)攻角指令維持飛行高度,使飛行器能在到達(dá)終端狀態(tài)時(shí)滿足高度約束,通過(guò)傾側(cè)角指令使飛行器能有效調(diào)節(jié)飛行時(shí)間,解決復(fù)雜條件下對(duì)飛行器飛行時(shí)間控制的問(wèn)題,使飛行器能在多約束條件下實(shí)現(xiàn)協(xié)同制導(dǎo)。最后,通過(guò)仿真實(shí)驗(yàn)驗(yàn)證了該方法的可行性。仿真結(jié)果表明,多個(gè)飛行器的到達(dá)時(shí)間與預(yù)期到達(dá)時(shí)間的差距均能控制在1 s內(nèi)。對(duì)比其他方法和思路,本文研究方法具有更優(yōu)異的表現(xiàn)。并且,通過(guò)蒙特卡羅仿真實(shí)驗(yàn)可得,即使在有初始條件偏差和過(guò)程擾動(dòng)的情況下,飛行器的時(shí)間偏差仍然不超過(guò)2 s,說(shuō)明本方法具有一定抗干擾性,為高升阻比飛行器協(xié)同制導(dǎo)問(wèn)題的解決提供了一種有效的技術(shù)手段。

    參考文獻(xiàn)

    [1]DING Y B, YUE X K, CHEN G S, et al. Review of control and guidance technology on hypersonic vehicle[J]. Chinese Journal of Aeronautics, 2022, 35(7): 1-18.

    [2]WU T C, WANG H L, LIU Y H, et al. Learning based interfered fluid avoidance guidance for hypersonic reentry vehicles with multiple constraints[J]. ISA Transactions, 2023, 139(8): 291-307.

    [3]JIA X, WU S T, WEN Y M, et al. A distributed decision method for missiles autonomous formation based on potential game[J]. Journal of Systems Engineering and Electronics, 2019, 30(4): 738-748.

    [4]LIANG Z X, REN Z, LI Q D. Evolved atmospheric entry corridor with safety factor[J]. Acta Astronautica, 2018, 143(2): 82-91.

    [5]ZHANG D, LIU L, WANG Y J. On line reentry guidance algorithm with both path and no fly zone constraints[J]. Acta Astronautica, 2015, 117(12): 243-253.

    [6]YOU H, CHANG X L, ZHAO J F, et al. Three dimensional impact angle constrained cooperative guidance strategy against maneuvering target[J]. ISA Transactions, 2023, 138(7): 262-280.

    [7]ZHAO J, ZHOU R, DONG Z N. Three dimensional cooperative guidance laws against stationary and maneuvering targets[J]. Chinese Journal of Aeronautics, 2015,28(4): 1104-1120.

    [8]SONG J H, SONG S M, XU S L. Three dimensional cooperative guidance law for multiple missiles with finite time convergence[J]. Aerospace Science and Technology, 2017, 67(8): 193-205.

    [9]WANG X H, LU X. Three dimensional impact angle constrai ned distributed guidance law design for cooperative attacks[J]. ISA Transactions, 2018, 73(2): 79-90.

    [10]LI J F, ZHOU J X, YAO J F, et al. Experimental observations of communication in blackout,topological waveguiding and Dirac zero index property in plasma sheath[J]. Nanophotonics, 2023, 12(10): 1847-1856.

    [11]LEE J I, JEON I S, TAHK M J. Guidance law to control impact time and angle[J]. IEEE Trans.on Aerospace and Electronic Systems, 2007, 43(1): 301-310.

    [12]KIM T H, LEE C H, JEON I S, et al. Augmented polynomial guidance with impact time and angle constraints[J]. IEEE Trans.on Aerospace and Electronic Systems, 2013, 49(4): 2806-2817.

    [13]JEON I S, LEE J I, TAHK M J. Impact time control guidance law for anti ship missiles[J]. IEEE Trans.on Control Systems Technology, 2006, 14(2): 260-266.

    [14]GUO Y H, LI X, ZHANG H J, et al. Entry guidance with terminal time control based on quasi equilibrium glide condition[J]. IEEE Trans.on Aerospace and Electronic Systems, 2019, 56(2): 887-896.

    [15]LI Z H, HE B, WANG M H, et al. Time coordination entry guidance for multi hypersonic vehicles[J]. Aerospace Science and Technology, 2019, 89(6): 123-135.

    [16]王肖, 郭杰, 唐勝景, 等. 基于解析剖面的時(shí)間協(xié)同再入制導(dǎo)[J]. 航空學(xué)報(bào), 2019, 40(3): 239-250.

    WANG X, GUO J, TANG S J, et al. Time cooperative entry guidance based on analytical profile[J]. Acta Aeronautica et Astronautica Sinica, 2019, 40(3): 239-250.

    [17]周宏宇, 王小剛, 單永志, 等. 基于改進(jìn)粒子群算法的飛行器協(xié)同軌跡規(guī)劃[J]. 自動(dòng)化學(xué)報(bào), 2022, 48(11): 2670-2676.

    ZHOU H Y, WANG X G, SHAN Y Z, et al. Synergistic path planning for multiple vehicles based on an improved particle swarm optimization method[J]. Acta Automatica Sinica, 2022, 48(11): 2670-2676.

    [18]YOU H, CHANG X L, ZHAO J F, et al. Three dimensional impact angle constrained fixed time cooperative guidance algorithm with adjustable impact time[J]. Aerospace Science and Technology, 2023, 141: 108574.

    [19]ZHANG W J, FU S N, LI W, et al. An impact angle constraint integral sliding mode guidance law for maneuvering targets interception[J]. Journal of Systems Engineering and Electronics, 2020, 31(1): 168-184.

    [20]ZHAO F J, YOU H. New three dimensional secondorder sliding mode guidance law with impact angle constraints[J]. The Aeronautical Journal, 2020, 124(1273): 368-384.

    [21]盛永智, 甘佳豪, 張成新. 彈道可調(diào)的落角約束分?jǐn)?shù)階滑模制導(dǎo)律設(shè)計(jì)[J]. 航空學(xué)報(bào), 2023, 44(7): 182-195.

    SHENG Y Z, GAN J H, ZHANG C X. Fractional order sliding mode guidance law design with trajectory adjustable and terminal angular constraint[J]. Acta Aeronautica et Astronautica Sinica, 2023, 44(7): 182-195.

    [22]LI W, WEN Q Q, HE L, et al. Three dimensional impact angle constrained distributed cooperative guidance law for anti ship missiles[J]. Journal of Systems Engineering and Electro nics, 2021, 32(2): 447-459.

    [23]LYU T, LI C J, GUO Y N, et al.Three dimensional finite time cooperative guidance for multiple missiles without radial velocity measurements[J]. Chinese Journal of Aeronautics, 2019, 32(5): 1294-1304.

    [24]王雨辰, 王偉, 林時(shí)堯, 等. 考慮攻擊時(shí)間及空間角度約束的三維自適應(yīng)滑模協(xié)同制導(dǎo)律設(shè)計(jì)[J]. 兵工學(xué)報(bào), 2023, 44(9): 2778-2790.

    WANG Y C, WANG W, LIN S Y, et al. Three dimensional adaptive sliding mode cooperative guidance law with impact time and angle constraints[J]. Acta Armamentarii, 2023, 44(9): 2778-2790.

    [25]BAO C Y, WANG P, TANG G J. Integrated method of gui dance,control and morphing for hypersonic morphing vehicle in glide phase[J]. Chinese Journal of Aeronautics, 2021, 34(5): 535-553.

    [26]SHENG Y Z, ZHANG Z, XIA L. Fractional order sliding mode control based guidance law with impact angle constraint[J]. Nonlinear Dynamics, 2021, 106(1): 425-444.

    [27]姜鵬, 郭棟, 韓亮, 等. 多飛行器再入段時(shí)間協(xié)同彈道規(guī)劃方法[J]. 航空學(xué)報(bào), 2020, 41(S1): 171-183.

    JIANG P, GUO D, HAN L, et al. Trajectory optimization for cooperative reentry of multiple hypersonic glide vehicle[J]. Acta Aeronautica et Astronautica Sinica, 2020, 41(S1): 171-183.

    [28]LIANG Z X, LV C, ZHU S. Lateral entry guidance with terminal time constraint[J]. IEEE Trans.on Aerospace and Electronic System, 2023, 59(3): 2544-2553.

    [29]JEON I S, LEE J I, TAHK M J. Homing guidance law for cooperative attack of multiple missiles[J]. Journal of Guidance, Control, and Dynamics, 2010, 33(1): 275-280.

    [30]CHO D, KIM H J, TAHK M J. Nonsingular sliding mode guidance for impact time control[J]. Journal of Guidance, Control, and Dynamics, 2016, 39(1): 61-68.

    [31]徐慧, 蔡光斌, 崔亞龍, 等. 高超聲速滑翔飛行器再入軌跡優(yōu)化[J]. 哈爾濱工業(yè)大學(xué)學(xué)報(bào), 2023, 55(4): 44-55.

    XU H, CAI G B, CUI Y L, et al. Reentry trajectory optimization method of hypersonic glide vehicle[J]. Journal of Harbin Institute of Technology, 2023,55(4): 44-55.

    作者簡(jiǎn)介

    郭 博(1999—),男,碩士研究生,主要研究方向?yàn)轱w行器制導(dǎo)與控制、飛行器智能決策規(guī)劃與仿真。

    鐵 鳴(1976—),男,研究員,博士,主要研究方向?yàn)橹悄茱w行器設(shè)計(jì)、復(fù)雜系統(tǒng)建模與仿真。

    范文慧(1966—),男,教授,博士,主要研究方向?yàn)閺?fù)雜系統(tǒng)建模與仿真。

    李傳旭(1997—),男,工程師,碩士,主要研究方向?yàn)轱w行器制導(dǎo)與控制、系統(tǒng)仿真。

    猜你喜歡
    滑??刂?/a>參數(shù)優(yōu)化穩(wěn)定性
    非線性中立型變延遲微分方程的長(zhǎng)時(shí)間穩(wěn)定性
    基于干擾觀測(cè)器的PID滑模變結(jié)構(gòu)控制
    基于神經(jīng)網(wǎng)絡(luò)的動(dòng)力電池組焊接參數(shù)優(yōu)化研究
    基于多算法的ROV遠(yuǎn)程協(xié)同控制器設(shè)計(jì)與實(shí)現(xiàn)
    半動(dòng)力系統(tǒng)中閉集的穩(wěn)定性和極限集映射的連續(xù)性
    研究LTE與WCDMA系統(tǒng)間小區(qū)互操作與參數(shù)優(yōu)化
    基于磁流變技術(shù)的汽車(chē)發(fā)動(dòng)機(jī)隔振系統(tǒng)的參數(shù)優(yōu)化
    科技視界(2016年23期)2016-11-04 08:17:36
    改進(jìn)的動(dòng)態(tài)面船舶航跡跟蹤控制研究
    上向進(jìn)路式尾砂膠結(jié)充填采礦法采場(chǎng)結(jié)構(gòu)參數(shù)優(yōu)化研究
    中國(guó)科技博覽(2016年5期)2016-04-23 05:45:45
    99国产精品一区二区三区| 国产精品,欧美在线| 国产野战对白在线观看| 亚洲国产中文字幕在线视频| 大型黄色视频在线免费观看| 国产成人aa在线观看| 美女午夜性视频免费| 国产探花在线观看一区二区| 久久九九热精品免费| 久久久久久久久久黄片| 黄色 视频免费看| 脱女人内裤的视频| 特大巨黑吊av在线直播| 在线观看日韩欧美| 午夜福利在线观看吧| 99国产综合亚洲精品| 真人一进一出gif抽搐免费| 亚洲男人的天堂狠狠| 韩国av一区二区三区四区| 999久久久精品免费观看国产| 欧美在线黄色| 精品免费久久久久久久清纯| 狂野欧美激情性xxxx| 一级毛片女人18水好多| 不卡一级毛片| 中文字幕人妻丝袜一区二区| 给我免费播放毛片高清在线观看| 欧美在线一区亚洲| 丰满的人妻完整版| 欧美绝顶高潮抽搐喷水| 中出人妻视频一区二区| 欧美在线一区亚洲| 欧美色视频一区免费| 亚洲第一电影网av| 亚洲欧美精品综合一区二区三区| 久久99热这里只有精品18| 一个人免费在线观看电影 | 国内精品一区二区在线观看| 亚洲国产精品999在线| 久久精品综合一区二区三区| 午夜精品在线福利| 亚洲欧美日韩高清在线视频| 身体一侧抽搐| 日韩三级视频一区二区三区| 国产欧美日韩一区二区精品| 精品国内亚洲2022精品成人| 一进一出抽搐gif免费好疼| 国产私拍福利视频在线观看| www.www免费av| 人妻久久中文字幕网| 国产熟女xx| 大型黄色视频在线免费观看| 国产精品一区二区三区四区久久| 好男人在线观看高清免费视频| 欧美大码av| 国产精品99久久99久久久不卡| 欧美av亚洲av综合av国产av| 搡老妇女老女人老熟妇| 欧美极品一区二区三区四区| 曰老女人黄片| 九九久久精品国产亚洲av麻豆 | 国产激情偷乱视频一区二区| 亚洲专区中文字幕在线| 一个人免费在线观看的高清视频| 日本撒尿小便嘘嘘汇集6| 成人欧美大片| 欧美+亚洲+日韩+国产| 麻豆国产97在线/欧美| 午夜视频精品福利| 搡老妇女老女人老熟妇| 免费在线观看成人毛片| netflix在线观看网站| av在线蜜桃| 久久这里只有精品中国| 欧美黑人巨大hd| 在线免费观看不下载黄p国产 | 成人特级av手机在线观看| 午夜精品在线福利| 啦啦啦观看免费观看视频高清| 在线a可以看的网站| 亚洲自拍偷在线| 精品久久久久久久毛片微露脸| 久久久国产成人免费| 久久国产精品人妻蜜桃| 国内精品美女久久久久久| 在线永久观看黄色视频| 久久久久国内视频| 国产高清三级在线| 国内精品美女久久久久久| 丰满人妻一区二区三区视频av | 成人性生交大片免费视频hd| 精品久久久久久久人妻蜜臀av| 不卡一级毛片| 亚洲av日韩精品久久久久久密| 亚洲欧美日韩高清专用| 99久久精品一区二区三区| av在线蜜桃| 欧美黑人欧美精品刺激| 可以在线观看的亚洲视频| 精品国产乱子伦一区二区三区| 中文亚洲av片在线观看爽| 国产免费av片在线观看野外av| 国产成人精品久久二区二区91| 国产日本99.免费观看| 亚洲国产日韩欧美精品在线观看 | 久久久成人免费电影| 亚洲国产欧美一区二区综合| 精品一区二区三区av网在线观看| 久久精品人妻少妇| 亚洲欧美日韩卡通动漫| 成人三级做爰电影| 成年女人看的毛片在线观看| 噜噜噜噜噜久久久久久91| 日本黄色片子视频| 在线观看免费视频日本深夜| 久久久久久久午夜电影| 精品99又大又爽又粗少妇毛片 | 国产美女午夜福利| 亚洲国产精品合色在线| 亚洲18禁久久av| 日韩有码中文字幕| 久久久久国产精品人妻aⅴ院| 一区福利在线观看| 麻豆成人午夜福利视频| 不卡av一区二区三区| 日本免费一区二区三区高清不卡| av在线天堂中文字幕| 日本一本二区三区精品| 久久久久久九九精品二区国产| 一个人免费在线观看的高清视频| 欧美日韩亚洲国产一区二区在线观看| 国产午夜福利久久久久久| 麻豆国产97在线/欧美| 久久久久性生活片| 黄色成人免费大全| 亚洲av日韩精品久久久久久密| 成年免费大片在线观看| 国产亚洲精品久久久久久毛片| 男女床上黄色一级片免费看| 久久久国产成人精品二区| 国产乱人伦免费视频| 日韩欧美三级三区| 国产精品亚洲av一区麻豆| 国产乱人伦免费视频| 亚洲av美国av| 美女高潮的动态| 成年版毛片免费区| 亚洲精品粉嫩美女一区| 人妻久久中文字幕网| 亚洲成人免费电影在线观看| 久久久精品欧美日韩精品| 亚洲第一欧美日韩一区二区三区| 一区二区三区激情视频| 国内毛片毛片毛片毛片毛片| 国产精品一区二区免费欧美| 久久久国产欧美日韩av| 亚洲五月婷婷丁香| 搡老熟女国产l中国老女人| 久久精品夜夜夜夜夜久久蜜豆| 日本一二三区视频观看| 免费看光身美女| 国产精品久久久久久精品电影| 毛片女人毛片| 老司机深夜福利视频在线观看| 免费av毛片视频| 国产淫片久久久久久久久 | 99久久精品国产亚洲精品| 国产精品国产高清国产av| 亚洲,欧美精品.| 很黄的视频免费| 成年版毛片免费区| 亚洲av电影在线进入| 色综合欧美亚洲国产小说| 女同久久另类99精品国产91| 久久中文字幕一级| 久久国产精品人妻蜜桃| 午夜福利18| 国产精品久久久久久人妻精品电影| 亚洲av成人精品一区久久| 日韩精品青青久久久久久| 亚洲av熟女| 久久久精品大字幕| 一级a爱片免费观看的视频| 十八禁人妻一区二区| 一区二区三区国产精品乱码| 久久精品国产清高在天天线| 女人被狂操c到高潮| 黑人巨大精品欧美一区二区mp4| 国产精品九九99| 成熟少妇高潮喷水视频| 欧美+亚洲+日韩+国产| 91久久精品国产一区二区成人 | 超碰成人久久| 哪里可以看免费的av片| 欧洲精品卡2卡3卡4卡5卡区| 欧美又色又爽又黄视频| 青草久久国产| 网址你懂的国产日韩在线| 天堂√8在线中文| 少妇的丰满在线观看| 亚洲熟女毛片儿| 最近视频中文字幕2019在线8| 亚洲avbb在线观看| 又粗又爽又猛毛片免费看| 国产精品久久久久久人妻精品电影| 精品国内亚洲2022精品成人| 色综合亚洲欧美另类图片| 可以在线观看的亚洲视频| 国产精品国产高清国产av| 国产av不卡久久| 婷婷丁香在线五月| 神马国产精品三级电影在线观看| 禁无遮挡网站| www日本在线高清视频| 长腿黑丝高跟| 一级作爱视频免费观看| 欧美最黄视频在线播放免费| 久久久久久人人人人人| 亚洲精品国产精品久久久不卡| 中文字幕精品亚洲无线码一区| 少妇熟女aⅴ在线视频| 精品乱码久久久久久99久播| 美女扒开内裤让男人捅视频| 亚洲真实伦在线观看| 久久久久国内视频| 日韩精品中文字幕看吧| 露出奶头的视频| 亚洲第一电影网av| www日本黄色视频网| 最近在线观看免费完整版| 国内精品一区二区在线观看| 嫩草影院精品99| 亚洲欧洲精品一区二区精品久久久| 日本五十路高清| a在线观看视频网站| 国产午夜精品论理片| 欧美黑人欧美精品刺激| 黄色丝袜av网址大全| 国产精品久久久久久亚洲av鲁大| 麻豆国产av国片精品| 国产精品国产高清国产av| 女生性感内裤真人,穿戴方法视频| 法律面前人人平等表现在哪些方面| 欧美+亚洲+日韩+国产| www.www免费av| 久久精品亚洲精品国产色婷小说| 久久久久久久午夜电影| 精品国产乱子伦一区二区三区| 99re在线观看精品视频| 美女高潮的动态| 亚洲午夜精品一区,二区,三区| 天堂动漫精品| 中文字幕最新亚洲高清| 亚洲片人在线观看| 欧美成人免费av一区二区三区| 欧美一区二区精品小视频在线| 日本 欧美在线| 亚洲成av人片在线播放无| 日韩欧美三级三区| 色av中文字幕| 91av网一区二区| 国产精品美女特级片免费视频播放器 | 亚洲国产高清在线一区二区三| 免费观看人在逋| 日日夜夜操网爽| 久久久成人免费电影| 午夜激情欧美在线| 日本黄色视频三级网站网址| 麻豆av在线久日| 国产成+人综合+亚洲专区| 老汉色av国产亚洲站长工具| 真实男女啪啪啪动态图| 国产精品久久久人人做人人爽| 桃色一区二区三区在线观看| 国产精品久久久久久亚洲av鲁大| 日韩有码中文字幕| 悠悠久久av| 18美女黄网站色大片免费观看| 无限看片的www在线观看| 19禁男女啪啪无遮挡网站| 丰满人妻熟妇乱又伦精品不卡| 日韩成人在线观看一区二区三区| aaaaa片日本免费| 1000部很黄的大片| 毛片女人毛片| 精品国产亚洲在线| 99riav亚洲国产免费| 岛国在线观看网站| 两个人视频免费观看高清| 日韩欧美在线乱码| 成人特级av手机在线观看| 男女那种视频在线观看| 久久亚洲精品不卡| 亚洲avbb在线观看| 午夜a级毛片| 女人被狂操c到高潮| 好男人电影高清在线观看| 国产亚洲精品久久久com| 色老头精品视频在线观看| 久久中文字幕一级| 一本综合久久免费| 叶爱在线成人免费视频播放| 亚洲国产日韩欧美精品在线观看 | 黄色片一级片一级黄色片| 国产一区二区三区视频了| 色精品久久人妻99蜜桃| 曰老女人黄片| 亚洲黑人精品在线| 国产精品国产高清国产av| 精品国产乱码久久久久久男人| 小蜜桃在线观看免费完整版高清| 亚洲中文av在线| 久久久久久久久免费视频了| 97超视频在线观看视频| 精品久久久久久久久久久久久| 精华霜和精华液先用哪个| 在线看三级毛片| 男人舔女人的私密视频| 日韩有码中文字幕| 亚洲 国产 在线| 99在线视频只有这里精品首页| 亚洲成人免费电影在线观看| 免费观看的影片在线观看| 一个人看的www免费观看视频| 搡老妇女老女人老熟妇| 亚洲最大成人中文| 99国产精品一区二区蜜桃av| 欧美日韩国产亚洲二区| 我要搜黄色片| 亚洲自偷自拍图片 自拍| av在线天堂中文字幕| 在线国产一区二区在线| 亚洲av电影在线进入| 天堂av国产一区二区熟女人妻| 岛国视频午夜一区免费看| 中文在线观看免费www的网站| av欧美777| 成人鲁丝片一二三区免费| 国产视频一区二区在线看| 舔av片在线| 校园春色视频在线观看| 国产午夜福利久久久久久| 国产精品爽爽va在线观看网站| 成熟少妇高潮喷水视频| 夜夜看夜夜爽夜夜摸| 九九热线精品视视频播放| 国产视频内射| 日本成人三级电影网站| 亚洲精品粉嫩美女一区| 国产精品久久电影中文字幕| 1000部很黄的大片| 男女之事视频高清在线观看| 最近最新中文字幕大全电影3| 好男人电影高清在线观看| 欧美在线一区亚洲| 国产av麻豆久久久久久久| 国产精品香港三级国产av潘金莲| 最近在线观看免费完整版| 亚洲欧美日韩东京热| 深夜精品福利| 最新在线观看一区二区三区| 夜夜看夜夜爽夜夜摸| 亚洲七黄色美女视频| 国产激情偷乱视频一区二区| 国产高清视频在线观看网站| 看免费av毛片| 成人特级av手机在线观看| 成人精品一区二区免费| 天天添夜夜摸| 一本综合久久免费| 亚洲欧美日韩卡通动漫| 亚洲精品一卡2卡三卡4卡5卡| 亚洲国产欧美一区二区综合| 亚洲中文字幕日韩| 国产淫片久久久久久久久 | 国产av麻豆久久久久久久| 变态另类丝袜制服| 亚洲av成人一区二区三| 香蕉丝袜av| 亚洲自偷自拍图片 自拍| 搡老妇女老女人老熟妇| 欧美精品啪啪一区二区三区| a级毛片在线看网站| 变态另类成人亚洲欧美熟女| 男插女下体视频免费在线播放| 日韩中文字幕欧美一区二区| 最近在线观看免费完整版| 日韩欧美 国产精品| 操出白浆在线播放| 精品人妻1区二区| 狂野欧美白嫩少妇大欣赏| 嫩草影视91久久| 老司机午夜十八禁免费视频| 嫩草影视91久久| 香蕉国产在线看| 亚洲色图av天堂| 亚洲狠狠婷婷综合久久图片| 制服人妻中文乱码| 亚洲自拍偷在线| 精品欧美国产一区二区三| 欧美乱色亚洲激情| 丰满的人妻完整版| 欧美黄色淫秽网站| 色哟哟哟哟哟哟| 夜夜爽天天搞| 99精品在免费线老司机午夜| 久久婷婷人人爽人人干人人爱| 久久久久性生活片| 精品一区二区三区视频在线观看免费| 免费看美女性在线毛片视频| 亚洲欧美日韩高清在线视频| 热99在线观看视频| 亚洲成a人片在线一区二区| 亚洲电影在线观看av| 亚洲精品中文字幕一二三四区| 99久久精品一区二区三区| 搡老岳熟女国产| 久99久视频精品免费| 国产精品1区2区在线观看.| 欧美性猛交╳xxx乱大交人| 99精品在免费线老司机午夜| 国产免费av片在线观看野外av| 国产亚洲av嫩草精品影院| 欧美高清成人免费视频www| 黄色女人牲交| 18禁裸乳无遮挡免费网站照片| 色综合欧美亚洲国产小说| 嫩草影视91久久| 久久九九热精品免费| 国语自产精品视频在线第100页| 99国产综合亚洲精品| 窝窝影院91人妻| 美女黄网站色视频| 精品午夜福利视频在线观看一区| 一区二区三区高清视频在线| 又大又爽又粗| www日本在线高清视频| 一个人免费在线观看的高清视频| 夜夜躁狠狠躁天天躁| 久久精品aⅴ一区二区三区四区| 国产精品美女特级片免费视频播放器 | 亚洲精品色激情综合| 成人鲁丝片一二三区免费| www.自偷自拍.com| 性欧美人与动物交配| 色噜噜av男人的天堂激情| 级片在线观看| 成年女人永久免费观看视频| 特级一级黄色大片| 丝袜人妻中文字幕| 大型黄色视频在线免费观看| 99久久综合精品五月天人人| 美女黄网站色视频| 中文字幕av在线有码专区| 国内毛片毛片毛片毛片毛片| 黄色成人免费大全| 在线观看日韩欧美| 成人三级黄色视频| 中文字幕最新亚洲高清| 亚洲av日韩精品久久久久久密| 国产免费男女视频| www.www免费av| 三级国产精品欧美在线观看 | 男人舔女人下体高潮全视频| 99热这里只有是精品50| 国产一区二区三区在线臀色熟女| 国产精品国产高清国产av| 桃红色精品国产亚洲av| 国产精品久久久久久亚洲av鲁大| 亚洲avbb在线观看| 天堂√8在线中文| 国产熟女xx| 日本 欧美在线| 黄色 视频免费看| 色视频www国产| 草草在线视频免费看| 免费观看精品视频网站| 日韩欧美免费精品| 99视频精品全部免费 在线 | 国产精品av视频在线免费观看| 亚洲在线自拍视频| 国产精品女同一区二区软件 | 亚洲国产色片| 亚洲欧美日韩无卡精品| 国产精品久久电影中文字幕| 亚洲精品美女久久av网站| 一本综合久久免费| 高清毛片免费观看视频网站| 欧美又色又爽又黄视频| 日本在线视频免费播放| 亚洲中文av在线| 欧美精品啪啪一区二区三区| 成年人黄色毛片网站| 最好的美女福利视频网| 午夜久久久久精精品| 亚洲中文av在线| 久久久久国产精品人妻aⅴ院| 国产亚洲av嫩草精品影院| 精品人妻1区二区| 成人18禁在线播放| 久久久久久大精品| 一个人看的www免费观看视频| 国产淫片久久久久久久久 | 757午夜福利合集在线观看| 国产亚洲av嫩草精品影院| 国产精品综合久久久久久久免费| 亚洲人成电影免费在线| 国产爱豆传媒在线观看| 在线观看免费午夜福利视频| 精品一区二区三区视频在线观看免费| 久久久久国内视频| 在线观看免费视频日本深夜| 欧美xxxx黑人xx丫x性爽| 国产成+人综合+亚洲专区| 叶爱在线成人免费视频播放| 波多野结衣巨乳人妻| 琪琪午夜伦伦电影理论片6080| 国产精品电影一区二区三区| 哪里可以看免费的av片| 人妻丰满熟妇av一区二区三区| 十八禁网站免费在线| 一级毛片女人18水好多| 长腿黑丝高跟| 国产av麻豆久久久久久久| 国产高清视频在线观看网站| 免费在线观看日本一区| 亚洲av成人av| 怎么达到女性高潮| 床上黄色一级片| 国产精品野战在线观看| 亚洲18禁久久av| 精品久久久久久久久久免费视频| 嫩草影院精品99| 在线观看免费视频日本深夜| 亚洲成a人片在线一区二区| 1024香蕉在线观看| 搡老熟女国产l中国老女人| 18禁裸乳无遮挡免费网站照片| 婷婷精品国产亚洲av在线| 国产av一区在线观看免费| 噜噜噜噜噜久久久久久91| 精品不卡国产一区二区三区| 在线看三级毛片| 美女黄网站色视频| 亚洲自偷自拍图片 自拍| 男人和女人高潮做爰伦理| 日韩欧美精品v在线| 欧美一级a爱片免费观看看| 国产黄色小视频在线观看| 国产一级毛片七仙女欲春2| 少妇的丰满在线观看| 一二三四社区在线视频社区8| 欧美精品啪啪一区二区三区| 亚洲人成电影免费在线| 婷婷六月久久综合丁香| 91在线观看av| 日本一二三区视频观看| 一级a爱片免费观看的视频| 99热精品在线国产| 亚洲黑人精品在线| 男人舔奶头视频| av黄色大香蕉| 少妇裸体淫交视频免费看高清| 亚洲 欧美 日韩 在线 免费| 亚洲精品456在线播放app | 99国产精品一区二区三区| 国产淫片久久久久久久久 | 久久国产精品影院| 桃色一区二区三区在线观看| 五月玫瑰六月丁香| 757午夜福利合集在线观看| 可以在线观看的亚洲视频| 黄色女人牲交| 一本综合久久免费| 99在线视频只有这里精品首页| 男女下面进入的视频免费午夜| 一级毛片高清免费大全| 俺也久久电影网| 97人妻精品一区二区三区麻豆| 亚洲 欧美一区二区三区| a级毛片在线看网站| 最好的美女福利视频网| 人妻丰满熟妇av一区二区三区| 色播亚洲综合网| 国产精品影院久久| 亚洲在线自拍视频| 免费搜索国产男女视频| 女人高潮潮喷娇喘18禁视频| 久久久精品大字幕| 亚洲黑人精品在线| 亚洲真实伦在线观看| 夜夜躁狠狠躁天天躁| 一个人观看的视频www高清免费观看 | 久久久久久国产a免费观看| 亚洲专区国产一区二区| 日本黄色片子视频| 变态另类丝袜制服| 日本成人三级电影网站| 国产单亲对白刺激| 色综合站精品国产| 亚洲电影在线观看av| 精品久久久久久成人av| 俄罗斯特黄特色一大片| 99久久久亚洲精品蜜臀av| 日本在线视频免费播放| 观看美女的网站| 精品电影一区二区在线| 级片在线观看| 最近最新中文字幕大全免费视频| 日韩精品青青久久久久久| 99热这里只有是精品50| 母亲3免费完整高清在线观看| 村上凉子中文字幕在线| 日日摸夜夜添夜夜添小说| 黄片大片在线免费观看|