一种多约束的风洞试验模型变形视频测量振动修正方法

摘要

本发明公开了一种多约束的风洞试验模型变形视频测量振动修正方法，包括如下步骤：布置视频观测相机；建立像点投影观测方程；建立刚体和极线约束方程；构建迭代方程；迭代求取方程组最优解；重构当前相机姿态和模型姿态。本发明的积极效果是：既利用了试验模型机身特征点的刚体约束条件、又利用了风洞试验段壁板特征点的刚体约束条件，进一步地增加了试验模型机翼特征点在多视测量几何中的共线约束条件，创造性地将所有的约束条件变为统一形式的迭代方程，具有求解迭代收敛快、求解精度高，可以同时求解得到试验模型的位置和姿态等优点。

说明书

技术领域

本发明涉及计算机视觉及摄影测量学领域，尤其涉及一种多约束的风洞试验模型变形视频测量振动修正方法，是机器视觉测量技术在高速风洞试验中的工程应用。 

背景技术

风洞试验过程中，试验模型机翼受气动载荷作用而发生弹性变形。对大展弦比机翼来说，这种变形将严重影响飞行器的气动特性，偏离设计需要，必须要对试验数据加以修正，因此需要在风洞试验的过程中准确地测量出试验模型的变形量大小。而利用非接触式的光学测量方法，对风洞流场干扰小、测量精度高，不改变试验模型的几何外形，是目前国际上应用最为广泛的风洞试验模型变形测量方法。 

在风洞中应用光学测量技术，特别是在高速风洞中应用视频测量技术，一个必须解决的关键技术问题是如何处理风洞现场振动导致相机姿态发生变化致使测量精度严重下降的问题。为了保证在整个风洞试验过程中测量数据的准确可靠，必须对试验过程中的每一帧图像都进行振动修正。文献“Demonstration of Optical Wing Deformation Measurements at the Arnold Engineering Development Center”（WimRuyten,Marvin Sellers.47th AIAA Aerospace Sciences Meeting including The New Horizons Forum and Aerospace Exposition.2009.）中阐述了一种利用立体视觉测量技术测量风洞 试验模型变形量的方法，为了修正风洞试验段振动造成的相机姿态变化，该方法在吹风过程中将试验模型的机身视为理想刚体，因此，当试验模型的姿态发生变化时，模型机身上的特征标记点空间相对位置应该保持固定不变，利用这种不变的约束关系以及模型在静止状态下标记点相对位置的初值，即可计算出振动导致的相机之间的位姿变化偏差。然而，仅仅利用模型机身的刚体约束条件，只能修正两台相机间的相对位置姿态，不能修正相机坐标系在试验段世界坐标系中的整体偏移，因而无法得到试验模型的当前姿态和位置。同时，由于试验模型的机身只占有效测量区域的一小部分，该方法并没有充分利用全部测量特征，且该方法对相机的内部参数十分敏感，任何微小的相机内部参数偏差都将导致修正结果产生较大的差异。文献“模型变形视频测量的相机位置坐标与姿态角确定”（罗川,张征宇,孙岩等.实验流体力学.2010,(6):88-91）阐述了使用角锥法来修正单个相机相对于风洞试验段世界坐标系的偏移量，该方法利用试验段壁板上特征标记点的图像坐标单独修正每个相机的振动偏差，但该方法并没有有效地利用各个相机和试验模型之间的隐含约束关系，因此在修正相机姿态的振动偏差后还需要额外的步骤来重构试验模型的位置和姿态。 

发明内容

为了克服现有技术的上述缺点，本发明提供了一种多约束的风洞试验模型变形视频测量振动修正方法，既利用了试验模型机身特征点的刚体约束条件、又利用了风洞试验段壁板特征点的刚体约束条件，更进一步地增加了试验模型机翼特征点在多视测量几何中的共线约 束条件，创造性地将所有的约束条件变为统一形式的迭代方程，具有求解迭代收敛快、求解精度高，可以同时求解得到试验模型的位置和姿态等优点。 

本发明解决其技术问题所采用的技术方案是：一种多约束的风洞试验模型变形视频测量振动修正方法，包括如下步骤： 

步骤一、布置视频观测相机： 

使用两台观测相机，对试验模型的测量部件进行交汇成像，试验时两台相机同步触发采集，两台相机能够同时观测到预先布置在风洞试验段壁板上和试验模型上的所有特征标记点； 

步骤二、建立像点投影观测方程； 

步骤三、建立刚体和极线约束方程： 

步骤四、构建迭代方程； 

步骤五、迭代求取方程组最优解； 

步骤六、重构当前相机姿态和模型姿态。 

与现有技术相比，本发明的积极效果是： 

（1）同时利用了试验模型机身点刚体约束、试验模型机翼点极线约束和试验段壁板点刚体约束三个约束条件，并创造性地将其变换为统一形式的多变量方程，以便于后续的迭代优化求解计算过程。可充分利用风洞试验段和试验模型的各种隐含约束条件，约束方程和求解参数可根据现场情况灵活选择；在实际工程应用过程中，可以根据风洞现场的情况和试验模型的几何外形灵活选择约束条件，简化迭代方程和计算过程。 

（2）可用于修正相机坐标系在风洞试验段坐标系中的振动偏差和两台相机相互之间的振动偏差，并可以同时计算得到试验模型在试验段坐标系中的实时姿态和位置。 

附图说明

本发明将通过例子并参照附图的方式说明，其中： 

图1是观测相机在试验段上的位姿布置示意图； 

图2是风洞吹风时观测相机在试验段中的振动情况示意图，图2中的符号：C1、C2为观测相机坐标系、定义C1与相机坐标系重合、W为试验段坐标系、M为试验模型坐标系。 

具体实施方式

本发明的技术方案是： 

第一、根据立体视觉的测量原理对风洞试验模型进行交汇成像；试验前对系统进行标定；在试验过程中，假定相机的焦距、畸变、图像主点等内部参数保持固定不变；同时假定受现场气流和振动影响，观测相机和试验模型的姿态均相对试验段坐标系发生偏移变化。 

第二、以静态时试验模型零状态下所有特征标记点的三维初始坐标作为参考，将试验模型机身和试验段壁板视为两个独立刚体，试验过程中其特征标记点的三维相对位置坐标保持固定不变；将试验模型机翼视为一个弹性体，根据立体视觉测量原理，试验过程中其特征点的图像坐标将始终满足极线约束条件。 

第三、对试验过程中的每一帧图像，以摄影测量共线方程作为观测方程，以模型机身刚体约束、试验段壁板刚体约束以及模型机翼极 线约束作为约束条件，以相机实时姿态和模型实时姿态作为未知变量，以零状态全部标记点三维坐标和相机初始姿态作为初值，用光束平差的方法求解上述方程的最优解。 

如图1和图2所示，一种多约束的风洞试验模型变形视频测量振动修正方法，具体包括如下步骤： 

步骤一、布置视频观测相机： 

使用两台观测相机（记为1#相机和2#相机），对试验模型的测量部件进行交汇成像，试验时两台相机同步触发采集，要求两台相机能够同时观测到预先布置在风洞试验段壁板上和试验模型上的所有特征标记点。特征标记点的布置方法如下： 

使用圆形特征标记点，将试验模型的机身视为一个刚体，在其上布置3个以上的特征标记点，标记点的位置分布不能全部在同一条直线上；将试验模型的机翼视为一个弹性体，在每个测量截面上成直线布置2个以上的特征标记点；将试验段壁板底部也视为一个刚体，在其上布置3个以上的特征标记点，标记点的位置分布位置也不能全部在同一条直线上。 

步骤二、建立像点投影观测方程： 

在风洞试验模型的零状态，采集并计算所有特征点的初值图像坐标和三维坐标，设特征点对应试验段坐标系下的三维坐标初值为两台相机之间的相对位姿关系矩阵为定义相机坐标系为1#相机坐标系，特征点在相机坐标系下的三维坐标初值为风洞试验过程中由于现场振动影响，相机在试验段中的位置已经偏离了初始值， 设第n帧图像相机之间新的位姿关系为特征点在相机坐标系下的三维坐标为对应试验段坐标系下的三维坐标为

对任一特征点P，设其在1#相机和2#相机坐标系下的三维坐标分别为P₁=[X₁ Y₁ Z₁]^T和P₂=[X₂ Y₂ Z₂]^T，其满足关系对应畸变矫正后的归一化像点坐标为（归一化公式为 其中x、y为特征点图像坐标，u、v为图像主点，f_x、f_y为等效焦距）。 

根据以上定义以及相机小孔成像原理，将特征点P的共线方程简化如下： 

$\left(\begin{array}{c}{\bar{x}}_1=X_1/Z_1\\ {\bar{y}}_1=Y_1/Z_1\\ {\bar{x}}_2=X_2/Z_2\\ {\bar{y}}_2=Y_2/Z_2\end{array}\right)$

改写为矩阵形式： 

$\left(\begin{array}{c}{M}_1{P}_1=0\\ M_2(R_c^n{P}_1+T_c^n)=0\end{array}\right)---\left(0.1\right)$

式中$M_1=\left(\begin{array}{ccc}1& 0& -{\bar{x}}_1\\ 0& 1& -{\bar{y}}_1\end{array}\right),$$M_2=\left(\begin{array}{ccc}1& 0& -{\bar{x}}_2\\ 0& 1& -{\bar{y}}_2\end{array}\right).$

上述方程对所有时刻所有特征点均成立： 

零状态时， 

$\left(\begin{array}{c}{M}_1^i{P}_1^0=0\\ M_2^i(R_c^0{P}_1^0+T_c^0)=0\end{array}\right)$

第n帧图像， 

$\left(\begin{array}{c}{M}_1^i{P}_1^n=0\\ M_2^i(R_c^n{P}_1^n+T_c^n)=0\end{array}\right)---\left(0.2\right)$

其中，$M_1^i=\left(\begin{array}{ccc}1& 0& -{\bar{x}}_1^i\\ 0& 1& -{\bar{y}}_1^i\end{array}\right),$$M_2^i=\left(\begin{array}{ccc}1& 0& -{\bar{x}}_2^i\\ 0& 1& -{\bar{y}}_2^i\end{array}\right).$

步骤三、建立刚体和极线约束方程： 

（1）对试验模型机身上特征标记点的三维坐标，其满足刚体约束条件，对一个刚体来说，一共只有旋转和平移6个自由度。设风洞吹风时模型机身在相机坐标系下观察到的位置和姿态变化为则可以得到如下刚体约束方程： 

$P_{\mathrm{ic}}^n=R_{\mathrm{mc}}^n{P}_{\mathrm{ic}}^0+T_{\mathrm{mc}}^n---\left(0.3\right)$

其中，i表示模型机身特征点，表示模型零状态时所有特征点在相机坐标系下的三维坐标值。 

将试验段坐标系下相机坐标系的位置和姿态变化表示为Rⁿ、Tⁿ，试验模型的位置和姿态变化表示为则有如下关系成立： 

$\left(\begin{array}{c}{R}_m^n=R^n{R}_{\mathrm{mc}}^n\\ T_m^n=R^n{T}_{\mathrm{mc}}^n+T^n\end{array}\right)---\left(0.4\right)$

（2）对试验段壁板特征标记点，其试验段坐标系下的三维坐标保持固定不变，相应相机坐标系下的三维坐标也满足刚体约束条件： 

$R^n{P}_{\mathrm{jc}}^0+T^n=P_{\mathrm{jc}}^n---\left(0.5\right)$

其中，j表示壁板特征点。 

（3）对试验模型机翼来说，刚度相对较弱，因其受气动载荷作用而发生弹性变形，其像点坐标满足极线约束条件： 

${\left(P_{k2}^n\right)}^T{\mathrm{EP}}_{k1}^n=0---\left(0.6\right)$

式中，k表示模型机翼特征点，为1#相机和2#相机的归一化像点的齐次坐标和

令$T_c^n=\left(\begin{array}{c}{t}_1\\ t_2\\ t_3\end{array}\right),$则$E=R_c^n\left(\begin{array}{ccc}0& {-t}_3& t_2\\ t_3& 0& {-t}_1\\ {-t}_2& t_1& 0\end{array}\right)$

步骤四、构建迭代方程： 

将约束条件(0.3)带入观测方程(0.2)，得到： 

$\left(\begin{array}{c}{M}_1^i(R_{\mathrm{mc}}^n{P}_{\mathrm{ic}}^0+T_{\mathrm{mc}}^n)=0\\ M_2^i(R_c^n(R_{\mathrm{mc}}^n{P}_{\mathrm{ic}}^0+T_{\mathrm{mc}}^n)+T_c^n)=0\end{array}\right)$

设函数$\left(\begin{array}{c}{F}_1=M_1^i(R_{\mathrm{mc}}^n{P}_{\mathrm{ic}}^0+T_{\mathrm{mc}}^n)=0\\ F_2=M_2^i(R_c^n(R_{\mathrm{mc}}^n{P}_{\mathrm{ic}}^0+T_{\mathrm{mc}}^n)+T_c^n)=0\end{array}\right)$

定义向量w=[a b c]^T的算子${\left|w\right|}_{\times }=\left(\begin{array}{ccc}0& -c& b\\ c& 0& -a\\ -b& a& 0\end{array}\right).$使用罗德里格斯变换公式R=(I+|w|_×)(I-|w|_×)^-1表示旋转矩阵，定义对应的参数为 $w_{\mathrm{mc}}=\left(\begin{array}{c}{a}_{\mathrm{mc}}\\ b_{\mathrm{mc}}\\ c_{\mathrm{mc}}\end{array}\right),$对应的参数为$w_c=\left(\begin{array}{c}{a}_c\\ b_c\\ c_c\end{array}\right).$将函数F₁和F₂进行一阶泰勒展开，得到误差公式： 

$\left(\begin{array}{c}\Delta F_1+M_1^i(-{\left|P_{\mathrm{ic}}^0\right|}_{\times }{\mathrm{\Delta w}}_{\mathrm{mc}}+{\mathrm{\Delta T}}_{\mathrm{mc}})+L_1=0\\ {\mathrm{\Delta F}}_2+M_2^i(-R_c^n{\left|P_{\mathrm{ic}}^0\right|}_{\times }{\mathrm{\Delta w}}_{\mathrm{mc}}+R_c^n{\mathrm{\Delta T}}_{\mathrm{mc}}-{|R_{\mathrm{mc}}^n{P}_{\mathrm{ic}}^0+T_{\mathrm{mc}}^n|}_{\times }{\mathrm{\Delta w}}_c+{\mathrm{\Delta T}}_c)+L_2=0\end{array}\right)$

使目标函数L₁和L₂的模为最小，可以得到法方程组，改写为矩阵形式： 

$\left(\begin{array}{cccc}{M}_1^i{J}_1& M_1^{i}I& 0& 0\\ M_2^i{J}_2& M_2^i{R}_c^n& M_2^i{J}_3& M_2^i\end{array}\right)\left(\begin{array}{c}{\mathrm{\Delta w}}_{\mathrm{mc}}\\ {\mathrm{\Delta T}}_{\mathrm{mc}}\\ {\mathrm{\Delta w}}_c\\ {\mathrm{\Delta T}}_c\end{array}\right)=\left(\begin{array}{c}{-\mathrm{\Delta F}}_1\\ -\Delta F_2\end{array}\right)---\left(0.7\right)$

上式对所有机身点i均成立，式中$J_3=-{|R_{\mathrm{mc}}^n{P}_{\mathrm{ic}}^0+T_{\mathrm{mc}}^n|}_{\times }.$

同理，对壁板特征点坐标，满足公式： 

$\left(\begin{array}{c}{F}_3=M_1^j(R^n{P}_{\mathrm{jc}}^0+T^n)=0\\ F_4=M_2^j(R_c^n(R^n{P}_{\mathrm{jc}}^0+T^n)+T_c^n)=0\end{array}\right).$定义Rⁿ的平差参数为$w_n=\left(\begin{array}{c}{a}_3\\ b_3\\ c_3\end{array}\right),$ΔTⁿ=ΔT_n，对上式进行一阶泰勒展开，得到误差公式： 

$\left(\begin{array}{c}{\mathrm{\Delta F}}_3+M_1^j(-{\left|P_{\mathrm{jc}}^0\right|}_{\times }\Delta w_n+{\mathrm{\Delta T}}_n)+L_3=0\\ {\mathrm{\Delta F}}_4+M_2^j(-R_c^n{\left|P_{\mathrm{jc}}^0\right|}_{\times }{\mathrm{\Delta w}}_n+R_c^n{\mathrm{\Delta T}}_n-{|R^n{P}_{\mathrm{jc}}^0+T^n|}_{\times }{\mathrm{\Delta w}}_c+{\mathrm{\Delta T}}_c)+L_4=0\end{array}\right)$

为使得目标函数L₃和L₄的模为最小，可以得到法方程组： 

$\left(\begin{array}{cccc}{M}_2^j{J}_5& M_2^j& M_2^j{J}_6& M_2^j{R}_c^n\\ 0& 0& M_1^j{J}_4& M_1^j\end{array}\right)\left(\begin{array}{c}{\mathrm{\Delta w}}_c\\ {\mathrm{\Delta T}}_c\\ {\mathrm{\Delta w}}_n\\ {\mathrm{\Delta T}}_n\end{array}\right)=\left(\begin{array}{c}-\Delta F_4\\ -\Delta F_3\end{array}\right)---\left(0.8\right)$

上式对所有壁板特征点j均成立，式中$J_6=-R_c^n{\left|P_{\mathrm{jc}}^0\right|}_{\times }.$

对模型机翼特征点，定义函数其中定义$P_{k1}^n={\left(\begin{array}{ccc}{x}_1& y_1& 1\end{array}\right)}^T,$$\left(\begin{array}{ccc}{t}_1& t_2& t_3\end{array}\right)={\left|T_c^n\right|}_{\times }=\left(\begin{array}{ccc}0& {-t}_3& t_2\\ t_3& 0& {-t}_1\\ -t_2& t_1& 0\end{array}\right),$函数f₅变形为 $f_5={\left(P_{k2}^n\right)}^T[x_1{R}_c^n{t}_1+y_1{R}_c^n{t}_2+R_c^n{t}_3]=0$

对函数f₅进行一阶泰勒展开，得到误差方程： 

${\mathrm{\Delta f}}_5+{\left(P_{k2}^n\right)}^T(J_7{\mathrm{\Delta w}}_c-R_c^n{\left|P_{k1}^n\right|}_{\times }{\mathrm{\Delta T}}_c)+l_k=0$

式中， 

$J_7=-x_1{\left|t_1\right|}_{\times }-y_1{\left|t_2\right|}_{\times }-{\left|t_3\right|}_{\times }=\left(\begin{array}{ccc}0& x_1{t}_2-y_1{t}_1& x_1{t}_3-t_1\\ -x_1{t}_2+y_1{t}_1& 0& y_1{t}_3-t_2\\ -x_1{t}_3+t_1& -y_1{t}_3+t_2& 0\end{array}\right)={\left(\begin{array}{c}{y}_1{t}_3-t_2\\ {-x}_1{t}_3+t_1\\ x_1{t}_2-y_1{t}_1\end{array}\right)}_{\times },$

改写为矩阵形式： 

$\left(\begin{array}{cc}{\left(P_{k2}^n\right)}^T{J}_7& -{\left(P_{k2}^n\right)}^T{R}_c^n{\left|P_{k1}^n\right|}_{\times }\end{array}\right)\left(\begin{array}{c}{\mathrm{\Delta w}}_c\\ {\mathrm{\Delta T}}_c\end{array}\right)=\left[{-\mathrm{\Delta f}}_5\right]---\left(0.9\right)$

上式对所有机翼特征点k均成立，则为使得目标函数$L_5=\left(\begin{array}{c}{l}_1\\ ·\\ ·\\ ·\\ l_k\end{array}\right)$模为最小，可以得到法方程组： 

$\left(\begin{array}{c}{\left(P_{12}^n\right)}^T{J_7}^1-{\left(P_{12}^n\right)}^T{R}_c^n{\left|P_{11}^n\right|}_{\times }\\ ·\\ ·\\ ·\\ {\left(P_{k2}^n\right)}^T{J_7}^k-{\left(P_{k2}^n\right)}^T{R}_c^n{\left|P_{k1}^n\right|}_{\times }\end{array}\right)\left(\begin{array}{c}{\mathrm{\Delta w}}_c\\ {\mathrm{\Delta T}}_c\end{array}\right)={\left(\begin{array}{c}-\Delta {f_5}^1\\ ·\\ ·\\ ·\\ -\Delta {f_5}^k\end{array}\right)}_{k\times 1}---\left(0.10\right)$

假定系统有i个模型机身特征点、j个试验段壁板点、k个模型机翼特征点，合并公式(0.7)、(0.8)和(0.10)，并假定所有的公式具有相同的权重，即可以得到含有18个平差变量的迭代方程： 

上式中${J_1}^i=-{\left|P_{\mathrm{ic}}^0\right|}_{\times },$${J_2}^i={-R}_c^n{\left|P_{\mathrm{ic}}^0\right|}_{\times },$${J_3}^i=-{|R_{\mathrm{mc}}^n{P}_{\mathrm{ic}}^0+T_{\mathrm{mc}}^n|}_{\times },$${J_4}^j=-{{|P}_{\mathrm{jc}}^0|}_{\times },$${J_5}^j=-{|R^n{P}_{\mathrm{jc}}^0+T^n|}_{\times },$${J_6}^j={-R}_c^n{\left|P_{\mathrm{jc}}^0\right|}_{\times },$${J_7}^k={\left(\begin{array}{c}{y_1}^k{t}_3-t_2\\ {{-x}_1}^k{t}_3+t_1\\ {x_1}^k{t}_2-{y_1}^k{t}_1\end{array}\right)}_{\times }.$

在实际工程应用的过程中，还可以根据风洞现场的实际试验条件，对上述方程进行简化。例如对于没有弹性体部件的导弹类试验模型来说，可以去掉上式中所有的和迭代方程同样可以成立；又例如在不考虑试验段壁板点约束的情况下，上述迭代方程的系数矩阵去掉最后的4j行并修改平差参数为${\left(\begin{array}{c}{\mathrm{\Delta w}}_{\mathrm{mc}}\\ {\mathrm{\Delta T}}_{\mathrm{mc}}\\ {\mathrm{\Delta w}}_c\\ \Delta T_c\end{array}\right)}_{12\times 1},$一样可以求解出相机的振动偏差，但此时就无法求解出模型的当前姿态。 

下面将上述迭代方程简写为： 

JΔD=ΔF   (0.11) 

步骤五、迭代求取方程组最优解： 

以系统零状态时相机之间的相对姿态（由系统标定得到）分解构成迭代初值$D_0=\left(\begin{array}{c}0\\ 0\\ w_c\\ T_c\\ 0\\ 0\end{array}\right),$以模型零状态下所有特征点的三维坐标和图像坐标作为以及图像帧n时刻机翼特征点像点归一化坐标构成偏差矩阵ΔF。利用公式(0.11)求得当前平差参数D的改正数ΔD=(J^T J)^-1J^TΔF，更新迭代参数D=D₀+ΔD。 

以新的迭代参数为初值，重新计算矩阵J和ΔF，计算新的平差参数D的改正数ΔD，如此反复迭代计算直至ΔF^T ΔF小于给定的阈值，迭代计算过程结束，此时得到的迭代参数D即为满足系统全部约束条件的最优解。 

步骤六、重构当前相机姿态和模型姿态。包括当前时刻相机新的位姿关系相机坐标系振动修正量Rⁿ、Tⁿ和模型姿态

将上述步骤计算得到的最优解D表示为：D=[a₁ b₁ c₁ t₁₁ t₁₂ t₁₃ a₂ b₂ c₂ t₂₁ t₂₂ t₂₃ a₃ b₃ c₃ t₃₁ t₃₂ t₃₃]^T，则系统当前时刻1#相机和2#相机之间的位姿关系为： 

$\left(\begin{array}{c}{R}_c^n=(I+S_2){(I-S_2)}^{-1}\\ T_c^n={\left(\begin{array}{ccc}{t}_{21}& t_{22}& t_{23}\end{array}\right)}^T\end{array}\right),$其中$S_2=\left(\begin{array}{ccc}0& {-c}_2& b_2\\ c_2& 0& {-a}_2\\ {-b}_2& a_2& 0\end{array}\right)$

相机坐标系振动修正量为： 

$\left(\begin{array}{c}{R}^n=(I+S_3){(I-S_3)}^{-1}\\ T^n={\left(\begin{array}{ccc}{t}_{31}& t_{32}& t_{33}\end{array}\right)}^T\end{array}\right),$其中$S_3=\left(\begin{array}{ccc}0& -c_3& b_3\\ c_3& 0& {-a}_3\\ -b_3& a_3& 0\end{array}\right)$

试验模型在试验段坐标系中的姿态为： 

$\left(\begin{array}{c}{R}_m^n=R^n(I+S_1){(I-S_1)}^{-1}\\ T_m^n=R^n\left(\begin{array}{c}{t}_{11}\\ t_{12}\\ t_{13}\end{array}\right)+T^n\end{array}\right),$其中$S_1=\left(\begin{array}{ccc}0& -c_1& b_1\\ c_1& 0& -a_1\\ {-b}_1& a_1& 0\end{array}\right)$ 。