-
在火星探测过程中,当飞行器从环绕停泊轨道至进入大气的过程之间,需执行连续推力机动实现离轨过程,因此,整个离轨过程从问题本质上属于多重弧段(即包括多段推力弧和滑行弧)有限推力轨道机动问题。本文考虑到实际任务时间限定和有利实际工程应用等客观因素,以一次“推–滑”两弧段转移形式为离轨过程飞行器的机动模式,并在此基础上给出最优制导方法。
尽管离轨机动本质属于有限推力机动过程,但部分约束形式等方面又具有其独特的一面。由于一般进入大气前的停泊轨道高度较低,因此,为保证制导算法在进行离轨轨迹求解时的实时性要求,其动力学方程的形式、状态变量以及协状态变量的求解等均会在一定的精度要求和假设前提下进行简化,所以动力学模型以及问题的描述形式也不同于一般的有限推力轨道机动问题。
-
通过火星惯性系下的动力学模型来描述离轨机动问题。其对应的无量纲动力学方程为
$$ \left\{ {\begin{array}{*{20}{l}} {{{\dot r}} = {{V}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;}\\ {{{\dot V}} = - \dfrac{1}{{{r^3}}}{{r}} + \dfrac{ T}{{m{{ g}_0}}}{{\upsilon }}}\\ {\dot m = - \dfrac{{\sqrt R { T}}}{{{I_{\rm sp}}{ g}_0^{1.5}}}\;\;\;\;\;\;\;\;} \end{array}} \right. $$ (1) 其中:r为火心惯性系下的无量纲位置矢量;V为火心惯性系下的无量纲速度矢量;m为有量纲的飞行器质量;T为飞行器发动机推力,Isp为飞行器发动机比冲;
${{ g}_0}$ 为火星表面引力加速度;${{\upsilon}}$ 为推力方向在火心惯性系下的单位矢量。在无量纲化过程中,距离的无量纲参数为火星半径R;速度的无量纲参数为$\sqrt {R{{ g}_0}} $ ;时间的无量纲参数为$\sqrt {R/{{ g}_0}} $ 。由于大气进入前的停泊轨道高度较低,且位置矢径大小对应的高度介于停泊轨道高度与大气边缘高度之间,因此,整个离轨过程位置矢径大小变化微小。基于该特性,在每次迭代制导的过程中,位置矢径可近似为[1]
$${ r} \approx {{ r}_0}$$ (2) 其中:r0为每次制导迭代初始时刻的无量纲位置矢径大小。因此,式(1)中动力学方程中速度矢量的微分项可重新表示为[1]
$${{\dot V}} = - {\varpi ^2}{{r}} + \dfrac{ T}{{m{{ g}_0}}}{{\upsilon }}$$ (3) 其中:为方便后续公式表示,取
$\varpi = 1/r_0^{1.5}$ 。 -
基于式(1)和式(3)的动力学模型,在给出最优性条件之前,首先给出协状态方程的表达式。假设速度矢量、位置矢量以及质量对应的协状态量分别为pr,pV和pm,基于最优控制理论[14],系统的哈密顿函数为
$$H = {{p}}_r^{\rm{T}}{{V}} - {\varpi ^2}{{p}}_V^{\rm{T}}{{r}} + \dfrac{ T}{{m{{ g}_0}}}{{p}}_V^{\rm{T}}{{\upsilon }}-\dfrac{{\sqrt R { T}}}{{{I_{\rm sp}}{ g}_0^{1.5}}}{p_m}$$ (4) 其中,协状态量pr和pV对应的微分方程为
$$ \left\{ {\begin{array}{*{20}{c}} {{{{{\dot p}}}_r} = - \partial H/\partial {{r}} = {\varpi ^2}{{{p}}_V}}\\ {{{{{\dot p}}}_V} = - \partial H/\partial {{V}} = - {{{p}}_r}\;} \end{array}} \right. $$ (5) 在离轨过程中,由于“推–滑”离轨机动模式下,飞行器发动机机动模式固定,所以此处控制变量只有“推”弧段的推力方向单位矢量。根据极大值原理,最优离轨过程对应的发动机推力方向单位矢量
${{\upsilon}}$ 为$${{{\upsilon}}^ * } = {{{p}}_V}/\left\| {{{{p}}_V}} \right\|$$ (6) 上述式(6)中给出的推力方向单位矢量
${{{\upsilon}}^ * }$ 是最优离轨机动过程的必要条件。 -
在如式(2)的假设前提下,离轨机动的动力学方程中部分项被线性化,而式(5)中的协状态微分方程直接具有线性微分方程的基本形式,这种特性给微分方程的解析求解提供了一定的便利。因此,本文采用一种将协态变量和状态变量变换的方法,求解得到状态方程和协态方程的解析解[1-2,12]。事实上,这种解析求解方法目前已经成熟应用于离轨制导和入轨制导等问题中。下面将给出该方法求解状态和协态变量解析解的具体过程。
首先,对协态变量进行变化,即引入一个新的量为
$${{\lambda}}\left( t \right) = {\left[ {{{p}}_V^{\rm{T}}\left( t \right), - {{p}}_r^{\rm{T}}\left( t \right)/\varpi } \right]^{\rm{T}}}$$ (7) 对于新变量,此处省略推导过程,直接给出对应的从式(5)变形得到的微分方程的解析解,即
$$ \begin{array}{l} {{\lambda}}\left( t \right) = {{\varPhi}} \left( {t - {t_0}} \right){{\lambda}}\left( {{t_0}} \right) \\ \;\;\;\;\;\; = \left[ {\begin{array}{*{20}{c}} {\cos \left[ {\varpi \left( {t - {t_0}} \right)} \right]{{{I}}_3}}&{\sin \left[ {\varpi \left( {t - {t_0}} \right)} \right]{{{I}}_3}} \\ { - \sin \left[ {\varpi \left( {t - {t_0}} \right)} \right]{{{I}}_3}}&{\cos \left[ {\varpi \left( {t - {t_0}} \right)} \right]{{{I}}_3}} \end{array}} \right]{{\lambda}}\left( {{t_0}} \right)\; \\ \end{array} $$ (8) 其中:t和t0分别是当前时间和初始时间,I3为3阶单位矩阵,
${{\lambda}}\left( {{t_0}} \right)$ 为初始时刻对应的值。其次,对于状态变量,此处同样需进行变换从而使得状态变量的微分方程可以解析得到,引入一个新的量为
$$ {{\varTheta }}\left( t \right) = {\left[ {{{{r}}^{\rm{T}}}\left( t \right),{{{V}}^{\rm{T}}}\left( t \right)/\varpi } \right]^{\rm{T}}} $$ (9) 根据参考文献[12]中的方法,在给出如式(9)的变换之后,从状态变量转换得到的新变量对应的微分方程具有如式(10)的解析解形式,即
$$ {{\varTheta }}\left( t \right) = {{\varPhi }}\left( {t - {t_0}} \right){{\varTheta }}\left( {{t_0}} \right) + {{\varXi }}\left( t \right){{\varGamma }}\left( {t,{t_0}} \right) $$ (10) 其中:
${{\varTheta }}\left( {{t_0}} \right)$ 为初始时刻对应的值,${{\varXi }}\left( t \right)$ 为$${{\varXi }}\left( t \right) = \left[ {\begin{array}{*{20}{c}} {\sin \left( {\varpi t} \right){{{I}}_3}}&{ - \cos \left( {\varpi t} \right){{{I}}_3}} \\ {\cos \left( {\varpi t} \right){{{I}}_3}}&{\sin \left( {\varpi t} \right){{{I}}_3}} \end{array}} \right]/\varpi $$ (11) 此外,式(10)中的向量
${{\varGamma }}\left( {t,{t_0}} \right)$ 的表达式为两个子向量的组合,其具体表达式为$$ {{\varGamma }}\left( {t,{t_0}} \right) = {\left[ {{{\varGamma }}_1^{\rm{T}}\left( {t,{t_0}} \right),{{\varGamma }}_2^{\rm{T}}\left( {t,{t_0}} \right)} \right]^{\rm{T}}} $$ (12) 其中,子向量的表达式分别为
$${{{\varGamma }}_1}\left( {t,{t_0}} \right) = \dfrac{ T}{{{{ g}_0}}}\int_{{t_0}}^t {\dfrac{{{{{\upsilon}}^ * }\left( \chi \right)}}{{m\left( \chi \right)}}} \cos \left( {\varpi \chi } \right){\rm{d}}\chi $$ (13) $${{{\varGamma }}_2}\left( {t,{t_0}} \right) = \dfrac{ T}{{{{ g}_0}}}\int_{{t_0}}^t {\dfrac{{{{{\upsilon}}^ * }\left( \chi \right)}}{{m\left( \chi \right)}}} \sin \left( {\varpi \chi } \right){\rm{d}}\chi $$ (14) 对于式(13)和式(14)中子向量的值,在每一时刻均采用数值积分求解,一般情况下,只要数据节点之间的间隔足够小,其积分精度能够得到保证。为了保证求解速度的要求,两个时间节点间隔之间只需进行一次积分即可,数值积分方法诸如4阶龙格库塔法以及米勒法等均能够以较高精度实现对上述积分函数的求解。具体求解公式此处不再赘述。
-
1.3小节给出了状态变量和协态变量的解析计算方法,对于离轨制导问题,从问题本质上属于典型的两点边值问题,因此,下面需要分别确定两点边值问题对应的变量及终端约束。对于离轨机动过程,一般初始状态是确定的,初始协态变量和总的离轨时间是可变的,且终端状态受约束。在离轨制导过程中,制导参数包括初始协态变量和总离轨时间,但是由于终端状态变量受约束,因此直接根据变分法得到的制导终端约束会包含拉格朗日乘子,因此下面将推导给出只包含状态变量和协态变量的制导终端约束。
对于火星离轨机动过程,一般包括如下3个终端状态约束,即
$$\dfrac{{{{r}}_{\rm f}^{\rm T}{{{r}}_{\rm f}}}}{2} - \dfrac{{r_{{\rm{EI}}}^2}}{2} = 0$$ (15) $$\dfrac{{{{V}}_{\rm f}^{\rm T}{{{V}}_{\rm f}}}}{2} - \dfrac{{V_{{\rm{EI}}}^2}}{2} = 0$$ (16) $${{V}}_{\rm f}^{\rm T}{{{r}}_{\rm f}} - {V_{{\rm{EI}}}}{r_{{\rm{EI}}}}\sin {\gamma _{{\rm{EI}}}} = 0$$ (17) 其中:rf和Vf分别是终端位置矢量和速度矢量,而rEI、VEI和
${\gamma _{{\rm{EI}}}}$ 分别为基于上述3个终端状态约束。对于最优离轨过程,其横截条件为$$\left\{ {\begin{array}{*{20}{c}} {{{{p}}_r}\left( {{t_{\rm f}}} \right) = {\mu _1}{{{r}}_{\rm f}} + {\mu _3}{{{V}}_{\rm f}}} \\ {{{{p}}_V}\left( {{t_{\rm f}}} \right) = {\mu _3}{{{r}}_{\rm f}} + {\mu _2}{{{V}}_{\rm f}}\;} \end{array}} \right.$$ (18) 其中:tf是终端时间;
${\mu _1}$ 、${\mu _2}$ 和${\mu _3}$ 是拉格朗日乘子。为了消去式(18)中的拉格朗日乘子,下面将结合式(15)~(17)与式(18)给出只包含状态和协状态的制导终端约束。首先对式(18)等式左右分别转置之后右乘位置矢量rf和速度矢量Vf,从而引入式(15)~(17)的等式关系,也即
$$\left\{ {\begin{array}{*{20}{c}} {{{p}}_r}^{\rm{T}}\left( {{t_{\rm f}}} \right){{{r}}_{\rm f}} = {\mu _1}r_{{\rm{EI}}}^2 + {\mu _3}{V_{{\rm{EI}}}}{r_{{\rm{EI}}}}\sin {\gamma _{{\rm{EI}}}} \\ {{{p}}_r}^{\rm{T}}\left( {{t_{\rm f}}} \right){{{V}}_{\rm f}} = {\mu _1}{V_{{\rm{EI}}}}{r_{{\rm{EI}}}}\sin {\gamma _{{\rm{EI}}}} + {\mu _3}V_{{\rm{EI}}}^2 \\ {{{{p}}_V}^{\rm{T}}\left( {{t_{\rm f}}} \right){{{r}}_{\rm f}} = {\mu _2}{V_{{\rm{EI}}}}{r_{{\rm{EI}}}}\sin {\gamma _{{\rm{EI}}}} + {\mu _3}r_{{\rm{EI}}}^2} \end{array}} \right.$$ (19) 通过式(19)中的3个公式,可以直接求解得到3个拉格朗日乘子
${\mu _1}$ 、${\mu _2}$ 和${\mu _3}$ ,然后将其带入式(18),可得到不含拉格朗日乘子的横截条件,其最终终端约束方程为$$ \left[ {\begin{array}{*{20}{c}} {{{{p}}_r}^{\rm{T}}\left( {{t_{\rm f}}} \right)}\\ {{{{p}}_V}^{\rm{T}}\left( {{t_{\rm f}}} \right)} \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{{{r}}_{\rm f}}}&{{{\bf{0}}_{3 \times 1}}}&{{{{V}}_{\rm f}}}\\ {{{\bf{0}}_{3 \times 1}}}&{{{{V}}_{\rm f}}}&{{{{r}}_{\rm f}}} \end{array}} \end{array}} \right]{{{A}}^{ - 1}}\left[ {\begin{array}{*{20}{c}} {{{{p}}_{\rm f}}^{\rm{T}}\left( {{t_{\rm f}}} \right){{{r}}_{\rm f}}}\\ {{{{p}}_r}^{\rm{T}}\left( {{t_{\rm f}}} \right){{{V}}_{\rm f}}}\\ {{{{p}}_V}^{\rm{T}}\left( {{t_{\rm f}}} \right){{{r}}_{\rm f}}} \end{array}} \right] = 0 $$ (20) 其中,中间矩阵A为
$${{A}} = \left[ {\begin{array}{*{20}{c}} {r_{{\rm{EI}}}^2}&0&{{V_{{\rm{EI}}}}{r_{{\rm{EI}}}}\sin {\gamma _{{\rm{EI}}}}} \\ {{V_{{\rm{EI}}}}{r_{{\rm{EI}}}}\sin {\gamma _{{\rm{EI}}}}}&0&{V_{{\rm{EI}}}^2} \\ 0&{{V_{{\rm{EI}}}}{r_{{\rm{EI}}}}\sin {\gamma _{{\rm{EI}}}}}&{r_{{\rm{EI}}}^2} \end{array}} \right]$$ (21) 式(20)即为最优“推–滑”离轨制导对应的由状态变量和协态变量构成的终端约束。此外,由于终端时间自由,因此,需要额外增加终端哈密顿函数为0的约束,即
$$H\left( {{t_f}} \right) = 0$$ (22) 至此,制导过程对应的变量和约束方程已给出。
在火星离轨制导过程中,为保证收敛性和求解速度每个制导周期内问题的求解均采用莱文贝格–马夸特(Levenberg-Marquardt)方法[15]。一般情况下,制导变量初值选取的优劣对首次制导方程求解的收敛速度有明显影响,因此,此处基于实际数值仿真结果,给出较为一般性得制导变量初值,即协态变量初值
${{{p}}_r}\left( {{t_0}} \right) = $ $ {\left[ {0,\;0,\;0} \right]^{\rm{T}}}$ ,${{{p}}_V}\left( {{t_0}} \right) = {{{V}}_0}/\left\| {{{{V}}_0}} \right\|$ 。此外,离轨总时间tf的初值选取对收敛速度影响不大,一般选取为初始轨道周期的一般较为合理。下文将给出离轨制导迭代流程。 -
为验证上述最优“推–滑”离轨制导方法的有效性和鲁棒性,本节以火星探测任务为背景,给出相关的仿真结果。首先,给定火星探测器在进入大气前所在的停泊轨道根数,其数值通过转换即可得到离轨机动的初始状态变量,如表1所示。
表 1离轨制导的初始轨道根数
Table 1.Initial orbit elements of deorbit guidance
根数 数值 半长轴/km 3 746 偏心率 0 轨道倾角/(°) 50 升交点赤/(°) 0 近心点幅角/(°) 0 真近点角/(°) 90 此外,给定火星半径R= 3 396 km,因此,初始轨道实际上为高度为350 km的圆轨道。在离轨机动过程中,飞行器的发动机推力T= 1 000 N,发动机比冲Isp= 311 s,飞行器初始质量假设为4 000 kg。由于整个离轨过程以飞行器到达大气进入边缘为终止位置,因此,此处给定大气边缘高度为128 km。此外,所有仿真均假定目标大气入口处速度约束为VEI= 3 516.2 m/s,飞行路径角
${\gamma _{{\rm{EI}}}}$ = –2.358°,而进入点的位置不进行约束。 -
基于上述离轨过程状态变量初值和发动机参数,本节给出了不考虑随机偏差影响的火星“推–滑”离轨机动最优解。图1给出了最优“推–滑”离轨机动对应的转移轨迹。从图中可以看出,整个过程机动弧段较短,最终实现的是类似于一个脉冲机动的效果,而滑行弧段较长。此外,从转移轨迹也能看出,在终端速度和飞行路径角约束下,整个机动几乎在同一轨道面内。
通过仿真可知,整个“推–滑”离轨机动最优解对应的离轨总时间为2 241.46 s,其中机动弧段总时间为295.31 s,而滑行弧段历时1 946.15 s。因此通过发动机参数推算,整个离轨机动所消耗的燃料为255.696 kg,对应的速度增量为76.29 m/s。对于机动弧段,图2给出了该阶段对应的发动机推力单位方向矢量在惯性系下的分布。为了更为直观展示推力方向与速度的夹角分布,图2还给出了推力方向角,此处推力方向角定义为推力方向与飞行器速度方向的夹角,从结果可以看出,整个离轨过程,推力几乎近似沿速度的反方向,也即减速制动从而使得飞行器降轨进入大气。此外,通过仿真计算,表2给出了最优解对应的初始协态变量的值。
表 2最优离轨过程对应的协态变量初值
Table 2.Initial costate value of optimal deorbit process
协态变量 数值 pr(t0) [1.256 76×10–5;–6.051 55×10–5;–7.211 96×10–5] pV(t0) [1.06×10–4;–1.448 1×10–7;–1.725 747×10–7] -
为了验证上述最优“推–滑”离轨制导方法的鲁棒性和精度,本节将在不确定性随机偏差存在的基础上进行蒙特卡罗仿真,并分析不确定性因素的存在对“推–滑”离轨制导过程的影响。在制导过程中,给定机动弧段制导周期为5 s,在每一次制导起始增加随机偏差,根据工程经验,此处给定位置偏差的3σ值为10 m,速度偏差的3σ值为0.1 m/s,蒙特卡罗打靶次数为1 000。
图3给出了蒙特卡罗仿真得到的最优离轨过程轨迹分布,从轨迹图可以看出尽管随机偏差存在,但整个轨迹散布范围很小。图4分别给出了终端状态和离轨总时间在随机偏差存在下的散布情况。从图中可以看出,进入点路径角的散布与标称值
${\gamma _{{\rm{EI}}}}$ 相差最大约为0.046°,而进入点速度的散布与标称值VEI相差最大约为0.5 m/s,因此,整个制导结果精度满足进入点状态偏差的精度的要求。从离轨总时间散布来看,其受随机偏差的影响较为明显,但其与最优解的值2 241.46 s相比最大偏差仅在8 s左右。此外,为反映机动过程推力方向的分布,图5给出了推力单位方向矢量分量与推力方向角在随机偏差影响下的散布情况,可以看出尽管随机偏差存在,但推力方向仍处于标称推力方向附近。
Research on Optimal Deorbit Guidance Method for Mars Exploration
-
摘要:火星的进入、下降与着陆(Entry,Descent,and Landing,EDL)是探测任务中的重要过程,而离轨制动是飞行器EDL过程前的关键机动过程,是保证飞行器大气进入过程顺利执行的决定性步骤。针对火星进入前的离轨机动问题,提出了一种“推–滑”最优离轨制导方法。该方法在状态方程局部线性化的基础上,给出了状态和协态变量的解析形式,并通过推导最优离轨机动的必要条件,给出了“推–滑”最优离轨制导律的求解方程。仿真结果表明:提出的制导策略在保证最优性和制导精度的前提下具有良好的鲁棒性,可为我国未来火星探测过程中的离轨机动提供理论依据和方法参考。Abstract:The EDL phases are of great importance in Mars exploration mission, while before which the deorbit braking is a key maneuver to ensure the smooth execution of the atmospheric entry process. Aiming at the problem of deorbit maneuver before Mars entry, a burn-coast optimal deorbit guidance method is presented. Based on the local linearization of the state equation, the analytical forms of state and costate variables are obtained. The solution equations of the optimal burn-coast deorbit guidance law are given through derivation of the necessary conditions for optimal deorbit maneuver. The simulation results show that the proposed guidance strategy has good robustness on the premise of guaranteeing optimality and guidance accuracy, providing theoretical reference for the future deorbit maneuver in Mars exploration.
-
Key words:
- Mars exploration/
- deorbit/
- optimal guidance/
- burn-coast
Highlights● A burn-coast optimal deorbit guidance method is presented for Mars deorbit maneuver. ● The deorbit guidance law are given based on local linearization and through derivation of the necessary conditions for optimal deorbit maneuver. ● Robustness and accuracy of the guidance law are proved by Monte-Carlo simulation, which shows the error of entry path angle is within 0.046°,while the error of entry velocity is within 0.5 m/s. -
表 1离轨制导的初始轨道根数
Table 1Initial orbit elements of deorbit guidance
根数 数值 半长轴/km 3 746 偏心率 0 轨道倾角/(°) 50 升交点赤/(°) 0 近心点幅角/(°) 0 真近点角/(°) 90 表 2最优离轨过程对应的协态变量初值
Table 2Initial costate value of optimal deorbit process
协态变量 数值 pr(t0) [1.256 76×10–5;–6.051 55×10–5;–7.211 96×10–5] pV(t0) [1.06×10–4;–1.448 1×10–7;–1.725 747×10–7] -
[1] BALDWIN M C,LU P. Optimal deorbit guidance[J]. Journal of Guidance,Control,and Dynamics,2012,35(1):93-103. [2] BALDWIN M C. Autonomous optimal deorbit guidance[D]. Iowa: Iowa State University, 2010. [3] GILLESPIE E M, SAHA K B. Crew return vehicle(CRV)deorbit opportunities[C]//Proc of AIAA Modeling and Simulation Technologies Conference. USA: AIAA, 2003. [4] PENZO P A. Orbit / deorbit analysis for a Mars rover and sample return mission[J]. Journal of Guidance,Control,and Dynamics,1990,27(4):425-432. [5] TUMINO G, ANGELINO E, LELEU F, et al. The IXV project the ESA re-entry system and technologies demonstrator paving the way to European autonomous space transportation and exploration endeavours: IAC-08-D2[R]. Glasgow, UK: European Space Agency, 2008. [6] GALMAN B A. Minimum energy deorbit[J]. Journal of Spacecraft and Rockets,1966,3(7):1030-1033.doi:10.2514/3.28592 [7] GALMAN B A. Minimum range-sensitivity deorbit[J]. Journal of Spacecraft and Rockets,1971,8(2):189-190.doi:10.2514/3.30243 [8] SUN F. Solution of the optimal two-impulse deorbit problem with the second impulse at the entry point[J]. Acta Astronautica,1982,9(1):1-13.doi:10.1016/0094-5765(82)90053-4 [9] 陈洪波,杨涤. 升力式再入飞行器离轨制动研究[J]. 飞行力学,2006,24(2):35-39.doi:10.3969/j.issn.1002-0853.2006.02.009CHEN H B,YANG D. De-orbit operations study of lifting reentry vehicle[J]. Flight Dynamics,2006,24(2):35-39.doi:10.3969/j.issn.1002-0853.2006.02.009 [10] BUSEMANN A,VINH N X. Optimum constrained disorbit by multiple impulses[J]. Journal of Optimization Theory and Applications,1968,2(1):40-64.doi:10.1007/BF00927162 [11] 汤国建. 载人飞船返回再入轨道和制导导航控制系统的研究与设计[D]. 长沙: 国防科技大学, 2000. [12] LU P,GRIFFIN B J. Rapid optimal multiburn ascent planning and guidance[J]. Journal of Guidance,Control,and Dynamics,2008,31(6):1656-1664. [13] ZHONG R,ZHU Z. Optimal current switch ing control of electrodynamic tethers for fast deorbit[J]. Journal of Guidance,Control,and Dynamics,2014,37:1501-1511. [14] BRYSON A E, HO Y C. Applied optimal control: optimization, estimation, and control[M]. USA: Halsted Press, 1975. [15] MORÉ J J. The Levenberg-marquardt algorithm:Implementation and theory[J]. Lecture Notes in Mathematics,1977,630:105-116.