CN104199024A - 甚短基线光学测量目标定位方法 - Google Patents

Abstract

本发明公开了一种甚短基线光学测量目标定位方法，属于航天测量与控制领域。该方法依次进行目标定位计算，速度分量计算和精度计算。依据全微分公式，根据测站站址差计算出测量元素，建立基于全微分方程的雷达-光电经纬仪联合定位模型，将经纬仪的高精度测角信息和雷达的高精度测距信息进行数据融合，确定飞行目标的空间位置。

Description

甚短基线光学测量目标定位方法

所属领域 

本发明专利属于航天测量与控制领域，涉及一种基于全微分算法的甚短基线光学测量目标定位方法。 

背景技术

目前，我国航天靶场试验中的光学跟踪测量设备，已逐渐利用多系统跟踪模式，这种多传感器跟踪可增强系统的可靠性和鲁棒性，为获得目标的高精度可靠性信息提供了有力的技术支持。伴随这种新测量体制的出现，相应地为我们的外弹道数据处理方法提出了新的挑战。为了适应后续测控网的发展，我们必须对这种测量体制的测元数据融合技术开展研究，根据多传感器的特性，结合外部环境条件，确定相应的算法和数学模型，完成目标状态的预测，寻找最佳的匹配、融合算法，实现高精度的外弹道数据处理目的。本发明依据全微分公式，根据测站站址差计算出测量元素，建立基于全微分方程的雷达-光电经纬仪联合定位模型，将经纬仪的高精度测角信息和雷达的高精度测距信息进行数据融合，确定飞行目标的空间位置。 

发明内容

在实验中，采用某光学测量设备及附近雷达站设备进行目标跟踪。利用这两台跟踪设备测量数据来完成运载火箭定位及精度结果的确定。 

本发明提出的甚短基线光学测量目标定位方法，具体过程包括如下步骤： 

步骤一：目标定位计算 

设经纬仪测得的方位角和俯仰角分别为A_1i、E_1i，雷达测距为R_2i，目标在发射坐标系中的标准位置坐标为x_i,y_i,z_i，x₀₁,y₀₁,z₀₁为第经纬仪站址坐标，x₀₂,y₀₂,z₀₂为雷达站址坐标，试求在目标在发射坐标系中的联合定位位置参数。 

由 

${\mathrm{tgA}}_{1i}=\frac{(z_i-z_{01})}{(x_i-x_{01})},R_{1i}=\sqrt{{(x_i-x_{01})}^2+{(y_{i1}-y_{01})}^2+{(z_i-z_{01})}^2}$

进行微分得 

$\mathrm{dA}=\frac{1}{R_{1i}\cos E_{1i}}(\sin A_{1i}d{x}_0-\cos A_{1i}{\mathrm{dz}}_0)---\left(1\right)$

由 

$\sin E_{1i}=(y_i-y_{01})/R_{\mathrm{li}},R_{\mathrm{li}}=\sqrt{{(x_i-x_{01})}^2+{(y_{i1}-y_{01})}^2+{(z_i-z_{01})}^2}$

进行微分得 

dE＝(sinE_1i(cosA_1idx₀+sinA_1idz₀)-cosE_1idy₀)/R_1i       (2) 

式中dx₀＝x₀₂-x₀₁，dy₀＝y₀₂-y₀₁，dz₀＝z₀₂-z₀₁

由(1)式和(2)式可得 

A'_2i＝A_1i+dA 

E'_2i＝E_1i+dE 

根据几何关系，可建立目标在发射坐标系的位置参数模型为： 

$\left[\begin{array}{c}x\\ y\\ z\end{array}\right]=\left[\begin{array}{c}{R}_{2i}\cos E_{2i}^{\′}\cos A_{2i}^{\′}\\ R_{2i}\sin E_{2i}^{\′}\\ R_{2i}\cos E_{2i}^{\′}\sin A_{2i}^{\′}\end{array}\right]+\left[\begin{array}{c}{x}_{02}\\ y_{02}\\ z_{02}\end{array}\right]---\left(3\right)$

步骤二：速度分量计算 

式(3)两侧对时间求导，并记R_2i＝R,E′_2i＝E,A′_2i＝A，可导出速度分量满足如下关系： 

$\left\{\begin{array}{c}\stackrel{.}{x}=\stackrel{.}{R}\cos E\cos A-R\stackrel{.}{E}\sin E\cos A-R\stackrel{.}{A}\sin E\sin A\\ \stackrel{.}{y}=\stackrel{.}{R}\sin E+R\stackrel{.}{E}\cos E\\ \stackrel{.}{z}=\stackrel{.}{R}\cos E\sin A-R\stackrel{.}{E}\sin E\sin A+R\stackrel{.}{A}\cos E\cos A\end{array}\right.---\left(4\right)$

或 

$\left(\begin{array}{c}\stackrel{.}{x}\\ \stackrel{.}{y}\\ \stackrel{.}{z}\end{array}\right)=\left(\begin{array}{ccc}\cos E\cos A& -R\sin E\sin A& -R\sin E\cos A\\ \sin E& 0& R\cos E\\ \cos E\sin A& R\cos E\cos A& -R\sin E\sin A\end{array}\right)\left(\begin{array}{c}\stackrel{.}{R}\\ \stackrel{.}{A}\\ \stackrel{.}{E}\end{array}\right)---\left(5\right)$

步骤三：精度计算 

由R,A,E的误差传播到测量系下弹道位置坐标上的误差为 

$\begin{array}{c}\mathrm{\Δ}{x}_c=\left(\frac{\∂x}{\∂R}\right)\mathrm{\ΔR}+\left(\frac{\∂x}{\∂A}\right)\mathrm{\ΔA}+\left(\frac{\∂x}{\∂E}\right)\mathrm{\ΔE}\\ =\mathrm{\Δ}R\cos A\cos E+\mathrm{\Δ}\mathrm{AR}\sin A\cos E-\mathrm{\Δ}\mathrm{ER}\cos A\sin E\end{array}---\left(4\right)$

$\begin{array}{c}\mathrm{\Δ}{y}_c=\left(\frac{\∂y}{\∂R}\right)\mathrm{\ΔR}+\left(\frac{\∂y}{\∂A}\right)\mathrm{\ΔA}+\left(\frac{\∂y}{\∂E}\right)\mathrm{\ΔE}\\ =\mathrm{\Δ}R\sin E+\mathrm{\Δ}\mathrm{ER}\cos E\end{array}---\left(5\right)$

$\begin{array}{c}\mathrm{\Δ}{z}_c=\left(\frac{\∂z}{\∂R}\right)\mathrm{\ΔR}+\left(\frac{\∂z}{\∂A}\right)\mathrm{\ΔA}+\left(\frac{\∂z}{\∂E}\right)\mathrm{\ΔE}\\ =\mathrm{\Δ}R\sin A\cos E+\mathrm{\Δ}\mathrm{AR}\cos A\cos E-\mathrm{\Δ}\mathrm{ER}\sin A\sin E\end{array}---\left(6\right)$

其中ΔR,ΔA,ΔE分别为R,A,E上的误差，记σ_R ²,σ_A ²,σ_E ²分别为R,A,E误差的方差.用矩阵表示为 

$\mathrm{\Δ}{X}_c=\left(\begin{array}{c}\mathrm{\Δ}{x}_c\\ \mathrm{\Δ}{y}_c\\ \mathrm{\Δ}{z}_c\end{array}\right)=\left(\begin{array}{ccc}\cos A\cos E& -R\sin A\cos E& R\cos A\sin E\\ \sin E& 0& R\cos E\\ \sin A\cos E& R\cos A\cos E& -R\sin E\sin A\end{array}\right)\left(\begin{array}{c}\mathrm{\ΔR}\\ \mathrm{\ΔA}\\ \mathrm{\ΔE}\end{array}\right)---\left(7\right)$

记测元误差传播到测量系下弹道位置参数的传播矩阵为 

$G=\left(\begin{array}{ccc}\cos A\cos E& -R\sin A\cos E& -R\cos A\sin E\\ \sin E& 0& R\cos E\\ \sin A\cos E& R\cos A\cos E& -R\sin A\sin E\end{array}\right)$

则测量系下弹道位置参数的协方差矩阵为 

$\mathrm{COV}\left(\mathrm{\Δ}{X}_c\right)=G\·\mathrm{diag}({{\mathrm{\σ}}_R}^2,{{\mathrm{\σ}}_A}^2,{{\mathrm{\σ}}_E}^2)\·G^T\widehat{=}\left(\begin{array}{ccc}{{\mathrm{\σ}}_{x_c}}^2& {\mathrm{\σ}}_{x_c{y}_c}& {\mathrm{\σ}}_{x_c{z}_c}\\ {\mathrm{\σ}}_{y_c{x}_c}& {{\mathrm{\σ}}_{y_c}}^2& {\mathrm{\σ}}_{y_c{z}_c}\\ {\mathrm{\σ}}_{z_c{x}_c}& {\mathrm{\σ}}_{z_c{y}_c}& {{\mathrm{\σ}}_{z_c}}^2\end{array}\right)---\left(8\right)$

其中 

$\left\{\begin{array}{c}{\mathrm{\σ}}_{x_c}^2={{\mathrm{\σ}}_R}^2{\cos }^{2}A{\cos }^{2}E+{{\mathrm{\σ}}_A}^2{R}^2{\sin }^{2}A{\cos }^{2}E+{\mathrm{\σ}}_E^2{R}^2{\cos }^{2}A{\sin }^{2}E\\ {\mathrm{\σ}}_{y_c}^2={{\mathrm{\σ}}_R}^2{\sin }^{2}E+{{\mathrm{\σ}}_E}^2{R}^2{\cos }^{2}E\\ {\mathrm{\σ}}_{z_c}^2={{\mathrm{\σ}}_R}^2{\sin }^{2}A{\cos }^{2}E+{{\mathrm{\σ}}_A}^2{R}^2{\cos }^{2}A{\cos }^{2}E+{{\mathrm{\σ}}_E}^2{R}^2{\sin }^{2}A{\sin }^{2}E\end{array}\right.$

$\left\{\begin{array}{c}{\mathrm{\σ}}_{x_c{y}_c}={\mathrm{\σ}}_R^2\cos A\cos E\sin E-{\mathrm{\σ}}_E^2{R}^2\cos A\sin E\cos E\\ {\mathrm{\σ}}_{x_c{z}_c}={\mathrm{\σ}}_R^2\sin A\cos A{\cos }^{2}E-{\mathrm{\σ}}_A^2{R}^2\sin A\cos A{\cos }^{2}E\\ +{\mathrm{\σ}}_E^2{R}^2\sin A\cos A{\sin }^{2}E\\ {\mathrm{\σ}}_{y_c{z}_c}={\mathrm{\σ}}_R^2\sin A\cos E\sin E-{\mathrm{\σ}}_E^2{R}^2\sin A\sin E\cos E\\ {\mathrm{\σ}}_{y_c{x}_c}={\mathrm{\σ}}_{x_c{y}_c}\\ {\mathrm{\σ}}_{z_c{x}_c}={\mathrm{\σ}}_{x_c{z}_c}\\ {\mathrm{\σ}}_{z_c{y}_c}={\mathrm{\σ}}_{y_c{z}_c}\end{array}\right.$

测元误差传播到发射系下位置参数的误差为 

$\mathrm{\ΔX}=\left(\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\\ \mathrm{\Δz}\end{array}\right)=\mathrm{M\Δ}{X}_c=\mathrm{MG}\left(\begin{array}{c}\mathrm{\ΔR}\\ \mathrm{\ΔA}\\ \mathrm{\ΔE}\end{array}\right)---\left(9\right)$

其中M＝(m_ij)_3*3为测量系到发射系的变换矩阵.则发射系下弹道位置参数的协方差为 

$\left\{\begin{array}{c}{{\mathrm{\σ}}_x}^2=m_{11}^2{\mathrm{\σ}}_{x_c}^2+m_{12}^2{\mathrm{\σ}}_{y_c}^2+m_{13}^2{\mathrm{\σ}}_{z_c}^2+2m_{11}{m}_{12}{\mathrm{\σ}}_{x_c{y}_c}+2m_{11}{m}_{13}{\mathrm{\σ}}_{x_c{z}_c}+2m_{12}{m}_{13}{\mathrm{\σ}}_{y_c{z}_c}\\ {{\mathrm{\σ}}_y}^2=m_{21}^2{\mathrm{\σ}}_{x_c}^2+m_{22}^2{\mathrm{\σ}}_{y_c}^2+m_{23}^2{\mathrm{\σ}}_{z_c}^2+2m_{21}{m}_{22}{\mathrm{\σ}}_{x_c{y}_c}+2m_{21}{m}_{23}{\mathrm{\σ}}_{x_c{z}_c}+2m_{22}{m}_{23}{\mathrm{\σ}}_{y_c{z}_c}\\ {{\mathrm{\σ}}_z}^2=m_{31}^2{\mathrm{\σ}}_{x_c}^2+m_{32}^2{\mathrm{\σ}}_{y_c}^2+m_{33}^2{\mathrm{\σ}}_{z_c}^2+2m_{31}{m}_{32}{\mathrm{\σ}}_{x_c{y}_c}+2m_{31}{m}_{33}{\mathrm{\σ}}_{x_c{z}_c}+2m_{32}{m}_{33}{\mathrm{\σ}}_{y_c{z}_c}\end{array}\right.---\left(10\right)$

具体实施方式

下面以具体实例说明，分别以某次试验任务某甚短基线光学测量设备及某雷达设备测量数据目标定位计算、速度分量计算和精度计算为例，给出具体实施过程。 

1.目标定位计算 

根据计算方法，得到计算结果如下(受篇幅限制，只给出部分计算结果,其中T为相对时，单位为秒，下同)： 

2.速度分量计算 

3.精度计算 

。

Claims (1)

1.甚短基线光学测量目标定位方法，其特征在于：具体过程包括如下步骤：

步骤一：目标定位计算

设经纬仪测得的方位角和俯仰角分别为A_1i、E_1i，雷达测距为R_2i，目标在发射坐标系中的标准位置坐标为x_i,y_i,z_i，x₀₁,y₀₁,z₀₁为第经纬仪站址坐标，x₀₂,y₀₂,z₀₂为雷达站址坐标，

由

${\mathrm{tgA}}_{1i}=\frac{(z_i-z_{01})}{(x_i-x_{01})},R_{1i}=\sqrt{{(x_i-x_{01})}^2+{(y_{i1}-y_{01})}^2+{(z_i-z_{01})}^2}$

进行微分得

$\mathrm{dA}=\frac{1}{R_{1i}\cos E_{1i}}(\sin A_{1i}d{x}_0-\cos A_{1i}{\mathrm{dz}}_0)---\left(1\right)$

由

$\sin E_{1i}=(y_i-y_{01})/R_{\mathrm{li}},R_{\mathrm{li}}=\sqrt{{(x_i-x_{01})}^2+{(y_{i1}-y_{01})}^2+{(z_i-z_{01})}^2}$

进行微分得

dE＝(sinE_1i(cosA_1idx₀+sinA_1idz₀)-cosE_1idy₀)/R_1i             (2)

式中dx₀＝x₀₂-x₀₁，dy₀＝y₀₂-y₀₁，dz₀＝z₀₂-z₀₁

由(1)式和(2)式可得

A'_2i＝A_1i+dA

E'_2i＝E_1i+dE

根据几何关系，可建立目标在发射坐标系的位置参数模型为：

$\left[\begin{array}{c}x\\ y\\ z\end{array}\right]=\left[\begin{array}{c}{R}_{2i}\cos E_{2i}^{\′}\cos A_{2i}^{\′}\\ R_{2i}\sin E_{2i}^{\′}\\ R_{2i}\cos E_{2i}^{\′}\sin A_{2i}^{\′}\end{array}\right]+\left[\begin{array}{c}{x}_{02}\\ y_{02}\\ z_{02}\end{array}\right]---\left(3\right)$

步骤二：速度分量计算

式(3)两侧对时间求导，并记R_2i＝R,E′_2i＝E,A′_2i＝A，可导出速度分量满足如下关系：

$\left\{\begin{array}{c}\stackrel{.}{x}=\stackrel{.}{R}\cos E\cos A-R\stackrel{.}{E}\sin E\cos A-R\stackrel{.}{A}\sin E\sin A\\ \stackrel{.}{y}=\stackrel{.}{R}\sin E+R\stackrel{.}{E}\cos E\\ \stackrel{.}{z}=\stackrel{.}{R}\cos E\sin A-R\stackrel{.}{E}\sin E\sin A+R\stackrel{.}{A}\cos E\cos A\end{array}\right.---\left(4\right)$

或

$\left(\begin{array}{c}\stackrel{.}{x}\\ \stackrel{.}{y}\\ \stackrel{.}{z}\end{array}\right)=\left(\begin{array}{ccc}\cos E\cos A& -R\sin E\sin A& -R\sin E\cos A\\ \sin E& 0& R\cos E\\ \cos E\sin A& R\cos E\cos A& -R\sin E\sin A\end{array}\right)\left(\begin{array}{c}\stackrel{.}{R}\\ \stackrel{.}{A}\\ \stackrel{.}{E}\end{array}\right)---\left(5\right)$

步骤三：精度计算

由R,A,E的误差传播到测量系下弹道位置坐标上的误差为

$\begin{array}{c}\mathrm{\Δ}{x}_c=\left(\frac{\∂x}{\∂R}\right)\mathrm{\ΔR}+\left(\frac{\∂x}{\∂A}\right)\mathrm{\ΔA}+\left(\frac{\∂x}{\∂E}\right)\mathrm{\ΔE}\\ =\mathrm{\Δ}R\cos A\cos E+\mathrm{\Δ}\mathrm{AR}\sin A\cos E-\mathrm{\Δ}\mathrm{ER}\cos A\sin E\end{array}---\left(4\right)$

$\begin{array}{c}\mathrm{\Δ}{y}_c=\left(\frac{\∂y}{\∂R}\right)\mathrm{\ΔR}+\left(\frac{\∂y}{\∂A}\right)\mathrm{\ΔA}+\left(\frac{\∂y}{\∂E}\right)\mathrm{\ΔE}\\ =\mathrm{\Δ}R\sin E+\mathrm{\Δ}\mathrm{ER}\cos E\end{array}---\left(5\right)$

$\begin{array}{c}\mathrm{\Δ}{z}_c=\left(\frac{\∂z}{\∂R}\right)\mathrm{\ΔR}+\left(\frac{\∂z}{\∂A}\right)\mathrm{\ΔA}+\left(\frac{\∂z}{\∂E}\right)\mathrm{\ΔE}\\ =\mathrm{\Δ}R\sin A\cos E+\mathrm{\Δ}\mathrm{AR}\cos A\cos E-\mathrm{\Δ}\mathrm{ER}\sin A\sin E\end{array}---\left(6\right)$

其中ΔR,ΔA,ΔE分别为R,A,E上的误差，记σ_R ²,σ_A ²,σ_E ²分别为R,A,E误差的方差.用矩阵表示为

$\mathrm{\Δ}{X}_c=\left(\begin{array}{c}\mathrm{\Δ}{x}_c\\ \mathrm{\Δ}{y}_c\\ \mathrm{\Δ}{z}_c\end{array}\right)=\left(\begin{array}{ccc}\cos A\cos E& -R\sin A\cos E& R\cos A\sin E\\ \sin E& 0& R\cos E\\ \sin A\cos E& R\cos A\cos E& -R\sin E\sin A\end{array}\right)\left(\begin{array}{c}\mathrm{\ΔR}\\ \mathrm{\ΔA}\\ \mathrm{\ΔE}\end{array}\right)---\left(7\right)$

记测元误差传播到测量系下弹道位置参数的传播矩阵为

$G=\left(\begin{array}{ccc}\cos A\cos E& -R\sin A\cos E& -R\cos A\sin E\\ \sin E& 0& R\cos E\\ \sin A\cos E& R\cos A\cos E& -R\sin A\sin E\end{array}\right)$

则测量系下弹道位置参数的协方差矩阵为

$\mathrm{COV}\left(\mathrm{\Δ}{X}_c\right)=G\·\mathrm{diag}({{\mathrm{\σ}}_R}^2,{{\mathrm{\σ}}_A}^2,{{\mathrm{\σ}}_E}^2)\·G^T\widehat{=}\left(\begin{array}{ccc}{{\mathrm{\σ}}_{x_c}}^2& {\mathrm{\σ}}_{x_c{y}_c}& {\mathrm{\σ}}_{x_c{z}_c}\\ {\mathrm{\σ}}_{y_c{x}_c}& {{\mathrm{\σ}}_{y_c}}^2& {\mathrm{\σ}}_{y_c{z}_c}\\ {\mathrm{\σ}}_{z_c{x}_c}& {\mathrm{\σ}}_{z_c{y}_c}& {{\mathrm{\σ}}_{z_c}}^2\end{array}\right)---\left(8\right)$

其中

$\left\{\begin{array}{c}{\mathrm{\σ}}_{x_c}^2={{\mathrm{\σ}}_R}^2{\cos }^{2}A{\cos }^{2}E+{{\mathrm{\σ}}_A}^2{R}^2{\sin }^{2}A{\cos }^{2}E+{\mathrm{\σ}}_E^2{R}^2{\cos }^{2}A{\sin }^{2}E\\ {\mathrm{\σ}}_{y_c}^2={{\mathrm{\σ}}_R}^2{\sin }^{2}E+{{\mathrm{\σ}}_E}^2{R}^2{\cos }^{2}E\\ {\mathrm{\σ}}_{z_c}^2={{\mathrm{\σ}}_R}^2{\sin }^{2}A{\cos }^{2}E+{{\mathrm{\σ}}_A}^2{R}^2{\cos }^{2}A{\cos }^{2}E+{{\mathrm{\σ}}_E}^2{R}^2{\sin }^{2}A{\sin }^{2}E\end{array}\right.$

$\left\{\begin{array}{c}{\mathrm{\σ}}_{x_c{y}_c}={\mathrm{\σ}}_R^2\cos A\cos E\sin E-{\mathrm{\σ}}_E^2{R}^2\cos A\sin E\cos E\\ {\mathrm{\σ}}_{x_c{z}_c}={\mathrm{\σ}}_R^2\sin A\cos A{\cos }^{2}E-{\mathrm{\σ}}_A^2{R}^2\sin A\cos A{\cos }^{2}E\\ +{\mathrm{\σ}}_E^2{R}^2\sin A\cos A{\sin }^{2}E\\ {\mathrm{\σ}}_{y_c{z}_c}={\mathrm{\σ}}_R^2\sin A\cos E\sin E-{\mathrm{\σ}}_E^2{R}^2\sin A\sin E\cos E\\ {\mathrm{\σ}}_{y_c{x}_c}={\mathrm{\σ}}_{x_c{y}_c}\\ {\mathrm{\σ}}_{z_c{x}_c}={\mathrm{\σ}}_{x_c{z}_c}\\ {\mathrm{\σ}}_{z_c{y}_c}={\mathrm{\σ}}_{y_c{z}_c}\end{array}\right.$

测元误差传播到发射系下位置参数的误差为

$\mathrm{\ΔX}=\left(\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\\ \mathrm{\Δz}\end{array}\right)=\mathrm{M\Δ}{X}_c=\mathrm{MG}\left(\begin{array}{c}\mathrm{\ΔR}\\ \mathrm{\ΔA}\\ \mathrm{\ΔE}\end{array}\right)---\left(9\right)$

其中M＝(m_ij)_3*3为测量系到发射系的变换矩阵.则发射系下弹道位置参数的协方差为

$\left\{\begin{array}{c}{{\mathrm{\σ}}_x}^2=m_{11}^2{\mathrm{\σ}}_{x_c}^2+m_{12}^2{\mathrm{\σ}}_{y_c}^2+m_{13}^2{\mathrm{\σ}}_{z_c}^2+2m_{11}{m}_{12}{\mathrm{\σ}}_{x_c{y}_c}+2m_{11}{m}_{13}{\mathrm{\σ}}_{x_c{z}_c}+2m_{12}{m}_{13}{\mathrm{\σ}}_{y_c{z}_c}\\ {{\mathrm{\σ}}_y}^2=m_{21}^2{\mathrm{\σ}}_{x_c}^2+m_{22}^2{\mathrm{\σ}}_{y_c}^2+m_{23}^2{\mathrm{\σ}}_{z_c}^2+2m_{21}{m}_{22}{\mathrm{\σ}}_{x_c{y}_c}+2m_{21}{m}_{23}{\mathrm{\σ}}_{x_c{z}_c}+2m_{22}{m}_{23}{\mathrm{\σ}}_{y_c{z}_c}\\ {{\mathrm{\σ}}_z}^2=m_{31}^2{\mathrm{\σ}}_{x_c}^2+m_{32}^2{\mathrm{\σ}}_{y_c}^2+m_{33}^2{\mathrm{\σ}}_{z_c}^2+2m_{31}{m}_{32}{\mathrm{\σ}}_{x_c{y}_c}+2m_{31}{m}_{33}{\mathrm{\σ}}_{x_c{z}_c}+2m_{32}{m}_{33}{\mathrm{\σ}}_{y_c{z}_c}\end{array}\right.---\left(10\right)$