问题
已知无人机在 NED 下的 ZYX 欧拉角为 [ ψ , ϕ , θ ] [\psi, \phi, \theta] [ψ,ϕ,θ],分别表示偏航角、滚转角和俯仰角。请问该姿态在 ENU 下的表示 [ ψ ′ , ϕ ′ , θ ′ ] [\psi', \phi', \theta'] [ψ′,ϕ′,θ′] 是多少。
求解
直接使用欧拉角难以描述连续两次刚体旋转,此处使用旋转矩阵推导。
从Body系到NED的旋转矩阵为
R b n = [ cos ψ cos θ cos ψ sin θ sin ϕ − sin ψ cos ϕ cos ψ sin θ cos ϕ + sin ψ sin ϕ sin ψ cos θ sin ψ sin θ sin ϕ + cos ψ cos ϕ sin ψ sin θ cos ϕ − cos ψ sin ϕ − sin θ cos θ sin ϕ cos θ cos ϕ ] R_{b}^{n}=\left[\begin{array}{ccc} \cos \psi \cos \theta & \cos \psi \sin \theta \sin \phi-\sin \psi \cos \phi & \cos \psi \sin \theta \cos \phi+\sin \psi \sin \phi \\ \sin \psi \cos \theta & \sin \psi \sin \theta \sin \phi+\cos \psi \cos \phi & \sin \psi \sin \theta \cos \phi-\cos \psi \sin \phi \\ -\sin \theta & \cos \theta \sin \phi & \cos \theta \cos \phi \end{array}\right] Rbn=⎣ ⎡cosψcosθsinψcosθ−sinθcosψsinθsinϕ−sinψcosϕsinψsinθsinϕ+cosψcosϕcosθsinϕcosψsinθcosϕ+sinψsinϕsinψsinθcosϕ−cosψsinϕcosθcosϕ⎦ ⎤
从 NED 到 ENU 的旋转矩阵为
R n e = [ 0 1 0 1 0 0 0 0 − 1 ] R_{n}^{e}=\left[\begin{array}{ccc} 0 &1 & 0\\ 1 & 0 & 0\\ 0 & 0 & -1 \end{array}\right] Rne=⎣ ⎡01010000−1⎦ ⎤
因此,从 Body 到 ENU 的旋转矩阵为
R b e = R n e R b n R_b^e = R_n^e R _b^n Rbe=RneRbn
最后再根据旋转矩阵计算姿态角。
ψ = a r c t a n 2 ( r 21 , r 11 ) ϕ = arctan ( r 32 r 33 ) θ = arcsin ( − r 31 ) \begin{aligned} \psi &= {\rm arctan2}\left({r_{21}},{r_{11}} \right)\\ \phi &= \arctan \left(\frac{r_{32}}{r_{33}} \right) \\ \theta &= \arcsin(-r_{31}) \end{aligned} ψϕθ=arctan2(r21,r11)=arctan(r33r32)=arcsin(−r31)
注意, a r c t a n 2 ( ) {\rm arctan2}() arctan2() 范围为 [ − π , π ] [-\pi, \pi] [−π,π], a r c s i n ( ) {\rm arcsin}() arcsin() 范围为 [ − π / 2 , π / 2 ] [-\pi/2, \pi/2] [−π/2,π/2], a r c t a n ( ) {\rm arctan}() arctan() 范围为 [ − π / 2 , π / 2 ] [-\pi/2, \pi/2] [−π/2,π/2]。
根据以上原理,使用MATLAB符号表达式求解,结果如下
ψ ′ = a r c t a n 2 ( 1 , tan ψ ) ϕ ′ = ϕ θ ′ = − θ \begin{aligned} \psi' &= {\rm arctan2}({1},{\tan\psi} ) \\ \phi' &= \phi \\ \theta' &= -\theta \end{aligned} ψ′ϕ′θ′=arctan2(1,tanψ)=ϕ=−θ
代码
%% 计算
syms psi phi theta real
% Body --> NED
Rbn = [cos(psi)*cos(theta) cos(psi)*sin(theta)*sin(phi)-sin(psi)*cos(phi) cos(psi)*sin(theta)*cos(phi)+sin(psi)*sin(phi)
sin(psi)*cos(theta) sin(psi)*sin(theta)*sin(phi)+cos(psi)*cos(phi) sin(psi)*sin(theta)*cos(phi)-cos(psi)*sin(phi)
-sin(theta) cos(theta)*sin(phi) cos(theta)*cos(phi)];
% NED --> ENU
Rne = [0 1 0
1 0 0
0 0 -1];
% Body --> ENU
Rbe = Rne * Rbn;
% Rotate Matrix --> Eular Angle
psi1 = simplify(atan2(Rbe(2,1), Rbe(1,1))) % psi in [-180°,180°]
phi1 = simplify(atan(Rbe(3,2) / Rbe(3,3)))
theta1 = simplify(asin(-Rbe(3,1)))
% 结果
% psi1 = angle(cos(theta)*(sin(psi) + cos(psi)*1i)) = atan2(1,tan(psi))
% phi1 = atan(tan(phi)) = phi
% theta1 = -asin(sin(theta)) = -theta
%% 验证
psi = -170 * pi/180 ;
phi = 20 * pi / 180;
theta = 10 * pi / 180;
disp(['NED psi = ' num2str(psi*180/pi), ', ENU psi = ' num2str(eval(psi1)*180/pi)])
disp(['NED phi = ' num2str(phi*180/pi), ', ENU phi = ' num2str(eval(phi1)*180/pi)])
disp(['NED theta = ' num2str(theta*180/pi), ', ENU theta = ' num2str(eval(theta1)*180/pi)])