CN112254743B - 一种基于星角距相减的星敏感器在轨标定方法 - Google Patents

Abstract

本发明涉及一种基于星角距相减的星敏感器在轨标定方法，属于天文导航领域。研究方法如下：基于场景物体成像到成像平面的物理过程建立针孔成像的相机模型；建立星敏感器角距模型；星敏感器相减方法；可观测性分析；星角距方法的改进；仿真实验。本发明通过改进的ADS（angular distance subtraction）算法，u ₀ 和v ₀ 比传统的AD（angular distance）算法精度分别提高了64.0%，21.7%，有效提高了主点标定精度。

Description

一种基于星角距相减的星敏感器在轨标定方法

技术领域

本本发明涉及天文导航领域，特别涉及一种标定方法，尤指一种基于星角距相减的星敏感器在轨标定方法。

背景技术

星敏感器是一种导航系统，通过对恒星的观测，获取载体的姿态信息。它是目前最精确的光学姿态传感器。由于其导航精度高、自主性强、无累积误差，受到航空航天工业的青睐。作为航天器的“眼睛”，星敏感器的精度直接决定了航天器的性能。然而，星敏感器是一种光学器件，其精度取决于成像质量和光学参数(包括焦距、主点和畸变)的精度。因此，标定是星敏感器的关键技术之一。

星敏感器在标定时需要进行实时标定以提高其导航精度，在所有标定参数中，主点的位置相比其他参数更容易受测量误差影响从而导致标定精度低，目前星敏感器在轨标定方法大多是以星间角距离作为标定参考，这种方法主点的可观测度较差，对主点的标定精度没有其他参数精度高。

发明内容

本发明涉及一种基于星角距相减的星敏感器在轨标定方法，解决了现有技术存在的上述问题。针对传统方法的不足，本发明提供了一种基于星角距相减的星敏感器在轨标定方法，在星敏感器标定中具有良好的应用效果，尤其适用于计算资源有限的星敏感器，针对星角距相减方法耗时长的问题，提出了改进星角距相减模型，与传统的星角距标定方法相比，改进ADS(angular distance subtraction)算法的u₀和v₀比AD(angular distance)法精度分别提高了64.0％，21.7％，有效提高了主点标定精度。

本发明的上述目的通过以下技术方案实现：

基于星角距相减的星敏感器在轨标定方法，包括如下步骤：

步骤1)、基于场景物体成像到成像平面的物理过程建立针孔成像的相机模型；

步骤2)、建立星角距相减模型；

步骤3)、可观测性分析；

步骤4)、星角距方法的改进。

步骤1)所述的基于场景物体成像到成像平面的物理过程建立针孔成像的相机模型是：(1)相机坐标系(O_c-X_cY_cZ_c)：以相机的光心为坐标原点，X轴和Y轴分别平行于图像坐标系的X轴和Y轴，相机的光轴为Z轴；

(2)物理图像坐标系(o’-xy)：以CCD图像平面与相机光轴的交点为坐标原点o’，X轴和Y轴分别平行于图像平面的两条垂直边，图像坐标系是用物理单位(例如毫米)表示像素在图像中的位置；

(3)像素坐标系(o-uv):以CCD图像平面的左上角顶点为原点，X轴和Y轴分别平行于图像坐标系的X轴和Y轴，像素坐标系就是以像素为单位的图像坐标系；

假设w＝[X，Y，Z]^T相机坐标系的任意星单位向量，其在针孔模型下理想的图像投影为p＝[x，y]^T，w和p之间的透视投影关系可以表示为：

其中,[u,v,1]^T是点p的齐次坐标，f_u和f_v分别为u轴和v轴方向的像素焦距；实际上镜头都存在不同大小的畸变，在考虑畸变的情况下，可以使用下述公式来描述相机非线性模型：

其中，(u，v)是等式中的无失真坐标，(u_d，v_d)是非线性模型下的图像坐标，即考虑镜头几何畸变时的图像坐标。δ_u(u，v)和δ_v(u，v)分别是u和v方向的畸变。相机的畸变主要分为径向畸变、偏心畸变以及薄棱镜畸变三种类型，由于径向畸变所带来的影响是最大的，并且高阶畸变可能导致数值不稳定，因此这里我们只考虑径向畸变的一阶和二阶，畸变方程为：

其中k₁，k₂是径向畸变系数。

步骤2)所述的建立星角距相减模型是：

目前星敏感器标定大多使用角距为标定参考,星敏感器的焦距为f，w和v分别为恒星在星敏感器坐标系中和天球坐标系中的方向矢量,恒星i在星敏感器成像平面坐标系下投影点中心坐标为(x_i，y_i)，则

其中，(x₀，y₀)主点坐标，α_i和δ_i分别表示第i颗星的赤经和赤纬；假设有i，j两颗星，根据星角距正交变换不变原理，在不考虑畸变和噪声的条件下，恒星i，j在星敏感器坐标系中的方向矢量w_i，w_j的夹角与对应的天球坐标系中的位置矢量v_i，v_j的夹角θ_ij相等，可以表示为：

cosθ_ij＝w_i ^Tw_j＝v_i ^Tv_j (7)

把公式(5)代入公式(7)中，得到

其中，

假设有i，j，k三颗星，利用公式(8)分别计算v_i ^Tv_j和v_j ^Tv_k，再将两式相减，用S表示，得到：

由于焦距f比CCD尺寸大很多，因此

D_iD_j≈D_jD_k (11)

将式(11)近似为D_aD_b，a，b为任意两颗星，a≠b；公式(10)可以写为：

即将分子的f²项消除，减小焦距的影响，放大对主点的计算，从而提高主点的可观测度，提高标定精度；使用扩展卡尔曼滤波的方法进行标定，对星点图像序列进行循环迭代，则可得到标定后的参数，状态方程为：

x_k＝I_2×2·x_k-1 (13)

其中x_k为需要标定的主点参数(x₀，y₀)，k-1和k分别代表第k-1和第k幅图像，I_2X2为单位矩阵，测量方程是：

z_k＝h(x_k)+n_c (14)

其中，z_k为星角距相减形成的矩阵，由天球坐标系中的位置矢量计算得到，h(x_k)为利用星敏感器标定模型和图像点求解星角距相减的过程，n_c是由噪声引起的测量误差，EKF预测方程为：

其中P_k ^-为k时刻的先验估计协方差，Q为系统过程的协方差矩阵，EKF更新方程为：

其中R是观测噪声的协方差矩阵，H_k是雅可比矩阵。

步骤3)所述的可观测性分析是：

可观测性可以反映状态可估计性的能力，是评价系统可行性的指标，即在不同的模型下，相同的输入偏差可能导致不同的输出偏差，如果输出偏差的幅度较大，即在相同的输出偏差下，输入偏差更小，则可观测性更好，系统更可行，反之亦然；根据可观测性的定义，我们可以得到：

δz_k＝H_kδx (20)

其中δx是输入偏差，δz_k输出偏差，H_k是雅可比矩阵，对雅克比矩阵进行可观测性分析，利用可观测矩阵的奇异值分解，公式(20)可表示为：

δz_k＝P_k∑_kQ_kδx (21)

其中P_k和Q_k分别是左奇异向量和右奇异向量的正交矩阵，∑_k为2×N对角矩阵，对角元素是非零奇异值σ_i(i＝1～2)；由于P_k和Q_k是正交矩阵，可以计算：

其中‖δx‖₂和‖δz_k‖₂分别为δx和δz_k的2范数，输出偏差的下确界为：

其中，σ_min是可观测矩阵的最小奇异值，σ_min越大，输出偏差的最小值越大，可观测性越好；根据以上分析，我们对式(16)中的雅克比矩阵H_k进行奇异值分解，计算每一帧的最小奇异值σ_min，对传统星角距方法和星角距相减方法进行可观测性分析比较；相比而言，星角距相减方法的最小奇异值比传统方法的高很多，这表示星角距相减方法的可观测性更好，但是由于使用角距相减组合过多不可避免会增加计算量，导致星敏感器工作效率降低，无法满足实时性要求，因此我们在此基础上对模型进行改进，在不降低标定精度的前提下减少星角距相减的数量，提高标定效率。

步骤4)所述的星角距方法的改进是：

为了提高星角距相减方法的效率，我们提出四种改进方法。假设一帧星点图有N颗星点，根据公式(7)，计算任意两颗星的角距离，可以得到一个N阶星角距对称方阵，我们只取上半部分三角形矩阵，给出如下四种相减方式方法1：每一行横向依次相减，即a₁₂-a₁₃,a₁₃-a₁₄,a₁₄-a₁₅……每行最后一个再和下一行第一个相减，即a_1N-a₂₃,a_2N-a₃₄……

方法2：同方法1，先每行横向依次相减，即a₁₂-a₁₃,a₁₃-a₁₄,a₁₄-a₁₅,……再每列依次相减，即a₁₃-a₂₃,a₁₄-a₂₄,a₂₄-a₃₄……

方法3：在方法2的基础上，增加对角线相减，即a₁₂-a₂₃,a₂₃-a₃₄,a₃₄-a₄₅……

方法4：在方法3的基础上，除了对角线相减外，增加所有斜向相减，即a₁₃-a₂₄,a₁₄-a₂₅,a₁₅-a₂₆……

从使用次数和均匀度综合来看，改进2较好。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，本发明的示意性实例及其说明用于解释本发明，并不构成对本发明的不当限定。

图1为摘要附图；

图2为相机坐标系与物理图像坐标系的关系；

图3为图像坐标系与物理图像坐标系的关系；

图4为星敏感器角距模型；

图5为传统星角距和星角距相减方法的最小奇异值比较；

图6为星角距矩阵；

图7为噪声标准差为0.5像素实验的姿态残差图。

具体实施方式

下面结合附图进一步说明本发明的详细内容及其具体实施方式。

1.建立星角距相减模型：

目前星敏感器标定大多使用角距为标定参考,星敏感器的焦距为f，w和v分别为恒星在星敏感器坐标系中和天球坐标系中的方向矢量,恒星i在星敏感器成像平面坐标系下投影点中心坐标为(x_i，y_i)，则

其中，(x₀，y₀)主点坐标，α_i和δ_i分别表示第i颗星的赤经和赤纬；假设有i，j两颗星，根据星角距正交变换不变原理，在不考虑畸变和噪声的条件下，恒星i，j在星敏感器坐标系中的方向矢量w_i，w_j的夹角与对应的天球坐标系中的位置矢量v_i，v_j的夹角θ_ij相等，可以表示为：

cosθ_ij＝w_i ^Tw_j＝v_i ^Tv_j (3)

把公式(5)代入公式(7)中，得到

其中，

假设有i，j，k三颗星，利用公式(4)分别计算v_i ^Tv_j和v_j ^Tv_k，再将两式相减，用S表示，得到：

由于焦距f比CCD尺寸大很多，因此

D_iD_j≈D_jD_k(7)

将式(11)近似为D_aD_b，a,b为任意两颗星，a≠b。公式(6)可以写为：

即将分子的f ²项消除，减小焦距的影响，放大对主点的计算，从而提高主点的可观测度，提高标定精度。

2.根据上述方法进行标定

使用扩展卡尔曼滤波，对星点图像序列进行循环迭代，则可得到标定后的参数，状态方程为：

x_k＝I_2×2·x_k-1 (9)

其中x_k为需要标定的主点参数(x₀,y₀)，k-1和k分别代表第k-1和第k幅图像，I_2X2为单位矩阵，测量方程是：

z_k＝h(x_k)+n_c (10)

其中，z_k为星角距相减形成的矩阵，由天球坐标系中的位置矢量计算得到，h(x_k)为利用星敏感器标定模型和图像点求解星角距相减的过程，n_c是由噪声引起的测量误差，EKF预测方程为：

其中P_k ^-为k时刻的先验估计协方差，Q为系统过程的协方差矩阵，EKF更新方程为：

其中R是观测噪声的协方差矩阵，H_k是雅可比矩阵。

3.可观测性分析

可观测性可以反映状态可估计性的能力，是评价系统可行性的指标，即在不同的模型下，相同的输入偏差可能导致不同的输出偏差。如果输出偏差的幅度较大，即在相同的输出偏差下，输入偏差更小，则可观测性更好，系统更可行，反之亦然。根据可观测性的定义，我们可以得到：

δz_k＝H_kδx (16)

其中δx是输入偏差，δz_k输出偏差，是雅可比矩阵，对雅克比矩阵进行可观测性分析，利用可观测矩阵的奇异值分解，公式(15)可表示为：

δz_k＝P_k∑_kQ_kδx (17)

其中P_k和Q_k分别是左奇异向量和右奇异向量的正交矩阵，∑_k为2×N对角矩阵，对角元素是非零奇异值σ_i(i＝1～2)。由于P_k和Q_k是正交矩阵，可以计算：

其中‖δx‖₂和‖δz_k‖₂分别为δx和δz_k的2范数。输出偏差的下确界为：

其中，σ_min是可观测矩阵的最小奇异值，σ_min越大，输出偏差的最小值越大，可观测性越好；根据以上分析，我们对式(11)中的雅克比矩阵H_k进行奇异值分解，计算每一帧的最小奇异值σ_min，对传统星角距方法和星角距相减方法进行可观测性分析比较；相比而言，星角距相减方法的最小奇异值比传统方法的高很多，这表示星角距相减方法的可观测性更好，但是由于使用角距相减组合过多不可避免会增加计算量，导致星敏感器工作效率降低，无法满足实时性要求，因此我们在此基础上对模型进行改进，在不降低标定精度的前提下减少星角距相减的数量，提高标定效率。

4.对星角距方法进行改进

为了提高星角距相减方法的效率，我们提出四种改进方法。假设一帧星点图有N颗星点，根据公式(3)，计算任意两颗星的角距离，可以得到一个N阶星角距对称方阵，我们只取上半部分三角形矩阵，给出如下四种相减方式，

方法1：每一行横向依次相减，即a₁₂-a₁₃,a₁₃-a₁₄,a₁₄-a₁₅……每行最后一个再和下一行第一个相减，即a_1N-a₂₃,a_2N-a₃₄……

方法2：同方法1，先每行横向依次相减，即a₁₂-a₂₃,a₁₃-a₁₄,a₁₄-a₁₅,……再每列依次相减，即a₁₃-a₂₃,a₁₄-a₂₄,a₂₄-a₃₄……

方法3：在方法2的基础上，增加对角线相减，即a₁₂-a₂₃,a₂₃-a₃₄,a₃₄-a₄₅……

方法4：在方法3的基础上，除了对角线相减外，增加所有斜向相减，即a₁₃-a₂₄,a₁₄-a₂₅,a₁₅-a₂₆……

表1：

对以上四种改进方法进行均匀度分析，以7颗星为例，每颗星的使用次数和使用率如表1所示；改进1方法每颗星的使用率最均匀，但是使用次数较少，可能会影响标定精度；改进3与改进4每颗星使用次数多，但是均匀度较差；因此从使用次数和均匀度综合来看，改进2较好。

5.仿真实验

采用19.14°×11.18°视场、1920×1080像素阵列星敏感器，以2Hz更新率进行模拟，模拟数据由三组数据组成：惯性坐标系中的3D星矢量、图像坐标系中相应的2D星坐标和具有正态分布噪声的2D星坐标；为了避免其他参数对主点精度的影响，我们将焦距和畸变参数都设为标准值，只对主点进行标定，参数和标准值见表2：

为了充分评价标定方法的性能，提出姿态残差作为评价标准：利用标定所使用的带噪声的仿真数据及主点标定结果计算姿态矩阵和光轴指向，与设定的标准参数计算出的姿态矩阵和光轴指向进行比较，计算每幅图像求得的两个光轴指向的夹角得到的即为姿态残差，最终评价指标是最后100幅图像中姿态残差的均值(μ_Att)；测试1传统方法与星角距相减方法比较，我们通过设置星敏感器的极限视觉星等，得到恒星数目不同的数据进行实验，AD方法的内参数和畸变参数同样使用表2的标准值，只对主点进行标定，三次实验的平均恒星数分别为7.7、13.6和19.1，极限视星等分别为4.6、5和5.5,2D星坐标添加均值为0，标准差为0.5的正态分布噪声，结果见表3：

表3中，Δu₀和Δv₀表示标定的主点结果与表1中标准值的差值，μ_Att是第四节提出的评价标准，即最后100幅图像中姿态残差的均值，T是处理一幅图像所用的时间，根据表3的结果，星角距相减方法的主点标定精度明显好于传统方法，但是和设想的一样，星角距相减方法耗时长，尤其是在星点数量较多时，这会严重影响星敏感器标定效率，无法满足实时性要求；在测试2中我们采用四种改进方案进行实验，为了验证3C节中对改进方法使用次数和均匀度的分析，进行如下实验。使用每帧平均恒星数为7.7颗星的数据进行实验，添加均值为0，标准差为0.5的正态分布噪声，实验结果如表4：

从表4可以看出，改进1虽然星点使用率比较均匀，但是星角距相减组合数较少，一定程度上导致标定效果不好，改进3和4的星点使用均匀度较差，这会导致标定结果不稳定，且改进4使用角距相减组合数量较多，虽然标定效果比改进1好，但是耗时明显增加，相比之下，改进2的姿态残差最小，精度比改进1提高了44.3％，比ADS方法提高了49.2％，同时耗时比ADS方法减少了12.9％，因此改进2标定效果更好，这与我们在3C节中的分析一致；在测试3中，我们将改进2与传统方法进行比较，使用每帧平均恒星数为7.7颗星，进行三组实验，每组实验噪声标准差分别为0.2像素，0.5像素，0.8像素。实验结果如表5所示：

从表5中可以看出，噪声标准差为0.2像素时改进方法的Δu₀比传统方法降低了约50.5％，但是Δv₀比传统方法精度提高了约61.1％，从姿态残差也可以看出改进方法对主点的标定精度有所提升；噪声标准差为0.5像素时，改进方法的Δu₀和Δv₀均好于传统方法；噪声标准差为0.8像素时改进方法的Δv₀与传统方法相差不大，但是Δu₀有很大幅度的提升；在测试4中，我们将改进2与传统方法全参数比较，为了充分分析IADS2方法的性能，我们对所有参数进行标定实验，实验条件同测试2，校准结果见表6:

	Δf(mm)	Δk1	Δk2	u0(pixel)	v0(pixel)	T(ms)
							AD	0.00002	-0.0007	0.0565	-2.060	2.639	126.446
改2	0.00058	-0.0087	0.3048	-1.329	1.926	88.561

我们发现，IADS2算法在其他参数标定效果没有AD方法好，但是在其他参数标定精度不高的前提下，主点的标定精度依然要好于传统方法，因此在标定主点前，可以先使用AD法标定其他参数，实现主点高精度标定。

Claims (1)

1.一种基于星角距相减的星敏感器在轨标定方法，其特征在于：包括如下步骤：

步骤1)、基于场景物体成像到成像平面的物理过程建立针孔成像的相机模型；

步骤2)、建立星敏感器角距相减模型；

步骤3)、可观测性分析；

步骤4)、星角距方法的改进；

具体内容如下：

其中，步骤1)所述的基于场景物体成像到成像平面的物理过程建立针孔成像的相机模型是：

(1)相机坐标系(O_c-X_cY_cZ_c)：以相机的光心为坐标原点，X轴和Y轴分别平行于图像坐标系的X轴和Y轴，相机的光轴为Z轴；

(2)物理图像坐标系(o’-xy)：以CCD图像平面与相机光轴的交点为坐标原点o’，X轴和Y轴分别平行于图像平面的两条垂直边，图像坐标系是用物理单位(例如毫米)表示像素在图像中的位置；

(3)像素坐标系(o-uv):以CCD图像平面的左上角顶点为原点，X轴和Y轴分别平行于图像坐标系的X轴和Y轴，像素坐标系就是以像素为单位的图像坐标系；

假设w＝[X，Y，Z]^T相机坐标系的任意星单位向量，其在针孔模型下理想的图像投影为p＝[x，y]^T，w和p之间的透视投影关系可以表示为：

其中,[u,v,1]^T是点p的齐次坐标，f_u和f_v分别为u轴和v轴方向的像素焦距；实际上镜头都存在不同大小的畸变，在考虑畸变的情况下，可以使用下述公式来描述相机非线性模型：

其中，(u，v)是等式中的无失真坐标，(u_d，v_d)是非线性模型下的图像坐标，即考虑镜头几何畸变时的图像坐标，δ_u(u，v)和δ_v(u，v)分别是u和v方向的畸变，相机的畸变主要分为径向畸变、偏心畸变以及薄棱镜畸变三种类型，由于径向畸变所带来的影响是最大的，并且高阶畸变可能导致数值不稳定，因此这里我们只考虑径向畸变的一阶和二阶，畸变方程为：

其中k₁，k₂是径向畸变系数；

其中，步骤2)所述的星敏感器角距相减方法模型是：

目前星敏感器标定大多使用角距为标定参考,星敏感器的焦距为f，w和v分别为恒星在星敏感器坐标系中和天球坐标系中的方向矢量,恒星i在星敏感器成像平面坐标系下投影点中心坐标为(x_i，y_i)，则

其中，(x₀，y₀)主点坐标，α_i和δ_i分别表示第i颗星的赤经和赤纬；假设有i，j两颗星，根据星角距正交变换不变原理，在不考虑畸变和噪声的条件下，恒星i，j在星敏感器坐标系中的方向矢量w_i，w_j的夹角与对应的天球坐标系中的位置矢量v_i，v_j的夹角θ_ij相等，可以表示为：

cosθ_ij＝w_i ^Tw_j＝v_i ^Tv_j (7)

把公式(5)代入公式(7)中，得到

其中，

假设有i，j，k三颗星，利用公式(8)分别计算v_i ^Tv_j和v_j ^Tv_k，再将两式相减，用S表示，得到：

由于焦距f比CCD尺寸大很多，因此

D_iD_j≈D_jD_k(11)

将式(11)近似为D_aD_b；a，b为任意两颗星，a≠b；公式(10)可以写为：

即将分子的f ²项消除，减小焦距的影响，放大对主点的计算，从而提高主点的可观测度，提高标定精度；使用扩展卡尔曼滤波的方法进行标定，对星点图像序列进行循环迭代，则可得到标定后的参数，状态方程为：

x_k＝I_2×2·x_k-1 (13)

其中x_k为需要标定的主点参数(x₀，y₀)，k-1和k分别代表第k-1和第k幅图像，I_2X2为单位矩阵，测量方程是：

z_k＝h(x_k)+n_c (14)

其中，z_k为星角距相减形成的矩阵，由天球坐标系中的位置矢量计算得到，h(x_k)为利用星敏感器标定模型和图像点求解星角距相减的过程，n_c是由噪声引起的测量误差，EKF预测方程为：

其中P_k ^-为k时刻的先验估计协方差，Q为系统过程的协方差矩阵，EKF更新方程为：

其中R是观测噪声的协方差矩阵，H_k是雅可比矩阵；

其中，步骤3)所述的可观测性分析是：

可观测性可以反映状态可估计性的能力，是评价系统可行性的指标，即在不同的模型下，相同的输入偏差可能导致不同的输出偏差，如果输出偏差的幅度较大，即在相同的输出偏差下，输入偏差更小，则可观测性更好，系统更可行，反之亦然；根据可观测性的定义，我们可以得到：

δz_k＝H_kδx (20)

其中δx是输入偏差，δz_k输出偏差，H_k是雅可比矩阵，对雅克比矩阵进行可观测性分析，利用可观测矩阵的奇异值分解，公式(20)可表示为：

δz_k＝P_k∑_k Q_kδx (21)

其中P_k和Q_k分别是左奇异向量和右奇异向量的正交矩阵，∑_k为2×N对角矩阵，对角元素是非零奇异值σ_i(i＝1～2)；由于P_k和Q_k是正交矩阵，可以计算：

其中‖δx‖₂和‖δz_k‖₂分别为δx和δz_k的2范数，输出偏差的下确界为：

其中，σ_min是可观测矩阵的最小奇异值，σ_min越大，输出偏差的最小值越大，可观测性越好；根据以上分析，我们对式(16)中的雅克比矩阵H_k进行奇异值分解，计算每一帧的最小奇异值σ_min，对传统星角距方法和星角距相减方法进行可观测性分析比较；相比而言，星角距相减方法的最小奇异值比传统方法的高很多，这表示星角距相减方法的可观测性更好，但是由于使用角距相减组合过多不可避免会增加计算量，导致星敏感器工作效率降低，无法满足实时性要求，因此我们在此基础上对模型进行改进，在不降低标定精度的前提下减少星角距相减的数量，提高标定效率；

其中，步骤4)所述的星角距方法的改进是：

假设一帧星点图有N颗星点，根据公式(7)，计算任意两颗星的角距离，可以得到一个N阶星角距对称方阵，我们只取上半部分三角形矩阵，给出如下四种相减方式：

方法1：每一行横向依次相减，即a₁₂-a₁₃,a₁₃-a₁₄,a₁₄-a₁₅……每行最后一个再和下一行第一个相减，即a_1N-a₂₃,a_2N-a₃₄……

方法2：同方法1，先每行横向依次相减，即a₁₂-a₁₃,a₁₃-a₁₄,a₁₄-a₁₅,……再每列依次相减，即a₁₃-a₂₃,a₁₄-a₂₄,a₂₄-a₃₄……

方法3：在方法2的基础上，增加对角线相减，即a₁₂-a₂₃,a₂₃-a₃₄,a₃₄-a₄₅……

方法4：在方法3的基础上，除了对角线相减外，增加所有斜向相减，即a₁₃-a₂₄,a₁₄-a₂₅,a₁₅-a₂₆……。