昆虫在飞行过程中，其翅膀通常会不同程度背离拍动平面^[1]，形成如“8”字形^[1-2]、椭圆形、“U”形^[2-3]及这些形状组合而成的扑翼轨迹^[4-7]。扑翼轨迹的多样性与不确定性给人们带来了很大的困惑，为了弄清楚扑翼轨迹对其气动性能的影响，人们进行了一系列的研究^[8-12]，结果表明，不同扑翼轨迹对昆虫的飞行气动力和能耗不尽相同但影响不大，基本在10%左右。除了气动力和能耗以外，飞行动稳定性对昆虫的飞行至关重要。Jankauski等^[13]指出抬升角能使昆虫在稳定飞行状态下产生转矩，从而调整其飞行姿态，对昆虫的飞行控制和机动飞行产生影响。人们针对昆虫飞行动稳定性开展的大量研究工作，都将昆虫翅膀的运动简化为只有拍动运动和翻转运动，忽略了抬升角^[14-19]。
本文在此前工作^{[17, 19]}的基础上，计算了有抬升角时(即食蚜蝇以“U”形扑翼轨迹飞行时)食蚜蝇的纵向与横向稳定性导数，求解了稳定性矩阵的特征值，并将这些结果与无抬升角时的动稳定性结果^{[17, 19]}进行比较，研究了抬升角对食蚜蝇飞行动稳定性的影响。
1 计算模型及方法 1.1 运动方程与分析方法 为方便描述食蚜蝇的运动，分别引入固定在食蚜蝇身体上随其一起运动的随体坐标系O_bx_by_bz_b和在地面上的固定坐标系OXYZ，如图 1所示。O_b位于身体的质心，当食蚜蝇平衡飞行时，x_b和y_b轴在水平面内，x_b轴水平向前，y_b轴指向食蚜蝇身体右侧，z_b轴垂直向下。u、v和w分别为昆虫质心在x_b、y_b和z_b轴方向上的平动速度，而p、q和r分别为昆虫身体绕x_b、y_b和z_b轴的转动角速度，θ_b和γ分别为食蚜蝇身体上俯和向右滚转时x_b和y_b轴与水平面的夹角。

图 1 食蚜蝇的身体坐标系和地面上的固定坐标系及各变量示意图 Fig. 1 Schematic diagram of true hoverfly body coordinate system and fixed ground coordinate system with variables

图选项

与此前研究昆虫飞行动稳定性的方法^{[17, 19-21]}相同，本文仍采用平均模型和线化理论分析方法。X、Y和Z分别为气动力一个拍动周期内的平均值在x_b、y_b和z_b轴方向的分量，L、M和N分别为绕x_b、y_b和z_b轴的平均气动力矩。食蚜蝇在悬停飞行时，u、v、w、p、q和r及其导数均为零，且X、Y、L、M和N也均等于零，而Z=-mg，m为食蚜蝇的总质量，g为重力加速度。运动方程线化之后，便可将其解耦为纵向和横向2种运动方程^[20](见式(1)和式(2))，分别进行研究。

(1)

(2)

式中：“·”表示各运动量关于时间t的导数；A₁⁺为纵向稳定性矩阵；A₂⁺为横向稳定性矩阵。

(3)

(4)

式中：X_u⁺、Z_u⁺、M_u⁺、X_w⁺、Z_w⁺、M_w⁺、X_q⁺、Z_q⁺和M_q⁺为食蚜蝇在扰动运动时的纵向稳定性导数；Y_v⁺、L_v⁺、N_v⁺、Y_p⁺、L_p⁺、N_p⁺、Y_r⁺、L_r⁺和N_r⁺为食蚜蝇在扰动运动时的横向稳定性导数；I_x⁺、I_y⁺、I_z⁺和I_xz⁺分别为身体坐标系O_bx_by_bz_b下的转动惯量和惯性积；“+”的变量和参数均用参考速度U(U=2Φnr₂为平均拍动速度，Φ为拍翅幅值，n为拍翅频率，r₂为翅膀面积二阶矩折合半径)、参考长度c(c为翅膀平均弦长)和参考时间T(T=1/n为拍翅周期)进行了无量纲化处理，各无量纲表达式如下：δu⁺=δu/U，δv⁺=δv/U，δw⁺=δw/U，δp⁺=δpT，δq⁺=δqT，δr⁺=δrT，m⁺=m/(0.5ρUS_tT)，I_x⁺=I_x/(0.5ρU²S_tcT²)，I_y⁺=I_y/(0.5ρU²S_tcT²)，I_z⁺=I_z/(0.5ρU²S_tcT²)，I_xz⁺=I_xz/(0.5ρU²S_tcT²)，g⁺=gT/U，X⁺=X/(0.5ρU²S_t)，Y⁺=Y/(0.5ρU²S_t)，Z⁺=Z/(0.5ρU²S_t)，L⁺=L/(0.5ρU²S_tc)，M⁺=M/(0.5ρU²S_tc)，N⁺=N/(0.5ρU²S_tc)，t⁺=t/T，ρ为空气密度，S_t为双翅面积。
本文所选用的食蚜蝇的形态学参数(m, I_x等)在文献[17, 19]中已经给出，因此只需计算得到有抬升角时的纵向和横向的气动导数即可确定稳定性矩阵A₁⁺和A₂⁺。基于稳定性矩阵可求出特征值^[17]，相应的特征值代表了相应的简单运动，该运动即为系统的特征模态；处于平衡飞行状态的食蚜蝇，受到小扰动后，其自由运动即为这些特征模态的组合。特征值为正数时(λ>0)表示不稳定发散模态，初始扰动会随时间增长，扰动振幅的倍幅期为0.693/|λ|。特征值为负数时(λ < 0)表示衰减模态，初始扰动会随时间减小最终消失，扰动振幅的半衰期为0.693/|λ|。特征值为复数时，如，扰动会随时间振荡，其振荡周期取决于特征值虚部，扰动的增长或衰减取决于特征值的实部。当时，扰动会增长，表示不稳定振荡模态；反之表示扰动会衰减，表示稳定的振荡模态，振荡模态的倍幅期(或半衰期)和周期分别由和给出。
1.2 流动方程求解与稳定性导数
1.2.1 翅膀运动 本文中食蚜蝇的模型翅膀与文献[3, 12, 17, 19]中用的模型翅膀相同，其计算网格沿展向和弦向的剖面如图 2所示。

图 2 食蚜蝇翅膀网格的展向和弦向剖面 Fig. 2 True hoverfly wing planform and portion of computational grid

图选项

翅膀的各运动学参数定义如图 3所示。图中：β为拍动平面的倾角(见图 3(a))，而翅膀运动的欧拉角?(拍动角)、ψ(翻转角)和θ(抬升角)定义在拍动平面坐标系(O_wX_wY_wZ_w)内如图 3(b)所示，其中O_wx_wy_wz_w为有抬升角θ时固定在食蚜蝇翅膀上的动坐标系，O_wx′_wy′_wz′_w为无抬升角θ时固定在其翅膀上的动坐标系。

图 3 拍动翅的运动学参数定义 Fig. 3 Definition of kinematic parameters of flapping wing

图选项

翅膀转动的3个欧拉角随时间的变化可以用不同的谐波函数表示，拍动角?可用下面的余弦函数表示：

(5)

式中：?为平均拍动角，?=(?_max+?_min)/2，Φ=?_max-?_min，?_max和?_min分别为拍动角的最大值和最小值。
对于翻转角，其在翅膀上/下拍动的中部不随时间改变近似为常数，而在翻转过程中随时间变化，如在第N个拍动周期翅膀由下拍变为上拍时ψ可表示为

(6)

式中：a为常数，Δt_s为翅膀上仰时间，a=(180°-α_u-α_d)/Δt_s，α_d为翅膀下拍中部攻角，α_u为上拍中部攻角；t₁为翅膀开始上仰的时刻，t₁=NT-0.5T-Δt_s/2。翅膀下俯运动时与以上定义类似。
实验观测结果^[3]显示，抬升角周期变化的频率约为拍翅频率的2倍，一般用如下谐波函数^[2-3]来拟合曲线：

(7)

式中: θ₀为偏离拍动平面的平均值; θ₁、θ₂和η₁、η₂分别为第一、二阶项的幅度和相位角。
根据实验测量结果^[3]，Φ=65.6°，n=162 Hz，Δt_s=0.35，θ₀=6°。加上抬升角后，本文通过力和力矩的平衡条件在无抬升角的基础^[19]上对上下拍攻角α_d、α_u和平均拍动角?进行了微调，新的平衡状态下α_d=18.9°，α_u=49°，?=16.9°与无抬升角^[19]时差别不大。

1.2.2 流动方程求解 稳定性导数需要通过求解Navier-Stokes方程来获得，Navier-Stokes方程中的无量纲参数为雷诺数Re(Re=Uc/ν，ν为运动黏性系数)，食蚜蝇悬停飞行时Re≈323^{[3, 19]}。本文求解Navier-Stokes方程的算法基于Rogers等^[22]发展的拟压缩性算法，计算程序与之前的研究所采用的相同^{[16-17, 19, 23]}，能够保证计算结果的可靠性和精度。
本文所使用的翅膀网格沿展向和弦向的剖面形状如图 2所示，网格近壁面法向第1层间距为平均弦长c的1/1 000，翅膀壁面沿法向到远场边界的距离约为20c，沿展向方向到远场边界的距离为6c。文献[3]中已经对网格密度、计算域、时间步长等进行了验证，进一步证明了本文所采用的计算参数和计算方法的可靠性。

1.2.3 稳定性导数计算 本文中，纵向和横向稳定性导数的计算方法与文献[15, 17, 19, 24]中对昆虫纵向和横向稳定性问题研究时所采用的方法相同。以横向稳定性导数Y_v、L_v和N_v为例，简要介绍导数计算。食蚜蝇悬停飞行，其各向力与力矩处于平衡状态时，保持p⁺、r⁺、γ均为0，对v⁺在0附近进行微调，计算得到一系列气动力和力矩相对悬停时的变化量(ΔY⁺、ΔL⁺和ΔN⁺)随v⁺的变化曲线，拟合该系列曲线，各曲线的斜率即为相应的稳定性导数。同理，可以计算得到u⁺、w⁺和q⁺系列等纵向稳定性导数与p⁺和r⁺系列等横向稳定性导数。计算得到稳定性导数后即可确定稳定性矩阵A₁⁺和A₂⁺，进而可以对有抬升角时食蚜蝇的飞行稳定性进行分析。
2 计算结果分析 2.1 稳定性导数比较 如1.2.1节所述，通过计算确定了有抬升角时食蚜蝇的平衡飞行状态，并在此平衡状态下对u⁺、v⁺、w⁺、p⁺、q⁺和r⁺等状态变量进行微调，分别得到了相应的气动力(X⁺、Y⁺、Z⁺)与气动力矩(L⁺、M⁺、N⁺)。这里将有抬升角和无抬升角时的各气动力与气动力矩随各状态变量变化曲线进行对比，结果如图 4和图 5所示。可以发现，有抬升角时各气动力和气动力矩系数随状态变量的变化线性较好，变化趋势也与无抬升角时大致相同。但值得注意的是，如图 5所示，有抬升角时ΔL⁺随横向速度v⁺的变化率比无抬升角时要大。为进行更深入的比较，根据图 4和图 5的变化曲线，计算出了纵向和横向稳定性导数(见表 1和表 2)。通过比较可以发现，有无抬升角2种情形下，所有纵向稳定性导数的差异都很小，且绝大多数的横向动稳定性导数也无明显差异，但是有抬升角时的L_v⁺导数比无抬升角时更小，数值上约为无抬升角时的2倍。有抬升角时L_v⁺为何会有如此明显的变化，其对横向动稳定性是否会产生影响值得进行深入的探讨。

图 4 纵向u+、w+和q+系列气动力和力矩系数 Fig. 4 Longitudinal u+, w+ and q+ series aerodynamic force and moment coefficients

图选项

图 5 横向v+、p+和r+系列气动力和力矩系数 Fig. 5 Lateral v+, p+ and r+ series aerodynamic force and moment coefficients

图选项

首先分析抬升角的存在为何使得L_v⁺导数进一步减小，在此之前需要明确有侧向来流时滚转力矩是如何产生的。侧向来流的存在引起了食蚜蝇左右翅气动力的不对称，这里将计算得到的两翅的气动力系数分别沿Y和Z轴投影，可得其侧向力系数Y⁺和举力系数Z⁺，通过分析比较左右翅Y⁺和Z⁺的差异可以得出滚转力矩系数ΔL⁺产生的原因。以横向扰动运动速度v⁺=0.15为例，将有抬升角和无抬升角2种情形下左右两翅的Y⁺、Z⁺和绕X轴的滚转力矩系数ΔL⁺进行了对比，如表 3所示。通过比较可以发现，2种情形下，右翅产生的沿Y轴负向的侧向力系数数值均大于左翅沿Y轴正向的侧向力系数，左右两翅侧向力系数的差值ΔY⁺分别为-0.210和-0.228，而食蚜蝇的形态学数据测量结果^[3]表明翅根位于身体质心的上方，因此ΔY⁺都会使食蚜蝇产生绕X轴负向的滚转力矩，且有抬升角时的左右两翅侧向力系数的差值ΔY⁺比无抬升角时略大，即有抬升角时绕X轴负向的滚转力矩比无抬升角时大。再者2种情形下，左翅沿Z轴的举力系数均大于右翅，故Z⁺会使食蚜蝇产生绕X轴正向的滚转力矩，但有抬升角时左右两翅举力系数的差值ΔZ⁺仅为无抬升角时的约1/3，这说明抬升角的存在使得左右两翅举力系数的差值减小，从而削弱了食蚜蝇正向的滚转力矩。综上，抬升角的存在使得左右两翅的侧向力差值ΔY⁺增大，增大了食蚜蝇的负向滚转力矩；却减小了左右两翅举力系数的差值ΔZ⁺，即减小了相应的正向滚转力矩，从而导致有抬升角时绕X轴负向的滚转力矩ΔL⁺是无抬升角时的2倍左右，最终使得有抬升角时的横向稳定性导数L_v⁺的数值约为无抬升角时的2倍。
表 1 纵向无量纲稳定性导数 Table 1 Longitudinal non-dimensional stability derivatives

有/无抬升角	Xu+	Zu+	Mu+	Xw+	Zw+	Mw+	Xq+	Zq+	Mq+
无抬升角	-1.998	-0.087	2.378	0.146	-1.777	0.475	-0.206	-0.220	-0.128
有抬升角	-1.519	-0.136	2.449	0.112	-1.829	0.630	-0.155	-0.188	-0.174

表选项






表 2 横向无量纲稳定性导数 Table 2 Lateral non-dimensional stability derivatives

有/无抬升角	Yv+	Lv+	Nv+	Yp+	Lp+	Np+	Yr+	Lr+	Nr+
无抬升角	-0.700	-0.530	0.473	-0.341	-3.727	0.249	0.066	0.805	-3.156
有抬升角	-0.759	-1.100	0.303	-0.495	-3.891	0.204	0.066	0.768	-3.155

表选项






表 3 v⁺=0.15时有无抬升角2种情形下左右两翅的侧向力系数、举力系数和滚转力矩系数 Table 3 Coefficients of lateral force, vertical force and rolling moment of left and right wing at v⁺=0.15 with and without stroke deviation

有/无抬升角	YL+	YR+	ΔY+	ZL+	ZR+	ΔZ+	LL+	LR+	ΔL+
无抬升角	0.093	-0.303	-0.210	1.773	1.690	0.083	2.843	-2.922	-0.079
有抬升角	0.326	-0.554	-0.228	1.780	1.749	0.031	2.961	-3.125	-0.164

表选项






2.2 稳定性分析 通过稳定性导数的比较可以发现，有无抬升角2种情形下除了L_v⁺导数的显著差异外其余稳定性导数差别不大，因此可以预见抬升角存在与否对食蚜蝇的纵向稳定性可能不会产生影响，但是有抬升角时的L_v⁺比无抬升角时小，此前的研究^[19]发现随着L_v⁺的变小，昆虫的横向不稳定性会变弱, L_v⁺的不同可能会对其横向稳定性产生影响。因此，为了进一步分析抬升角对食蚜蝇稳定性的影响，将得到的纵向和横向稳定性导数代入式(3)和式(4)中的纵向稳定性矩阵A₁⁺和横向稳定性矩阵A₂⁺中，并用MATLAB软件求解得到了A₁⁺和A₂⁺的特征值，这里将有抬升角和无抬升角时特征值的对比结果分别列于表 4和表 5中。
表 4 纵向稳定性矩阵A₁⁺的特征值 Table 4 Eigenvalues of longitudinal system matrix A₁⁺

有/无抬升角	模态1λ1	模态2λ2	模态3λ3, 4
无抬升角	-0.169	-0.028	0.065±0.136i
有抬升角	-0.169	-0.028	0.067±0.138i

表选项






表 5 横向稳定性矩阵A₂⁺的特征值 Table 5 Eigenvalues of lateral system matrix A₂⁺

有/无抬升角	模态1λ1	模态2λ2	模态3λ3, 4
无抬升角	-1.269	-0.153	-0.006±0.058i
有抬升角	-1.330	-0.157	-0.003±0.089i

表选项






通过比较表 4中的数据可以发现，有抬升角时食蚜蝇的纵向运动的特征值(特征模态)与无抬升角^[17]时相似，也分别为一个数值较大的负实根λ₁(代表快衰减模态)、一个数值较小的负实根λ₂(代表慢衰减模态)和一对具有正实部的共轭复根λ_{3, 4}(代表不稳定振荡模态)。这与文献[14-15, 25]中对蝇类、蜂类、蛾类等昆虫的纵向稳定性研究结果一致。而由表 5则可以发现，有抬升角时食蚜蝇的横向运动的特征值(特征模态)与无抬升角^[19]时也相似，分别为一个数值较大的负实根λ₁(代表快衰减模态)、一个数值较小的负实根λ₂(代表慢衰减模态)和一对具有极小负实部的共轭复根λ_{3, 4}(代表弱稳定振荡模态)。综上可以发现，抬升角的存在与否除了对食蚜蝇悬停时的气动力有微小影响外，其对食蚜蝇的纵向与横向扰动运动特征模态没有影响，这也说明此前在研究食蚜蝇悬停飞行的气动力与动稳定性问题时只考虑了翅膀的拍动角和攻角而忽略了抬升角这一做法是合理的。
3 结论 1) 有抬升角时，食蚜蝇的各气动力和气动力矩系数随各状态变量的变化规律线性较好，变化率与无抬升角时基本相同。
2) 除L_v⁺导数外，其余稳定性导数与无抬升角时的差异很小，甚至相同；有抬升角时的L_v⁺比无抬升角时更小，数值上约为无抬升角时的2倍；该值变小是由于抬升角的存在使得有侧向来流时左右翅的举力系数差别减小，即减小了左右翅举力不同产生的正滚转力矩；而左右翅侧向力的差值比无抬升角时略有增大，即有抬升角时侧向力产生的负滚转力矩数值略有增大，从而使得总的负滚转力矩变大。
3) 抬升角的存在并未改变食蚜蝇的纵向和横向动稳定性，其特征模态与无抬升角时相同，只是特征值的数值略有差异。因此，研究昆虫悬停飞行的动稳定性问题时，当抬升角在一定范围内只考虑翅膀的拍动角和攻角而忽略了抬升角这一做法是合理的。

参考文献

[1]	ELLINGTON C. The aerodynamics of hovering insect flight.III.Kinematics[J]. Philosophical Transactions of the Royal Society of London.Series B, Biological Sciences, 1984, 305(1122): 41-78.
[2]	LIU Y, SUN M. Wing kinematics measurement and aerodynamics of hovering droneflies[J]. Journal of Experimental Biology, 2008, 211(13): 2014-2025. DOI:10.1242/jeb.016931
[3]	MOU X L, LIU Y P, SUN M. Wing motion measurement and aerodynamics of hovering true hoverflies[J]. Journal of Experimental Biology, 2011, 214(17): 2832-2844. DOI:10.1242/jeb.054874
[4]	ENNOS A R. The kinematics and aerodynamics of the free flight of some diptera[J]. Journal of Experimental Biology, 1989, 142(1): 49-85. DOI:10.1242/jeb.142.1.49
[5]	WILLMOTT A P, ELLINGTON C P. The mechanics of flight in the hawkmoth Manduca sexta.I.Kinematics of hovering and forward flight[J]. Journal of Experimental Biology, 1997, 200(21): 2705-2722. DOI:10.1242/jeb.200.21.2705
[6]	FRY S N, SAYAMAN R, DICKINSON M H. The aerodynamics of free-flight maneuvers in Drosophila[J]. Science, 2003, 300(5618): 495-498. DOI:10.1126/science.1081944
[7]	FRY S N, SAYAMAN R, DICKINSON M H. The aerodynamics of hovering flight in Drosophila[J]. Journal of Experimental Biology, 2005, 208(12): 2303-2318. DOI:10.1242/jeb.01612
[8]	WANG C, ZHOU C Y, XIE P. Numerical investigation into the effects of stroke trajectory on the aerodynamic performance of insect hovering flight[J]. Journal of Mechanical Science and Technology, 2016, 30(4): 1659-1669. DOI:10.1007/s12206-016-0322-3
[9]	LUO G Y, DU G, SUN M. Effects of stroke deviation on aerodynamic force production of a flapping wing[J]. AIAA Journal, 2018, 56(1): 25-35. DOI:10.2514/1.J055739
[10]	HU F J, LIU X M. Effects of stroke deviation on hovering aerodynamic performance of flapping wings[J]. Physics of Fluids, 2019, 31(11): 111901. DOI:10.1063/1.5124916
[11]	LEI M L, LI C Y. The aerodynamic performance of passive wing pitch in hovering flight[J]. Physics of Fluids, 2020, 32(5): 051902. DOI:10.1063/5.0006902
[12]	牟晓蕾, 许娜. 翅尖轨迹对食蚜蝇悬停时气动特性的影响[J]. 北京航空航天大学学报, 2016, 42(12): 2603-2609. MOU X L, XU N. Effect of wing-tip trajectory on aerodynamics of hovering true hoverfly[J]. Journal of Beijing University of Aeronautics and Astronautics, 2016, 42(12): 2603-2609. (in Chinese)
[13]	JANKAUSKI M, DANIEL T L, SHEN I Y. Asymmetries in wing inertial and aerodynamic torques contribute to steering in flying insects[J]. Bioinspiration & Biomimetics, 2017, 12(4): 046001.
[14]	SUN M, XIONG Y. Dynamic flight stability of a hovering bumblebee[J]. Journal of Experimental Biology, 2005, 208(3): 447-59. DOI:10.1242/jeb.01407
[15]	SUN M, WANG J K, XIONG Y. Dynamic flight stability of hovering insects[J]. Acta Mechanica Sinica, 2007, 23(3): 231-246. DOI:10.1007/s10409-007-0068-3
[16]	ZHANG Y, SUN M. Dynamic flight stability of a hovering model insect: Lateral motion[J]. Acta Mechanica Sinica, 2010, 26(2): 175-190. DOI:10.1007/s10409-009-0303-1
[17]	MOU X L, SUN M. Dynamic flight stability of a model hoverfly in inclined-stroke-plane hovering[J]. Journal of Bionic Engineering, 2012, 9(3): 294-303. DOI:10.1016/S1672-6529(11)60123-6
[18]	XU N, SUN M. Lateral dynamic flight stability of a model bumblebee in hovering and forward flight[J]. Journal of Theoretical Biology, 2013, 319: 102-115. DOI:10.1016/j.jtbi.2012.11.033
[19]	XU N, SUN M. Lateral dynamic flight stability of a model hoverfly in normal and inclined stroke-plane hovering[J]. Bioinspiration & Biomimetics, 2014, 9(3): 036019.
[20]	ETKIN B, REID L D. Dynamics of flight: Stability and control[M]. New York: John Wiley and Sons, Inc, 1996.
[21]	TAYLOR G K, THOMAS A L R. Animal flight dynamics.II.Longitudinal stability in flapping flight[J]. Journal of Theoretical Biology, 2002, 214(3): 351-370. DOI:10.1006/jtbi.2001.2470
[22]	ROGERS S E, KWAK D, KIRIS C. Numerical solution of the incompressible Navier-Stokes equations for steady-state and time-dependent problems[C]//27th Aerospace Sciences Meeting, 1989: 0463.
[23]	SUN M, TANG J. Lift and power requirements of hovering flight in Drosophila virilis[J]. Journal of Experimental Biology, 2002, 205(16): 2413-2427. DOI:10.1242/jeb.205.16.2413
[24]	ZHU H J, MENG X G, SUN M. Forward flight stability in a drone-fly[J]. Scientific Reports, 2020, 10(1): 12. DOI:10.1038/s41598-019-55410-5
[25]	CHENG B, DENG X. Translational and rotational damping of flapping flight and its dynamics and stability at hovering[J]. IEEE Transactions on Robotics, 2011, 27(5): 849-864.