CN102901459B - 一种基于准直光路的角秒级三维变形测量装置与方法 - Google Patents

Abstract

本发明公开了一种基于准直光路的角秒级三维变形测量装置及方法。装置由发射模块、接收模块和处理模块构成，发射模块由光源、光阑和发射光学系统构成，光源和光阑的中心都在发射光学系统的光轴上；光阑的透光部分采用直线结构；接收模块由接收光学系统和面阵探测器构成，面阵探测器位于接收光学系统的焦平面上，其中心位于接收光学系统的光轴上；处理模块中安装了变形解算软件。测量方法是发射光学系统将透过光阑的光束准直发射到接收模块，接收光学系统将光束会聚到面阵探测器上，面阵探测器将光束转换成光阑图像，变形解算软件根据光阑图像计算三维变形。采用本发明能同时测量方位角、俯仰角和横滚角变形，精度高，能实现角秒级三维变形测量。

Description

一种基于准直光路的角秒级三维变形测量装置与方法

技术领域

本发明属于光学变形测量技术领域，特别涉及一种基于准直光路的角秒级三维变形测量装置与方法。

背景技术

光学变形测量技术可广泛用于大型运动平台变形测量和高精度姿态基准传递等方面，具有测量精度高、直接测量、稳定性好等优点。目前的光学变形测量大多采用准直光路方法实现。

中国专利CN97251378.7、美国专利US7227627、《中国惯性技术学报》杂志2006年第3期中“一种舰船变形测量方法的研究与应用”和《长春理工大学学报》杂志2006年第3期中“车载平台变形的激光自准直测量方法研究”论文公开了一类基于自准直光路、平面反射镜和面阵CCD（Charge Coupled Device，电荷耦合器件）的二维变形测量装置，该装置如图1所示，主要由光源1、光阑2、分光镜3、面阵CCD探测器4、光学系统5和平面反射镜6构成。该装置可以实时测量两个载体之间的相对变形，为了方便起见将其中一个载体称为待测物，另外一个载体称为参考基准。可将该测量装置的平面反射镜6与待测物固连，其余部分与参考基准固连，反之亦然。这种测量装置的测量光学变形的方法是利用光源1照射光阑2，透过光阑2的光部分透过分光镜3，然后通过光学系统5对光阑2进行成像，所成的光阑像经过平面反射镜6、分光镜3和光学系统5后成像在面阵CCD探测器4上，根据面阵CCD探测器的光阑图像位置解算变形。该装置可以测量与光轴垂直方向的二维变形，不能测量绕光轴方向的旋转角变形。该装置的二维变形测量精度可以达到角秒量级。

美国专利US7697127公开了一种测量绕光轴旋转角的装置，该装置如图2所示，主要由光源1、光阑2和面阵探测器8构成。该装置的光阑2采用周期结构的光阑，由光阑2和面阵探测器8实现绕光轴7旋转角测量，具体测量方法：在光源1照射下光阑2直接投影到面阵探测器8上，利用图像处理方法处理面阵探测器8探测光阑图像得到旋转角变化。该装置可以测量绕光轴7的旋转角，但不能测量与光轴垂直方向的角度变化。由于没有准直光路和直接投影的图像质量限制，该装置的旋转角测量精度难以达到角秒量级。

中国专利CN200410010856.3、《计量技术》杂志2005年第4期中“干涉法在船体角度变形测量中的应用”、《光学技术》杂志2005年第5期中“船体三维角度变形的自准直干涉测量方法”和《飞行器测控学报》2009年第3期中“基于莫尔条纹测量法设计的船体变形测量系统”论文公开了一类干涉类型的三维变形测量装置，该装置组成与图1所示装置类似，如图3所示，由第一光源1、光阑2、第一分光棱镜3、第二分光棱镜12、第三分光棱镜7、光学系统5、反射镜6、第一面阵CCD探测器41、第二光源9、自准直光栅10、瞄准光栅11和第二面阵CCD探测器42构成。该装置的测量方法是：利用第一光源1、光阑2、第一分光棱镜3、第二分光棱镜12、第三分光棱镜7、光学系统5、反射镜6、第一面阵CCD探测器41构成一套自准直光路，采用与图1类似的方法测量与光轴垂直方向的二维变形角；利用第二光源9、自准直光栅10、第二分光棱镜12、第三分光棱镜7、光学系统5、反射镜6、瞄准光栅11和第二面阵CCD探测器42构成一套光栅干涉光路，采用莫尔条纹方法测量绕光轴方向的旋转角变形。该方法测量的与光轴垂直方向的二维变形（方位和俯仰方向变形）精度较高，能达到角秒量级；绕光轴方向的变形（横滚方向变形）测量精度相对较低，难以达到角秒量级。该装置采用部件较多且结构较复杂，具有制作难度大、稳定性相对较低和成本高等特点。

《光学精密工程》杂志2008年第11期中“船体变形测量的新技术”和《长春理工大学学报》杂志2006年第3期中“基于测量船角变形光电测量系统的研究”论文介绍了一种采用大钢管的三维变形测量装置，该装置的主要结构如图4所示，主要由第一发射光管13、第二发射光管14、第一接收光管19、第二接收光管18、大钢管16、钢管安装座24、第一测扭光管20、第一反射镜15、第二测扭光管21和第二反射镜17构成，用于测量待测物22与参考基准23之间的三维变形。其中，第一发射光管13、第二发射光管14、第一测扭光管20、第二测扭光管21和第二反射镜17都与待测物22固连；第一接收光管19、第二接收光管18、第二测扭光管21都与参考基准23固连；第一反射镜15与大钢管16的一端固连，第二反射镜17与大钢管16的另一端固连；大钢管16采用特殊材料设计加工，具有温度稳定性好和横滚方向变形小等特点；大钢管16由钢管安装座24安装，钢管安装座24能够支撑大钢管16并且不会使大钢管16在横滚方向产生变形。在第一接收光管19、第二接收光管18、第一测扭光管20和第二测扭光管21内部都分别安装了线阵CCD探测器。该装置通过第一发射光管13和第一接收光管19构成一套自准直光路，测量待测物22与参考基准之间的俯仰方向变形；通过第二发射光管14和第二接收光管18构成一套自准直光路，测量待测物22与参考基准之间的方位方向变形；通过第一测扭光管20和第一反射镜15测量待测物22与大钢管16一端之间在横滚方向的变形，通过第二测扭光管21和第二反射镜17测量参考基准23与大钢管16另一端在横滚方向的变形，由此计算待测物22与参考基准23在横滚方向的变形。该装置能够实时测量待测物22与参考基准23之间的三维变形，方位、俯仰和横滚方向变形的测量精度都可以达到角秒量级，但该装置具有结构复杂、成本高、加工难度大和占用空间大等不足。

如何既能实时测量待测物和参考基准之间的三维变形（方位、俯仰和横滚方向变形），又具有测量精度高、结构简单、占用空间小、体积小和成本低等特点，是本领域技术人员极为关注的技术问题。

发明内容

本发明要解决的技术问题是提出一种基于准直光路的角秒级三维变形测量装置与方法，可以同时测量三个方向（方位、俯仰和横滚方向）变形，并且具有测量精度高、结构简单、占用空间小和成本低等优点。

本发明三维变形测量装置由发射模块、接收模块和处理模块三大部分构成。发射模块与参考基准固连，接收模块与待测物固连。发射模块与接收模块之间能够通视，发射模块发射的光束直接进入接收模块。该装置能实时测量待测物与参考基准的三维变形。

发射模块由光源、光阑和发射光学系统构成。采用光学领域常用的光学固定座和调整架将光源、光阑和发射光学系统固定在参考基准上。光阑位于光源和发射光学系统之间，且位于发射光学系统的焦平面上，光源和光阑的中心都在发射光学系统的光轴上，由此构成发射准直光路。光源照射光阑，发射光学系统将透过光阑的光束准直发射到接收模块中。光源可以是半导体激光器或LED或灯泡等。发射光学系统可以采用双胶合透镜或双分离透镜或多片透镜组等类型光学系统。光阑是在玻璃等透明材料上利用光刻等工艺加工而成的，加工后的光阑2有些部分透光，有些部分不透光。光阑2的透光部分采用L_C（L_C为大于1的常数）条直线结构，各直线呈放射状分布，直线结构为透光部分，其余部分不透光。两直线结构光阑（如图6(a)所示）的透光部分由两条直线构成，两直线可互相垂直或不垂直。每条直线由两个线段构成，其中线段L1a、L1b构成直线L1，线段L2a、L2b构成直线L2。采用多直线结构光阑可以提高变形测量精度，三条直线结构光阑如6(b)所示，四条直线结构光阑如6(c)所示，四条以上的直线结构光阑依次类推，直线越多，测量精度越高，但横滚方向的形变测量范围会减小。

接收模块由接收光学系统和面阵探测器构成。采用光学领域常用的光学固定座和调整架将接收光学系统和面阵探测器固定在待测物上。面阵探测器位于接收光学系统的焦平面上，且其中心位于接收光学系统的光轴上，由此构成接收准直光路。接收光学系统可采用双胶合透镜或双分离透镜或多片透镜组等类型光学系统，面阵探测器可采用灰度或彩色的面阵CCD或CMOS等类型探测器。发射模块发射的准直光束经过接收光学系统会聚后成像在面阵探测器上，面阵探测器探测的是光阑的直线结构光阑图像。

处理模块可采用工控机、嵌入式计算机、台式计算机或笔记本电脑等设备，在处理模块中安装了变形解算软件。处理模块通过千兆以太网、高速USB或1394等接口与面阵探测器连接，实时采集面阵探测器探测的光阑图像，变形解算软件计算三维变形。

按图5所示方式定义坐标系，即发射模块光轴方向为Y轴，与Y轴垂直且朝上的方向为Z轴，X轴与Y轴、Z轴垂直。当参考基准和待测物之间发生变形时，发射模块与接收模块之间会发生整体变形。相对于发射模块，接收模块整体绕Z轴旋转产生方位角变形，绕Y轴旋转产生横滚角变形，绕X轴旋转产生俯仰角变形。当发射模块的光轴与接收模块的光轴平行且发射模块中光阑的水平直线与接收模块的面阵探测器的水平中心线平行时，发射模块与接收模块之间的方位、俯仰和横滚方向变形都为零。此时，面阵探测器探测的光阑图像中心与面阵探测器中心重合，光阑图像中直线结构的某一直线（称为水平直线）与面阵探测器的水平中心线平行，如图7所示（图中以两直线结构光阑为例，与面阵探测器的水平中心线平行的直线为水平直线）。

当发射模块与接收模块之间存在方位和俯仰方向变形时，面阵探测器探测的光阑图像发生平移，光阑图像中心偏离面阵探测器中心，如图8所示（图中以两直线结构光阑为例）。将光阑图像中心与面阵探测器中心在X轴方向的偏移量记为ΔX，Z轴方向的偏移量记为ΔZ；将接收模块光学系统的焦距记为f。根据ΔX、ΔZ和f可以计算方位方向变形ψ、俯仰方向变形θ。当发射模块与接收模块之间存在横滚方向变形时，面阵探测器的光阑图像发生旋转，如图9所示（图中两直线结构光阑为例）。光阑图像中直线方向与面阵探测器中心线存在夹角该角度为横滚方向变形。

采用本发明所述装置进行变形测量的具体方法如下：

第一步，发射模块的光源照射光阑，发射光学系统将透过光阑的光束准直发射到接收模块。

第二步，接收模块的接收光学系统接收发射模块发射的光束，并将其会聚到面阵探测器上。

第三步，面阵探测器将接收的光束转换成光阑图像，此光阑图像为数字化的彩色图像或灰度图像。

第四步，处理模块从面阵探测器接收光阑图像。

第五步，处理模块中变形解算软件根据光阑图像按以下步骤计算三维变形：

5.1若面阵探测器是彩色类型探测器，则将彩色图像转换成灰度图像，可参照RafaelC.Gonzalez等编著、阮秋琦等翻译、电子工业出版社2005年出版的《数字图像处理（MATLAB版）》中151～160页介绍的方法进行转换；若面阵探测器是灰度类型探测器，则直接进入步骤5.2。

5.2对光阑图像进行图像校正，可采用R.Y.Tsai在1986年IEEE Proceedings ofConference on Computer Vision and Pattern Recognition发表的An Efficient and AccurateCamera Calibration Technique For 3D Machine Vision论文介绍的传统标定或张正友在2000年第11期IEEE Transactions on Pattern Analysis and Machine Intelligence发表的A FlexibleNew Technique for Camera Calibration论文介绍的自标定方法对光阑图像进行校正。发射光学系统畸变、接收光学系统畸变和面阵探测器畸变等会导致光阑图像产生畸变，对光阑图像的畸变进行校正可提高变形测量精度。

5.3对校正后的光阑图像进行亮度分析，并滤除背景噪声。对光阑图像的所有像素根据亮度进行排序，根据光阑的直线结构的直线数量L_C选取L3个（L3通常为L_C的10～100倍）最亮像素，根据选择的最亮像素计算平均亮度，将该平均亮度作为光阑图像的最大亮度I_MAX。根据最大亮度I_MAX和面阵探测器噪声大小设置亮度阈值I_TH，亮度阈值I_TH大于面阵探测器噪声，且小于最大亮度I_MAX的20%。根据I_TH滤除背景噪声，将小于I_TH的所有像素的亮度置为零。

5.4计算光阑图像的各像素行的各直线中心点坐标。滤除背景噪声后，光阑图像的每条直线在包含该直线的图像像素行上对应一段非零的亮度数据，根据该亮度数据可以计算该直线在图像像素行上的中心点坐标。将面阵探测器的光阑图像行序号记为m，1≤m≤M，M为光阑图像的像素总行数，从下向上依次由小到大编号；将光阑图像列序号记为n，1≤n≤N，N为光阑图像的像素总列数，从左向右依次由小到大编号。计算第1至第M行光阑图像的各直线中心点坐标，具体方法为：

5.4.1令m=1；

5.4.2读取第m行光阑图像数据，根据各像素的亮度大小和位置计算所有非零区段；根据各非零区段长度判断区段的有效性，若非零区段的像素长度大于5且小于图像行总长的1/4则该非零区段有效，否则该非零区段无效；

5.4.3计算第m行光阑图像的有效非零区段总数K_m。建立面阵探测器坐标系，坐标系原点O为面阵探测器中心，X轴和Z轴分别沿面阵探测器的水平中心线和垂直中心线。将第m行、n列图像像素的灰度记为I(m,n)。计算第m行光阑图像的第1至第K_m个有效非零区段的中心坐标，具体方法为：

5.4.3.1令k_m=1；

5.4.3.2采用行亮度重心方法或行曲线拟合方法计算第k_m个有效非零区段的中心坐标；

行亮度重心方法为：按式(一)和式(二)分别计算光阑图像第m行图像第k_m个有效非零区段的中心水平坐标X_mkm和垂直坐标Z_mkm，其中n_kms为第k_m个有效非零区段的起始图像列序号，n_kme为第k_m个有效非零区段的结束图像列序号。由于k_m作为X_m的下标后k_m的下标中的“m”特别小，导致显示不清楚，故在用k_m作为下标时均写成“km”，如“X_mkm”所示，以下所有以k_m作为下标的标识都采用这种办法。

$X_{\mathrm{mkm}}=\frac{\underset{n=n_{\mathrm{kms}}}{\overset{n_{\mathrm{kme}}}{\mathrm{\Σ}}}\mathrm{nI}(m,n)}{\underset{n=n_{\mathrm{kms}}}{\overset{n_{\mathrm{kme}}}{\mathrm{\Σ}}}I(m,n)}-\frac{N}{2}$ (一)

$Z_{\mathrm{mkm}}=m-\frac{M}{2}$ (二)

行曲线拟合方法为：按式(三)对光阑图像第m行图像第k_m个有效非零区段的亮度进行P阶多项式拟合，J_kmp为P阶多项式的第p个系数（1≤p≤P）。阶数P为大于1的常数，满足使拟合曲线计算的亮度与第k_m个有效非零区段的亮度的相对误差小于5%。

$I_{\mathrm{km}}^{\′}(m,n)=\underset{p=0}{\overset{P}{\mathrm{\Σ}}}{J}_{\mathrm{kmp}}{n}^p$ (n_kms≤n≤n_kme)    (三)

采用式(四)计算多项式系数

$\left[\begin{array}{c}{J}_{\mathrm{km}1}\\ J_{\mathrm{km}2}\\ \·\\ \·\\ \·\\ J_{\mathrm{kmP}}\end{array}\right]={\left[N_{\mathrm{mkm}}^T{N}_{\mathrm{mkm}}\right]}^{-1}{N}_{\mathrm{mkm}}^T{I}_{\mathrm{mkm}}$ (四)

其中

$N_{\mathrm{mkm}}=\left[\begin{array}{cccc}{n}_{\mathrm{kms}}& n_{\mathrm{kms}}^2& \·\·\·& n_{\mathrm{kms}}^P\\ n_{\mathrm{kms}}+1& {(n_{\mathrm{kms}}+1)}^2& \·\·\·& {(n_{\mathrm{kms}}+1)}^P\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ n_{\mathrm{kme}}& n_{\mathrm{kme}}^2& \·\·\·& n_{\mathrm{kme}}^P\end{array}\right]$ $I_{\mathrm{mkm}}=\left[\begin{array}{c}I(m,n_{\mathrm{kms}})\\ I(m,n_{\mathrm{kms}}+1)\\ \·\\ \·\\ \·\\ I(m,n_{\mathrm{kme}})\end{array}\right]$

矩阵为N_mkm的转置矩阵。采用上述方法求解出第m行图像第k_m个有效非零区段的亮度分布曲线后，采用式(五)计算该有效非零区段的中心水平坐标X_mkm，按式(二)计算垂直坐标Z_mkm。

$X_{\mathrm{mkm}}=\frac{{\∫}_{n_{\mathrm{kms}}}^{n_{\mathrm{kme}}}{\mathrm{nI}}_{\mathrm{km}}^{\′}(m,n)\mathrm{dn}}{{\∫}_{n_{\mathrm{kms}}}^{n_{\mathrm{kme}}}{I}_{\mathrm{km}}^{\′}(m,n)\mathrm{dn}}-\frac{N}{2}$ (五)

5.4.3.3 k_m＝k_m+1；若k_m大于K_m，则转入步骤5.4.4，否则转入步骤5.4.3.2；注：若第k_m次使用行亮度重心方法，则第k_m+1次也使用行亮度重心方法，若第k_m次使用列曲线拟合方法，则第k_m+1次也使用列曲线拟合方法；

5.4.4 m=m+1；若m大于M，则转入步骤5.5，否则转入步骤5.4.2；

5.5将所有图像行的有效非零区段的中心坐标的总数记为G_M，对第1个至第G_M个中心坐标进行分组，总共分为L_C组，对应光阑图像的L_C条直线。具体方法为：

5.5.1建立2L_C个二维数组，用于保存各组中心坐标，每个二维数组的长度为M，将第L4个二维数组记为DM1_L4（1≤L4≤2L_C）；建立2L_C个整数，用于表示每组中心坐标的有效中心坐标数量，将第L4个整数记为EM1_L4（1≤L4≤2L_C），将所有二维数组DM1_L4和整数EM1_L4清零；

5.5.2令g_M=1；

5.5.3对第g_M个中心坐标进行分组，具体方法为：

5.5.3.1令L4=1；

5.5.3.2若EM1_L4=0，则将第g_M个中心坐标复制到DM1_L4的第1个坐标数据中，EM1_L4=1，进入步骤5.5.4；否则进入步骤5.5.3.3；

5.5.3.3分别计算第g_M个中心坐标与DM1_L4的第1个至第EM1_L4个坐标数据之间的最小距离，若最小距离大于5个像素则进入步骤5.5.3.4；否则将第g_M个中心坐标复制到DM1_L4的第EM1_L4+1个坐标数据中，EM1_L4＝EM1_L4+1，进入步骤5.5.4；

5.5.3.4 L4=L4+1，若L4大于2L_C，则进入步骤5.5.4；否则进入步骤5.5.3.2；

5.5.4 g_M=g_M+1，若g_M大于G_M，进入步骤5.5.5；否则进入步骤5.5.3；

5.5.5对2L_C组中心坐标进行合并，将光阑图像中每条直线的两段合并为一条直线，合并后的中心坐标组数为L_C组，具体方法为：

5.5.5.1建立L_C个二维数组，用于保存各组中心坐标，每个二维数组的长度为2M，将第L个二维数组记为DM2_L（1≤L≤L_C）；建立L_C个整数，用于表示每组中心坐标的有效中心坐标数量，将第L个整数记为EM2_L（1≤L≤L_C）；建立2L_C个整数，用于表示每组中心坐标的合并状态（0表示未合并，1表示已合并），将第L4个整数记为FM2_L4（1≤L4≤2L_C）；将所有二维数组DM2_L、整数EM2_L和整数FM2_L4清零；

5.5.5.2令L4=1，L=1；

5.5.5.3若FM2_L4=1，则进入步骤5.5.5.9；否则将二维数组DM1_L4的所有中心坐标数据复制到二维数组DM2_L的第1至EM1_L4个数据中，然后令EM2_L=EM1_L4，FM2_L4=1；

5.5.5.4令L5=L4+1；

5.5.5.6若FM2_L5=1，则进入步骤5.5.5.8；否则采用最小二乘法分别计算二维数组DM1_L4和DM1_L5的直线方程，然后计算两直线之间的最小距离；

5.5.5.7若最小距离大于5个像素则进入步骤5.5.5.8；否则将二维数组DM1_L5的所有中心坐标数据复制到二维数组DM2_L的第EM2_L+1至EM2_L+EM1_L5个数据中，然后令EM2_L=EM2_L+EM1_L5，FM2_L5=1，L=L+1，进入步骤5.5.5.9；

5.5.5.8 L5=L5+1；若L5大于2L_C，则进入步骤5.5.5.9；否则进入步骤5.5.5.6；

5.5.5.9 L4=L4+1；

5.5.5.10若L4大于2L_C，进入步骤5.6；否则进入步骤5.5.5.3；

5.6计算第1组至第L_C组中心坐标的直线方程，具体方法为：

5.6.1令L=1；

5.6.2将第L组中心坐标分别记为X_MLa、Z_MLa，1≤a≤A_L，A_L为第L组中心坐标总数。采用式(六)表示第L组中心坐标的直线方程，其中B_ML、C_ML为直线方程参数；按式(七)计算直线方程的各系数。

X=B_MLZ+C_ML    (六)

$\left[\begin{array}{c}{B}_{\mathrm{ML}}\\ C_{\mathrm{ML}}\end{array}\right]={\left[Z_{\mathrm{ML}}^T{Z}_{\mathrm{ML}}\right]}^{-1}{Z}_{\mathrm{ML}}^T{X}_{\mathrm{ML}}$ (七)

其中

$X_{\mathrm{ML}}=\left[\begin{array}{c}{X}_{\mathrm{ML}1}\\ X_{\mathrm{ML}2}\\ \·\\ \·\\ \·\\ X_{{\mathrm{MLA}}_L}\end{array}\right]$ $Z_{\mathrm{ML}}=\left[\begin{array}{cc}{Z}_{\mathrm{ML}1}& 1\\ Z_{\mathrm{ML}2}& 1\\ \·& \·\\ \·& \·\\ \·& \·\\ Z_{{\mathrm{MLA}}_L}& 1\end{array}\right]$

矩阵为Z_ML的转置矩阵。

5.6.3按式(八)计算第L个直线的倾斜角

(八)

5.6.4 L=L+1，若L大于L_C，则进入步骤5.7；否则进入步骤5.6.2；

5.7根据第1至第L_C个直线倾斜角计算横滚方向变形按式(九)进行计算；

(九)

5.8计算光阑图像的各像素列的各直线中心点坐标。光阑图像的每条直线在包含该直线的图像像素列上对应一段非零的亮度数据，根据该亮度数据可以计算该直线在图像像素列上的中心点坐标。计算第1至第N列光阑图像的各直线中心点坐标，具体方法为：

5.8.1令n=1；

5.8.2读取第n列光阑图像数据，根据各像素的亮度大小和位置计算所有非零区段；根据各非零区段长度判断区段的有效性，若非零区段的像素长度大于5且小于图像列总长的1/4则该非零区段有效，否则该非零区段无效；

5.8.3计算第n列光阑图像的有效非零区段总数K_n。计算第n列光阑图像的第1至第K_n个有效非零区段的中心坐标，具体方法为：

5.8.3.1令k_n=1；

5.8.3.2采用列亮度重心方法或列曲线拟合方法计算第k_n个有效非零区段的中心坐标；

列亮度重心方法为：按式(十)和式(十一)分别计算光阑图像第n列图像第k_n个有效非零区段的中心水平坐标X_nkn和垂直坐标Z_nkn，其中m_kns为第k_n个有效非零区段的起始图像行序号，m_kne为第k_n个有效非零区段的结束图像行序号。由于k_n作为X_n的下标后k_n的下标中的“n”特别小，导致显示不清楚，故在用k_n作为下标时均写成“kn”，如“X_nkn”所示，以下所有以k_n作为下标的标识都采用这种办法。

$X_{\mathrm{nkn}}=n-\frac{N}{2}$ (十)

$Z_{\mathrm{nkn}}=\frac{\underset{m=m_{\mathrm{kns}}}{\overset{m_{\mathrm{kne}}}{\mathrm{\Σ}}}\mathrm{mI}(m,n)}{\underset{m=m_{\mathrm{kns}}}{\overset{m_{\mathrm{kne}}}{\mathrm{\Σ}}}I(m,n)}-\frac{M}{2}$ (十一)

列曲线拟合方法为：按式(十二)对光阑图像第n列图像第k_n个有效非零区段的亮度进行P阶多项式拟合，其中J_knp为P阶多项式的第p个系数（1≤p≤P）。阶数P为大于1的常数，满足使拟合曲线计算的亮度与第k_n个有效非零区段的亮度的相对误差小于5%。

$I_{\mathrm{kn}}^{\′}(m,n)=\underset{p=0}{\overset{P}{\mathrm{\Σ}}}{J}_{\mathrm{knp}}{m}^p$ (m_kns≤m≤m_kne)    (十二)

采用式(十三)计算多项式系数

$\left[\begin{array}{c}{J}_{\mathrm{kn}1}\\ J_{\mathrm{kn}2}\\ \·\\ \·\\ \·\\ J_{\mathrm{knP}}\end{array}\right]={\left[N_{\mathrm{nkn}}^T{N}_{\mathrm{nkn}}\right]}^{-1}{N}_{\mathrm{nkn}}^T{I}_{\mathrm{nkn}}$ (十三)

其中

$N_{\mathrm{nkn}}=\left[\begin{array}{cccc}{m}_{\mathrm{kns}}& m_{\mathrm{kns}}^2& \·\·\·& m_{\mathrm{kns}}^P\\ m_{\mathrm{kns}}+1& {(m_{\mathrm{kns}}+1)}^2& \·\·\·& {(m_{\mathrm{kns}}+1)}^P\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ m_{\mathrm{kne}}& m_{\mathrm{kne}}^2& \·\·\·& m_{\mathrm{kne}}^P\end{array}\right]$ $I_{\mathrm{nkn}}=\left[\begin{array}{c}I(m_{\mathrm{kns}},n)\\ I(m_{\mathrm{kns}}+1,n)\\ \·\\ \·\\ \·\\ I(m_{\mathrm{kne}},n)\end{array}\right]$

矩阵为N_nkn的转置矩阵。采用上述方法求解出第n列图像第k_n个有效非零区段的亮度分布曲线后，采用式(十四)计算该有效非零区段的中心垂直坐标Z_nkn，按式(十)计算水平坐标X_nkn。

$Z_{\mathrm{nkn}}=\frac{{\∫}_{m_{\mathrm{kns}}}^{m_{\mathrm{kne}}}{\mathrm{mI}}^{\′}(m,n)\mathrm{dm}}{{\∫}_{m_{\mathrm{kns}}}^{m_{\mathrm{kne}}}{I}^{\′}(m,n)\mathrm{dm}}-\frac{M}{2}$ (十四)

5.8.3.3 k_n=k_n+1；若k_n大于K_n，则转入步骤5.8.4，否则转入步骤5.8.3.2；注：若第k_n次使用列亮度重心方法，则第k_n+1次也使用列亮度重心方法，若第k_n次使用列曲线拟合方法，则第k_n+1次也使用列曲线拟合方法；

5.8.4 n=n+1；若n大于N，则转入步骤5.9，否则转入步骤5.8.2；

5.9将所有图像列的有效非零区段的中心坐标的总数记为G_N，对第1个至第G_N个中心坐标进行分组，总共分为L_C组，对应光阑图像的L_C条直线。具体方法为：

5.9.1建立2L_C个二维数组，用于保存各组中心坐标，每个二维数组的长度为N，将第L4个二维数组记为DN1_L4（1≤L4≤2L_C）；建立2L_C个整数，用于表示每组中心坐标的有效中心坐标数量，将第L4个整数记为EN1_L4（1≤L4≤2L_C），将所有二维数组DN1_L4和整数EN1_L4清零；

5.9.2令g_N=1；

5.9.3对第g_N个中心坐标进行分组，具体方法为：

5.9.3.1令L4=1；

5.9.3.2若EN1_L4=0，则将第g_N个中心坐标复制到DN1_L4的第1个坐标数据中，EN1_L4=1，进入步骤5.9.4；否则进入步骤5.9.3.3；

5.9.3.3分别计算第g_N个中心坐标与DN1_L4的第1个至第EN1_L4个坐标数据之间的最小距离，若最小距离大于5个像素则进入步骤5.9.3.4；否则将第g_N个中心坐标复制到DN1_L4的第EN1_L4+1个坐标数据中，EN1_L4＝EN1_L4+1，进入步骤5.9.4；

5.9.3.4 L4=L4+1，若L4大于2L_C，则进入步骤5.9.4；否则进入步骤5.9.3.2；

5.9.4 g_N=g_N+1，若gN大于G_N，进入步骤5.9.5；否则进入步骤5.9.3；

5.9.5对2L_C组中心坐标进行合并，将光阑图像中每条直线的两段合并为一条直线，合并后的中心坐标组数为L_C组，具体方法为：

5.9.5.1建立L_C个二维数组，用于保存各组中心坐标，每个二维数组的长度为2N，将第L个二维数组记为DN2_L（1≤L≤L_C）；建立L_C个整数，用于表示每组中心坐标的有效中心坐标数量，将第L个整数记为EN2_L（1≤L≤L_C）；建立2L_C个整数，用于表示每组中心坐标的合并状态（0表示未合并，1表示已合并），将第L4个整数记为FN2_L4（1≤L4≤2L_C）；将所有二维数组DN2_L、整数EN2_L和整数FN2_L4清零；

5.9.5.2令L4=1，L=1；

5.9.5.3若FN2_L4=1，则进入步骤5.9.5.9；否则将二维数组DN1_L4的所有中心坐标数据复制到二维数组DN2_L的第1至EN1_L4个数据中，然后令EN2_L=EM1_L4，FN2_L4=1；

5.9.5.4令L5=L4+1；

5.9.5.6若FN2_L5=1，则进入步骤5.9.5.8；否则采用最小二乘法分别计算二维数组DN1_L4和DN1_L5的直线方程，然后计算两直线之间的最小距离；

5.9.5.7若最小距离大于5个像素则进入步骤5.9.5.8；否则将二维数组DN1_L5的所有中心坐标数据复制到二维数组DN_2L的第EN2L+1至EN2_L+EN1_L5个数据中，然后令EN2_L=EN2_L+EN1_L5，FN2_L5=1，L=L+1，进入步骤5.9.5.9；

5.9.5.8 L5=L5+1；若L5大于2L_C，则进入步骤5.9.5.9；否则进入步骤5.9.5.6；

5.9.5.9 L4=L4+1；

5.9.5.10若L4大于2L_C，进入步骤5.10；否则进入步骤5.9.5.3；

5.10计算第1组至第L_C组中心坐标的直线方程，具体方法为：

5.10.1令L=1；

5.10.2将第L组中心坐标分别记为X_NLb、Z_NLb，1≤b≤B_L，B_L为第L组中心坐标总数。采用式(十五)表示第L组中心坐标的直线方程，其中B_NL、C_NL为直线方程参数；按式(十六)计算直线方程的各系数。

X=B_NLZ+C_NL    (十五)

$\left[\begin{array}{c}{B}_{\mathrm{NL}}\\ C_{\mathrm{NL}}\end{array}\right]={\left[Z_{\mathrm{NL}}^T{Z}_{\mathrm{NL}}\right]}^{-1}{Z}_{\mathrm{NL}}^T{X}_{\mathrm{NL}}$ (十六)

其中

$X_{\mathrm{NL}}=\left[\begin{array}{c}{X}_{\mathrm{NL}1}\\ X_{\mathrm{NL}2}\\ \·\\ \·\\ \·\\ X_{{\mathrm{NLB}}_L}\end{array}\right]$ $Z_{\mathrm{NL}}=\left[\begin{array}{cc}{Z}_{\mathrm{NL}1}& 1\\ Z_{\mathrm{NL}2}& 1\\ \·& \·\\ \·& \·\\ \·& \·\\ Z_{{\mathrm{NLB}}_L}& 1\end{array}\right]$

矩阵为Z_NL的转置矩阵。

5.10.3按式(十七)计算第L个直线的倾斜角

十七

5.10.4 L=L+1，若L大于L_C，则进入步骤5.11；否则进入步骤5.10.2；

5.11根据第1至第L_C个直线倾斜角计算横滚方向变形按式(十八)进行计算；(十八)

(十八)

5.12计算最终的横滚方向变形按式(十九)其中为横滚角零位。

(十九)

5.13根据光阑图像数据计算方位和俯仰方向变形。具体计算方法为：首先计算整个光阑图像的中心坐标，采用式(二十)计算整个光阑图像中心的水平坐标采用式(二十一)计算整个光阑图像中心的垂直坐标然后采用式(二十二)计算方位方向变形ψ，采用式(二十三)计算俯仰方向变形θ，其中f为接收光学系统焦距，μ_X为面阵探测器水平方向的像元间距，μ_Z为面阵探测器垂直方向的像元间距。

$\stackrel{\‾}{X}=\frac{\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\mathrm{nI}(m,n)}{\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}I(m,n)}-\frac{N}{2}$ (二十)

$\stackrel{\‾}{Z}=\frac{\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\mathrm{mI}(m,n)}{\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}I(m,n)}-\frac{M}{2}$ (二十一)

$\mathrm{\ψ}=\frac{\stackrel{\‾}{X}{\mathrm{\μ}}_X}{f}$ (二十二)

$\mathrm{\θ}=\frac{\stackrel{\‾}{Z}{\mathrm{\μ}}_Z}{f}$ (二十三)

本发明可以达到以下技术效果：

(1)本发明装置采用发射光路和接收光路分离的准直光路结构，可以有效克服发射模块和接收模块之间的位置变化引起的变形测量误差。

(2)本发明装置采用一套发射准直光路、一套接收准直光路、一个面阵探测器即可实现三维变形测量，大幅度降低了整个测量装置的复杂性、成本和占用体积，明显提高了系统的可靠性和稳定性；

(3)采用两条或多条直线结构光阑和相应的变形解算方法，能同时测量方位角、俯仰角和横滚角变形，能实现角秒级三维变形测量，具有测量精度高、结构简单和易于实施等优点。

附图说明

图1是背景技术中国专利CN97251378.7、美国专利US7227627、《中国惯性技术学报》杂志2006年第3期中“一种舰船变形测量方法的研究与应用”和《长春理工大学学报》杂志2006年第3期中“车载平台变形的激光自准直测量方法研究”论文公布的二维变形测量装置示意图。

图2是背景技术美国专利US7697127公布的旋转角测量装置示意图。

图3是背景技术中国专利CN200410010856.3、《计量技术》杂志2005年第4期中“干涉法在船体角度变形测量中的应用”、《光学技术》杂志2005年第5期中“船体三维角度变形的自准直干涉测量方法”和《飞行器测控学报》2009年第3期中“基于莫尔条纹测量法设计的船体变形测量系统”论文公布的干涉型三维变形测量装置示意图。

图4是背景技术《光学精密工程》杂志2008年第11期中“船体变形测量的新技术”和《长春理工大学学报》杂志2006年第3期中“基于测量船角变形光电测量系统的研究”论文公布的大钢管型三维变形测量装置示意图。

图5是本发明三维变形测量装置结构示意图。

图6是本发明采用的光阑的直线结构示意图。

图7是无变形时面阵探测器探测的光阑图像。

图8是有方位和俯仰方向变形时面阵探测器探测的光阑图像。

图9是有横滚方向变形时面阵探测器探测的光阑图像。

图10为本发明变形测量方法整体流程图。

图11是光阑图像的某行图像数据示意图。

具体实施方式

本发明三维变形测量装置由发射模块、接收模块和处理模块三大部分构成，主要结构如图5所示。发射模块33与接收模块34之间能够通视。该装置能实时测量参考基准23与待测物22之间的三维变形。

发射模块33由光源1、光阑2和发射光学系统31构成。采用光学固定座和调节架等将光源1、光阑2和发射光学系统31固定在参考基准23。光阑2位于光源1和发射光学系统31之间，且位于发射光学系统31的焦平面上，光源1和光阑2的中心都在发射光学系统31的光轴上，由此构成发射准直光路。光源1照射光阑2，发射光学系统31将透过光阑2的光束准直发射到接收模块34中。光源1可以是半导体激光器或高亮LED或高亮灯泡等。发射光学系统31可以采用双胶合透镜或双分离透镜或多片透镜组等类型光学系统。光阑2是在玻璃等透明材料上通过光刻等工艺加工而成的，加工后的光阑2有些部分透光，有些部分不透光。光阑2的透光部分采用多条直线结构。

接收模块34由接收光学系统32和面阵探测器8构成。采用光学固定座和调节架等将接收光学系统32和面阵探测器8固定在待测物22上。面阵探测器8位于接收光学系统4的焦平面上，且其中心位于接收光学系统32光轴上，由此构成接收准直光路。接收光学系统32可采用双胶合透镜或双分离透镜或多片透镜组等类型光学系统，面阵探测器8可采用灰度或彩色的面阵CCD或CMOS等类型探测器。发射模块33发射的准直光束经过接收光学系统32会聚后成像在面阵探测器8上，面阵探测器8探测的是光阑2的直线结构光阑图像。

处理模块35可采用工控机、嵌入式计算机、台式计算机或笔记本电脑等设备，在处理模块35中包含变形解算软件。处理模块35通过千兆以太网、高速USB或1394等接口与接收模块34中的面阵探测器8连接，实时采集面阵探测器8探测的光阑图像，变形解算软件计算三维变形。

图6是本发明采用的光阑的直线结构示意图。光阑2是在玻璃等透明材料上利用光刻等工艺加工而成的，加工后的光阑2有些部分透光，有些部分不透光。光阑2的透光部分采用L_C（L_C为大于1的常数）条直线结构，各直线呈放射状分布，直线结构为透光部分，其余部分不透光。两直线结构光阑（如图6(a)所示）的透光部分由两条直线构成，两直线可互相垂直或不垂直。每条直线由两个线段构成，其中线段L1a、L1b构成直线L1，线段L2a、L2b构成直线L2。采用多直线结构光阑可以提高变形测量精度，三条直线结构光阑如6(b)所示，四条直线结构光阑如6(c)所示，四条以上的直线结构光阑依次类推。

图7是无变形时面阵探测器探测的光阑图像。图中以两直线结构的光阑为例，光阑图像中包含两条直线，与光阑2的直线结构对应，且光阑图像中心与面阵探测器8的中心重合。

图8是有方位和俯仰方向变形时面阵探测器探测的光阑图像。图中以两直线结构的光阑为例，当发射模块33与接收模块34之间存在方位和俯仰方向变形时，面阵探测器8探测的光阑图像发生平移，光阑图像中心偏离面阵探测器8中心。将光阑图像中心与面阵探测器8中心在X轴方向的偏移量记为ΔX，Z轴方向的偏移量记为ΔZ。

图9是有横滚方向变形时面阵探测器探测的光阑图像。图中以两直线结构的光阑为例，当发射模块33与接收模块34之间存在横滚方向变形时，面阵探测器8的光阑图像发生旋转。光阑图像中直线方向与面阵探测器8中心线存在夹角该角度为横滚方向变形。

图10为本发明变形测量方法整体流程图，主要包含五个步骤：

第一步，发射模块33的光源1照射光阑2，发射光学系统31将透过光阑2的光束准直发射到接收模块。

第二步，接收模块34的接收光学系统32接收发射模块发射的光束，并将其会聚到面阵探测器8上。

第三步，面阵探测器8将接收的光束转换成光阑图像，此光阑图像为数字化的彩色图像或灰度图像。

第四步，处理模块35从面阵探测器8接收光阑图像。

第五步，处理模块35中变形解算软件根据光阑图像计算三维变形。

图11是光阑图像的第m行图像数据示意图。滤除背景噪声后，光阑图像的每条直线在包含该直线的像素行上对应一段非零的亮度数据。该图像像素行中总共包含K_m段非零亮度数据，编号依次从1至K_m，k_m为某段非零亮度数据的编号（1≤k_m≤K_m）。

Claims (5)

1.一种基于准直光路的角秒级三维变形测量装置，其特征在于三维变形测量装置由发射模块(33)、接收模块(34)和处理模块(35)三大部分构成，发射模块(33)与参考基准(23)固连，接收模块(34)与待测物(22)固连，发射模块(33)与接收模块(34)之间能够通视，发射模块(33)发射的光束直接进入接收模块(34)；

发射模块(33)由光源(1)、光阑(2)和发射光学系统(31)构成，光源(1)、光阑(2)和发射光学系统(31)固定在参考基准(23)上；光阑(2)位于光源(1)和发射光学系统(31)之间，且位于发射光学系统(31)的焦平面上，光源(1)和光阑(2)的中心都在发射光学系统(31)的光轴上，由此构成发射准直光路；光源(1)照射光阑(2)，发射光学系统(31)将透过光阑(2)的光束准直发射到接收模块(34)中；光源(1)是半导体激光器或LED或灯泡；发射光学系统(31)采用双胶合透镜或双分离透镜或多片透镜组；光阑(2)是在玻璃等透明材料上利用光刻工艺加工而成的，加工后的光阑(2)有些部分透光，有些部分不透光，光阑(2)的透光部分采用L_C条直线结构，各直线呈放射状分布，直线结构为透光部分，其余部分不透光，L_C为大于1的常数；每条直线由两个线段构成，直线越多，测量精度越高，但横滚方向的形变测量范围会减小；

接收模块(34)由接收光学系统(32)和面阵探测器(8)构成，接收光学系统(32)和面阵探测器(8)固定在待测物(22)上；面阵探测器(8)位于接收光学系统(32)的焦平面上，且其中心位于接收光学系统(32)的光轴上，由此构成接收准直光路；接收光学系统(32)采用双胶合透镜或双分离透镜或多片透镜组，面阵探测器(8)采用灰度或彩色的面阵CCD或CMOS探测器；发射模块(33)发射的准直光束经过接收光学系统(32)会聚后成像在面阵探测器(8)上，面阵探测器(8)探测光阑(2)的直线结构光阑图像；

处理模块(35)采用工控机、嵌入式计算机、台式计算机或笔记本电脑，在处理模块(35)中安装了变形解算软件，处理模块(35)通过千兆以太网、高速USB或1394接口与面阵探测器(8)连接，实时采集面阵探测器(8)探测的光阑图像，变形解算软件计算参考基准(23)与待测物(22)之间的三维变形，包括方位角、俯仰角和横滚角变形。

2.一种采用如权利要求1所述基于准直光路的角秒级三维变形测量装置进行三维变形测量的方法，其特征在于包括以下步骤：

第一步，发射模块(33)的光源(1)照射光阑(2)，发射光学系统(31)将透过光阑(2)的光束准直发射到接收模块(34)；

第二步，接收模块(34)的接收光学系统(32)接收发射模块发射的光束，并将其会聚到面阵探测器(8)上；

第三步，面阵探测器(8)将接收的光束转换成光阑图像，此光阑图像为数字化的彩色图像或灰度图像；

第四步，处理模块(35)从面阵探测器(8)接收光阑图像；

第五步，处理模块(35)中变形解算软件根据光阑图像按以下步骤计算三维变形：

5.1 若面阵探测器(8)是彩色类型探测器，则将彩色图像转换成灰度图像；若面阵探测器(8)是灰度类型探测器，则直接进入步骤5.2；

5.2 采用传统标定或自标定方法对发射光学系统畸变、接收光学系统畸变和面阵探测器畸变引起的光阑图像畸变进行校正；

5.3 对校正后的光阑图像进行亮度分析，并滤除背景噪声；

5.4 将面阵探测器(8)的光阑图像行序号记为m，1≤m≤M，M为光阑图像的像素总行数，从下向上依次由小到大编号；将光阑图像列序号记为n，1≤n≤N，N为光阑图像的像素总列数，从左向右依次由小到大编号，计算第1至第M行光阑图像的各直线中心点坐标，具体方法为：

5.4.1 令m＝1；

5.4.2 读取第m行光阑图像数据，根据各像素的亮度大小和位置计算所有非零区段；根据各非零区段长度判断区段的有效性，若非零区段的像素长度大于5且小于图像行总长的1/4则该非零区段有效，否则该非零区段无效；

5.4.3 统计第m行光阑图像的有效非零区段总数K_m；建立面阵探测器坐标系，坐标系原点O为面阵探测器(8)中心，X轴和Z轴分别沿面阵探测器的水平中心线和垂直中心线，将第m行、n列图像像素的灰度记为I(m,n)，计算第m行光阑图像的第1至第K_m个有效非零区段的中心坐标，具体方法为：

5.4.3.1 令k_m＝1；

5.4.3.2 采用行亮度重心方法计算第k_m个有效非零区段的中心坐标，行亮度重心方法为：按式(一)和式(二)分别计算光阑图像第m行图像第k_m个有效非零区段的中心水平坐标X_mkm和垂直坐标Z_mkm，其中n_kms为第k_m个有效非零区段的起始图像列序号，n_kme为第k_m个有效非零区段的结束图像列序号：

$X_{\mathrm{mkm}}=\frac{\underset{n=n_{\mathrm{kms}}}{\overset{n_{\mathrm{kme}}}{\mathrm{\Σ}}}\mathrm{nI}(m,n)}{\underset{n=n_{\mathrm{kms}}}{\overset{n_{\mathrm{kme}}}{\mathrm{\Σ}}}I(m,n)}-\frac{N}{2}$    (一)

$Z_{\mathrm{mkm}}=m-\frac{M}{2}$    (二)；

5.4.3.3 k_m＝k_m+1；若k_m大于K_m，则转入步骤5.4.4，否则转入步骤5.4.3.2，若第k_m次使用行亮度重心方法，则第k_m+1次也使用行亮度重心方法；

5.4.4 m＝m+1；若m大于M，转入步骤5.5，否则转入步骤5.4.2；

5.5 将所有图像行的有效非零区段的中心坐标的总数记为G_M，对第1个至第G_M个中心坐标进行分组，总共分为L_C组，对应光阑图像的L_C条直线，具体方法为：

5.5.1 建立2L_C个二维数组，用于保存各组中心坐标，每个二维数组的长度为M，将第L4个二维数组记为DM1_L4，1≤L4≤2L_C；建立2L_C个整数，用于表示每组中心坐标的有效中心坐标数量，将第L4个整数记为EM1_L4，将所有二维数组DM1_L4和整数EM1_L4清零；

5.5.2 令g_M＝1；

5.5.3 对第g_M个中心坐标进行分组，具体方法为：

5.5.3.1 令L4＝1；

5.5.3.2 若EM1_L4＝0，则将第g_M个中心坐标复制到DM1_L4的第1个坐标数据中，EM1_L4＝1，进入步骤5.5.4；否则进入步骤5.5.3.3；

5.5.3.3 分别计算第g_M个中心坐标与DM1_L4的第1个至第EM1_L4个坐标数据之间的最小距离，若最小距离大于5个像素则进入步骤5.5.3.4；否则将第g_M个中心坐标复制到DM1_L4的第EM1_L4+1个坐标数据中，EM1_L4＝EM1_L4+1，进入步骤5.5.4；

5.5.3.4L4＝L4+1，若L4大于2L_C，则进入步骤5.5.4；否则进入步骤5.5.3.2；

5.5.4g_M＝g_M+1，若g_M大于G_M，进入步骤5.5.5；否则进入步骤5.5.3；

5.5.5 对2L_C组中心坐标进行合并，将光阑图像中每条直线的两段合并为一条直线，合并后的中心坐标组数为L_C组，具体方法为：

5.5.5.1 建立L_C个二维数组，用于保存各组中心坐标，每个二维数组的长度为2M，将第L个二维数组记为DM2_L，1≤L≤L_C；建立L_C个整数，用于表示每组中心坐标的有效中心坐标数量，将第L个整数记为EM2_L；建立2L_C个整数，用于表示每组中心坐标的合并状态，0表示未合并，1表示已合并，将第L4个整数记为FM2_L4，1≤L4≤2L_C；将所有二维数组DM2_L、整数EM2_L和整数FM2_L4清零；

5.5.5.2 令L4＝1，L＝1；

5.5.5.3 若FM2_L4＝1，则进入步骤5.5.5.9；否则将二维数组DM1_L4的所有中心坐标数据复制到二维数组DM2_L的第1至EM1_L4个数据中，然后令EM2_L＝EM1_L4，FM2_L4＝1；

5.5.5.4 令L5＝L4+1；

5.5.5.6 若FM2_L5＝1，则进入步骤5.5.5.8；否则采用最小二乘法分别计算二维数组DM1_L4和DM1_L5的直线方程，然后计算两直线之间的最小距离；

5.5.5.7 若最小距离大于5个像素则进入步骤5.5.5.8；否则将二维数组DM1_L5的所有中心坐标数据复制到二维数组DM2_L的第EM2_L+1至EM2_L+EM1_L5个数据中，然后令EM2_L＝EM2_L+EM1_L5，FM2_L5＝1，L＝L+1，进入步骤5.5.5.9；

5.5.5.8 L5＝L5+1；若L5大于2L_C，则进入步骤5.5.5.9；否则进入步骤5.5.5.6；

5.5.5.9 L4＝L4+1；

5.5.5.10 若L4大于2L_C，进入步骤5.6；否则进入步骤5.5.5.3；

5.6 计算第1组至第L_C组中心坐标的直线方程，具体方法为：

5.6.1 令L＝1；

5.6.2 将第L组中心坐标分别记为X_MLa、Z_MLa，1≤a≤A_L，A_L为第L组中心坐标总数，采用式(六)表示第L组中心坐标的直线方程，其中B_ML、C_ML为直线方程参数；按式(七)计算直线方程的各系数：

X＝B_MLZ+C_ML   (六)

$\left[\begin{array}{c}{B}_{\mathrm{ML}}\\ C_{\mathrm{ML}}\end{array}\right]={\left[\begin{array}{cc}{Z}_{\mathrm{ML}}^T& Z_{\mathrm{ML}}\end{array}\right]}^{-1}{Z}_{\mathrm{ML}}^T{X}_{\mathrm{ML}}$    (七)

其中

$X_{\mathrm{ML}}=\left[\begin{array}{c}{X}_{\mathrm{ML}1}\\ X_{\mathrm{ML}2}\\ .\\ .\\ .\\ X_{\mathrm{ML}{A}_L}\end{array}\right]$ $Z_{\mathrm{ML}}=\left[\begin{array}{cc}{Z}_{\mathrm{ML}1}& 1\\ Z_{\mathrm{ML}2}& 1\\ .& .\\ .& .\\ .& .\\ Z_{\mathrm{ML}{A}_L}& 1\end{array}\right]$

矩阵为Z_ML的转置矩阵；

5.6.3 按式(八)计算第L个直线的倾斜角

   (八)

5.6.4 L＝L+1，若L大于L_C，则进入步骤5.7；否则进入步骤5.6.2；

5.7 根据第1至第L_C个直线倾斜角计算横滚方向变形按式(九)进行计算；

   (九)

5.8 计算第1至第N列光阑图像的各直线中心点坐标，具体方法为：

5.8.1 令n＝1；

5.8.2 读取第n列光阑图像数据，根据各像素的亮度大小和位置计算所有非零区段；根据各非零区段长度判断区段的有效性，若非零区段的像素长度大于5且小于图像列总长的1/4则该非零区段有效，否则该非零区段无效；

5.8.3 统计第n列光阑图像的有效非零区段总数K_n；计算第n列光阑图像的第1至第K_n个有效非零区段的中心坐标，具体方法为：

5.8.3.1 令k_n＝1；

5.8.3.2 采用列亮度重心方法计算第k_n个有效非零区段的中心坐标，列亮度重心方法为：按式(十)和式(十一)分别计算光阑图像第n列图像第k_n个有效非零区段的中心水平坐标X_nkn和垂直坐标Z_nkn，其中m_kns为第k_n个有效非零区段的起始图像行序号，m_kne为第k_n个有效非零区段的结束图像行序号：

$X_{\mathrm{nkn}}=n-\frac{N}{2}$    (十)

$Z_{\mathrm{nkn}}=\frac{\underset{m=m_{\mathrm{kns}}}{\overset{m_{\mathrm{kne}}}{\mathrm{\Σ}}}\mathrm{mI}(m,n)}{\underset{m=m_{\mathrm{kns}}}{\overset{m_{\mathrm{kne}}}{\mathrm{\Σ}}}I(m,n)}-\frac{M}{2}$    (十一)

5.8.3.3 k_n＝k_n+1；若k_n大于K_n，则转入步骤5.8.4，否则转入步骤5.8.3.2，若第k_n次使用列亮度重心方法，则第k_n+1次也使用列亮度重心方法；

5.8.4 n＝n+1；若n大于N，则转入步骤5.9，否则转入步骤5.8.2；

5.9 将所有图像列的有效非零区段的中心坐标的总数记为G_N，对第1个至第G_N个中心坐标进行分组，总共分为L_C组，对应光阑图像的L_C条直线，具体方法为：

5.9.1 建立2L_C个二维数组，用于保存各组中心坐标，每个二维数组的长度为N，将第L4个二维数组记为DN1_L4；建立2L_C个整数，用于表示每组中心坐标的有效中心坐标数量，将第L4个整数记为EN1_L4，将所有二维数组DN1_L4和整数EN1_L4清零；

5.9.2 令g_N＝1；

5.9.3 对第g_N个中心坐标进行分组，具体方法为：

5.9.3.1 令L4＝1；

5.9.3.2 若EN1_L4＝0，则将第g_N个中心坐标复制到DN1_L4的第1个坐标数据中，EN1_L4＝1，进入步骤5.9.4；否则进入步骤5.9.3.3；

5.9.3.3 分别计算第g_N个中心坐标与DN1_L4的第1个至第EN1_L4个坐标数据之间的最小距离，若最小距离大于5个像素则进入步骤5.9.3.4；否则将第g_N个中心坐标复制到DN1_L4的第EN1_L4+1个坐标数据中，EN1_L4＝EN1_L4+1，进入步骤5.9.4；

5.9.3.4 L4＝L4+1，若L4大于2L_C，则进入步骤5.9.4；否则进入步骤5.9.3.2；

5.9.4 g_N＝g_N+1，若g_N大于G_N，进入步骤5.9.5；否则进入步骤5.9.3；

5.9.5 对2L_C组中心坐标进行合并，将光阑图像中每条直线的两段合并为一条直线，合并后的中心坐标组数为L_C组，具体方法为：

5.9.5.1 建立L_C个二维数组，用于保存各组中心坐标，每个二维数组的长度为2N，将第L个二维数组记为DN2_L；建立L_C个整数，用于表示每组中心坐标的有效中心坐标数量，将第L个整数记为EN2_L；建立2L_C个整数，用于表示每组中心坐标的合并状态，0表示未合并，1表示已合并，将第L4个整数记为FN2_L4；将所有二维数组DN2_L、整数EN2_L和整数FN2_L4清零；

5.9.5.2 令L4＝1，L＝1；

5.9.5.3 若FN2_L4＝1，则进入步骤5.9.5.9；否则将二维数组DN1_L4的所有中心坐标数据复制到二维数组DN2_L的第1至EN1_L4个数据中，然后令EN2_L＝EM1_L4，FN2_L4＝1；

5.9.5.4 令L5＝L4+1；

5.9.5.6 若FN2_L5＝1，则进入步骤5.9.5.8；否则采用最小二乘法分别计算二维数组DN1_L4和DN1_L5的直线方程，然后计算两直线之间的最小距离；

5.9.5.7 若最小距离大于5个像素则进入步骤5.9.5.8；否则将二维数组DN1_L5的所有中心坐标数据复制到二维数组DN2_L的第EN2_L+1至EN2_L+EN1_L5个数据中，然后令EN2_L＝EN2_L+EN1_L5，FN2_L5＝1，L＝L+1，进入步骤5.9.5.9；

5.9.5.8 L5＝L5+1；若L5大于2L_C，则进入步骤5.9.5.9；否则进入步骤5.9.5.6；

5.9.5.9 L4＝L4+1；

5.9.5.10 若L4大于2L_C，进入步骤5.10；否则进入步骤5.9.5.3；

5.10 计算第1组至第L_C组中心坐标的直线方程，具体方法为：

5.10.1 令L＝1；

5.10.2 将第L组中心坐标分别记为X_NLb、Z_NLb，1≤b≤B_L，B_L为第L组中心坐标总数，采用式(十五)表示第L组中心坐标的直线方程，其中B_NL、C_NL为直线方程参数；按式(十六)计算直线方程的各系数；

X＝B_NLZ+C_NL   (十五)

$\left[\begin{array}{c}{B}_{\mathrm{NL}}\\ C_{\mathrm{NL}}\end{array}\right]={\left[\begin{array}{cc}{Z}_{\mathrm{NL}}^T& Z_{\mathrm{NL}}\end{array}\right]}^{-1}{Z}_{\mathrm{NL}}^T{X}_{\mathrm{NL}}$    (十六)

其中

$X_{\mathrm{NL}}=\left[\begin{array}{c}{X}_{\mathrm{NL}1}\\ X_{\mathrm{NL}2}\\ .\\ .\\ .\\ X_{\mathrm{NL}{B}_L}\end{array}\right]$ $Z_{\mathrm{NL}}=\left[\begin{array}{cc}{Z}_{\mathrm{NL}1}& 1\\ Z_{\mathrm{NL}2}& 1\\ .& .\\ .& .\\ .& .\\ Z_{\mathrm{NL}{B}_L}& 1\end{array}\right]$

矩阵为Z_NL的转置矩阵；

5.10.3 按式(十七)计算第L个直线的倾斜角

   (十七)

5.10.4 L＝L+1，若L大于L_C，则进入步骤5.11；否则进入步骤5.10.2；

5.11 根据第1至第L_C个直线倾斜角计算横滚方向变形按式(十八)进行计算；

   (十八)

5.12 计算最终的横滚方向变形按式(十九)其中为横滚角零位：

   (十九)

5.13 根据光阑图像数据计算方位和俯仰方向变形，具体计算方法为：首先计算整个光阑图像的中心坐标，采用式(二十)计算整个光阑图像中心的水平坐标采用式(二十一)计算整个光阑图像中心的垂直坐标然后采用式(二十二)计算方位方向变形ψ，采用式(二十三)计算俯仰方向变形θ，其中f为接收光学系统焦距，μ_X为面阵探测器水平方向的像元间距，μ_Z为面阵探测器垂直方向的像元间距：

$\stackrel{\‾}{X}=\frac{\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\mathrm{nI}(m,n)}{\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}I(m,n)}-\frac{N}{2}$    (二十)

$\stackrel{\‾}{Z}=\frac{\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\mathrm{mI}(m,n)}{\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}I(m,n)}-\frac{M}{2}$    (二十一)

$\mathrm{\ψ}=\frac{\stackrel{\‾}{X}{\mathrm{\μ}}_X}{f}$    (二十二)

$\mathrm{\θ}=\frac{\stackrel{\‾}{Z}{\mathrm{\μ}}_Z}{f}$    (二十三)。

3.如权利要求2所述的采用如权利要求1所述基于准直光路的角秒级三维变形测量装置进行三维变形测量的方法，其特征在于5.3步对校正后的光阑图像进行亮度分析，并滤除背景噪声的方法是：对光阑图像的所有像素根据亮度进行排序，根据光阑(2)的直线结构的直线数量L_C选取L3个最亮像素，根据选择的最亮像素计算平均亮度，将该平均亮度作为光阑图像的最大亮度I_MAX；根据I_MAX和面阵探测器噪声大小设置亮度阈值I_TH，根据I_TH滤除背景噪声，将小于I_TH的所有像素的亮度置为零，I_TH大于面阵探测器噪声且小于I_MAX的20％，L3为L_C的10～100倍。

4.如权利要求2所述的采用如权利要求1所述基于准直光路的角秒级三维变形测量装置进行三维变形测量的方法，其特征在于5.4.3.2步采用行曲线拟合方法计算第k_m个有效非零区段的中心坐标，若第k_m次使用行曲线拟合方法，则5.4.3.3步第k_m+1次也使用行曲线拟合方法；行曲线拟合方法为：按式(三)对光阑图像第m行图像第k_m个有效非零区段的亮度进行P阶多项式拟合，其中J_kmp为P阶多项式的第p个系数，1≤p≤P，阶数P满足使拟合曲线计算的亮度与第k_m个有效非零区段的亮度的相对误差小于5％；

$I_{\mathrm{km}}^{\′}(m,n)=\underset{p=0}{\overset{P}{\mathrm{\Σ}}}{J}_{\mathrm{kmp}}{n}^p,n_{\mathrm{kms}}\≤n\≤n_{\mathrm{kme}}$    (三)

采用式(四)计算多项式系数

$\left[\begin{array}{c}{J}_{\mathrm{km}1}\\ J_{\mathrm{km}2}\\ .\\ .\\ .\\ J_{\mathrm{kmP}}\end{array}\right]={\left[N_{\mathrm{mkm}}^T{N}_{\mathrm{mkm}}\right]}^{-1}{N}_{\mathrm{mkm}}^T{I}_{\mathrm{mkm}}$    (四)

其中

$N_{\mathrm{mkm}}=\left[\begin{array}{cccc}{n}_{\mathrm{kms}}& n_{\mathrm{kms}}^2& ...& n_{\mathrm{kms}}^P\\ n_{\mathrm{kms}}+1& {(n_{\mathrm{kms}}+1)}^2& ...& {(n_{\mathrm{kms}}+1)}^P\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ n_{\mathrm{kme}}& n_{\mathrm{kme}}^2& ...& n_{\mathrm{kme}}^P\end{array}\right]$ $I_{\mathrm{mkm}}=\left[\begin{array}{c}I(m,n_{\mathrm{kms}})\\ I(m,n_{\mathrm{kms}}+1)\\ .\\ .\\ .\\ I(m,n_{\mathrm{kme}})\end{array}\right]$

矩阵为N_mkm的转置矩阵；

采用式(五)计算该有效非零区段的中心水平坐标X_mkm，按式(二)计算垂直坐标Z_mkm：

$X_{\mathrm{mkm}}=\frac{{\∫}_{n_{\mathrm{kms}}}^{n_{\mathrm{kme}}}n{I}_{\mathrm{km}}^{\′}(m,n)\mathrm{dn}}{{\∫}_{n_{\mathrm{kms}}}^{n_{\mathrm{kme}}}{I}_{\mathrm{km}}^{\′}(m,n)\mathrm{dn}}-\frac{N}{2}$    (五)。

5.如权利要求2所述的采用如权利要求1所述基于准直光路的角秒级三维变形测量装置进行三维变形测量的方法，其特征在于5.8.3.2步采用列曲线拟合方法计算第k_n个有效非零区段的中心坐标，若第k_n次使用列曲线拟合方法，则5.8.3.3步第k_n+1次也使用列曲线拟合方法；列曲线拟合方法为：按式(十二)对光阑图像第n列图像第k_n个有效非零区段的亮度进行P阶多项式拟合，其中J_knp为P阶多项式的第p个系数，1≤p≤P，阶数P满足使拟合曲线计算的亮度与第k_n个有效非零区段的亮度的相对误差小于5％；

$I_{\mathrm{kn}}^{\′}(m,n)=\underset{p=0}{\overset{P}{\mathrm{\Σ}}}{J}_{\mathrm{knp}}{m}^p,m_{\mathrm{kns}}\≤m\≤m_{\mathrm{kne}}$    (十二)

采用式(十三)计算多项式系数

$\left[\begin{array}{c}{J}_{\mathrm{kn}1}\\ J_{\mathrm{kn}2}\\ .\\ .\\ .\\ J_{\mathrm{knP}}\end{array}\right]={\left[N_{\mathrm{nkn}}^T{N}_{\mathrm{nkn}}\right]}^{-1}{N}_{\mathrm{nkn}}^T{I}_{\mathrm{nkn}}$    (十三)

其中

$N_{\mathrm{nkn}}=\left[\begin{array}{cccc}{m}_{\mathrm{kns}}& m_{\mathrm{kns}}^2& ...& m_{\mathrm{kns}}^P\\ m_{\mathrm{kns}}+1& {(m_{\mathrm{kns}}+1)}^2& ...& {(m_{\mathrm{kns}}+1)}^P\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ m_{\mathrm{kne}}& m_{\mathrm{kne}}^2& ...& m_{\mathrm{kne}}^P\end{array}\right]$ $I_{\mathrm{nkn}}=\left[\begin{array}{c}I(m_{\mathrm{kns}},n)\\ I(m_{\mathrm{kns}}+1,n)\\ .\\ .\\ .\\ I(m_{\mathrm{kne}},n)\end{array}\right]$

矩阵为N_nkn的转置矩阵；采用式(十四)计算该有效非零区段的中心垂直坐标Z_nkn，按式(十)计算水平坐标X_nkn：

$Z_{\mathrm{nkn}}=\frac{{\∫}_{m_{\mathrm{kns}}}^{m_{\mathrm{kne}}}m{I}^{\′}(m,n)\mathrm{dm}}{{\∫}_{m_{\mathrm{kns}}}^{m_{\mathrm{kne}}}{I}^{\′}(m,n)\mathrm{dm}}-\frac{M}{2}$    (十四)。