留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

基于2.5D标定扇的线阵相机内参数标定方法

张旭 毛庆洲 时春霖 胡庆武 金光 周昊 谢亦

张旭, 毛庆洲, 时春霖, 胡庆武, 金光, 周昊, 谢亦. 基于2.5D标定扇的线阵相机内参数标定方法[J]. 红外与激光工程, 2024, 53(4): 20230670. doi: 10.3788/IRLA20230670
引用本文: 张旭, 毛庆洲, 时春霖, 胡庆武, 金光, 周昊, 谢亦. 基于2.5D标定扇的线阵相机内参数标定方法[J]. 红外与激光工程, 2024, 53(4): 20230670. doi: 10.3788/IRLA20230670
Zhang Xu, Mao Qingzhou, Shi Chunlin, Hu Qingwu, Jin Guang, Zhou Hao, Xie Yi. Internal parameter calibration method of line-scan camera based on 2.5D calibration fan[J]. Infrared and Laser Engineering, 2024, 53(4): 20230670. doi: 10.3788/IRLA20230670
Citation: Zhang Xu, Mao Qingzhou, Shi Chunlin, Hu Qingwu, Jin Guang, Zhou Hao, Xie Yi. Internal parameter calibration method of line-scan camera based on 2.5D calibration fan[J]. Infrared and Laser Engineering, 2024, 53(4): 20230670. doi: 10.3788/IRLA20230670

基于2.5D标定扇的线阵相机内参数标定方法

doi: 10.3788/IRLA20230670
基金项目: 国家重点研发计划项目(2023YFC3009400, 2023YFB2603702)
详细信息
    作者简介:

    张旭,男,博士生,主要从事多传感器标定与数据融合方面的研究

  • 中图分类号: TP391.41

Internal parameter calibration method of line-scan camera based on 2.5D calibration fan

Funds: National Key Research and Development Program of China (2023YFC3009400, 2023YFB2603702)
  • 摘要: 针对工业生产线上或集成设备中线阵相机内参数定期检校困难且成本高的问题,提出一种基于2.5D标定扇的线阵相机内参数标定方法。设计了一款特定的2.5D标定扇,构建包含镜头畸变的线阵相机内参数标定模型,顾及了标靶相对相机的两个姿态角,通过方程线性变换法解算模型参数初值,利用改进的Levenberg-Marquardt (L-M)算法加速优化相机参数的求解过程,给出了详细的解算步骤和数据处理过程,并通过仿真和实测数据验证了方法的可行性。结果表明:该线阵相机标定方法简单灵活,可得到大量分布规律的特征点对,参数标定精度不受相机移动精度和特定移动方向的限制;在扇骨面与靶面间夹角小于10°时,可以获得高精度、高一致性的相机内参数,特征点重投影误差的最大值优于0.28 pixel,平均RMSE为0.11 pixel,标准差仅0.01 pixel。
  • 图  1  标定扇的设计

    Figure  1.  Design of the calibration fan

    图  2  采集标定扇的图像

    Figure  2.  Image acquisition of calibration fan

    图  3  线阵相机内参数标定模型

    Figure  3.  Internal parameter calibration model of line-scan camera

    图  4  参数解算与优化流程

    Figure  4.  The process of parameter calculation and optimization

    图  5  相机内参数误差的期望、方差与绝对误差

    Figure  5.  Expectation, variance and absolute error of camera intrinsic parameter error

    图  6  参数解不稳定的原因分析

    Figure  6.  Cause analysis of parameter solution instability

    图  7  实验场景

    Figure  7.  Test scene

    图  8  13根扇骨解算$ {y_0} $与$ f $的结果。(a) $ {y_0} $的解算结果;(b) $ f $的解算结果

    Figure  8.  The results of 13 fan bones solving $ {y_0} $ and $ f $. (a) The result of $ {y_0} $; (b) The result of $ f $

    图  9  重投影误差分布占比

    Figure  9.  The proportion of re-projection error distribution

    表  1  不同标定方法性能比较

    Table  1.   Performance comparison of different calibration methods

    MethodDimension of the
    calibration object
    Auxiliary equipmentAlgorithm difficultyRMSE of
    reprojection error
    Maximum reprojection
    error
    Contains lens
    distortion
    Horaud[6]2DPlatformLowNAN.a.No
    Luna[7]3DNoMediumNA0.71No
    Li[8]3DNoMedium0.3481.58Yes
    Sun[14]2DArea-array cameraHigh0.150NANo
    Song[9]3DNoHigh0.223N.a.No
    Yao[15]2DNoHigh0.460NANo
    Niu[16]3DNoMedium0.160-0.360NANo
    Su[17]3DNoMedium0.3771.76Yes
    Liao[10]3DNoMedium0.120-0.310NANo
    Hui[12]3DPlatformHigh0.310NAYes
    Hui[13]3DArea-array cameraHigh0.2200.36Yes
    Wang[18]2DLaserLow0.370.48No
    下载: 导出CSV

    表  2  仿真实验中相机的内参数和外参数

    Table  2.   The internal and external parameters of the camera in the simulation experiment

    Parameter Value Parameter Value
    y0/pixel 2048 $ f $/pixel 25/0.007
    $ a $ −1e−4 $ b $ 1e−10
    $ \alpha $/(°) 0 $ \varphi $/(°) 0
    $ {D_x} $/mm 1000 $ {D_y} $/mm 300
    下载: 导出CSV

    表  3  第7根扇骨的参数标定结果

    Table  3.   Parameter calibration results of the 7th fan bone

    Parameter Value Parameter Value
    y0/pixel 2055.35 $ f $/pixel 3571.62
    $ a $ 3.01e−5 $ b $ −1.67e−8
    $ \alpha $/(°) 2.80 $ \varphi $/(°) 1.32
    $ {D_x} $/mm 1449.50 $ {D_y} $/mm 58.68
    下载: 导出CSV

    表  4  重投影误差分析

    Table  4.   Analysis of reprojection error

    NumberMax/pixelMin/pixelRMSE/pixel
    40.1570.0020.086
    50.2410.0030.120
    60.1840.0260.104
    70.2620.0080.122
    80.2240.0270.124
    90.1940.0110.102
    100.2760.0020.129
    下载: 导出CSV

    表  5  参数优化方法对比

    Table  5.   Comparison of parameter optimization methods

    Number Standard L-M Improved L-M
    Times RMSE/pixel Times RMSE/pixel
    4 18 0.086 10 0.086
    5 21 0.122 12 0.120
    6 15 0.105 7 0.104
    7 16 0.126 8 0.122
    8 18 0.125 11 0.124
    9 19 0.102 8 0.102
    10 22 0.131 14 0.129
    下载: 导出CSV
  • [1] Wang Jing, Wei Liang, Xiang Wenhao, et al. High-precision camera calibration method considering projected circular edge blur and eccentricity error [J]. Infrared and Laser Engineering, 2021, 50(12): 20210130. (in Chinese)
    [2] Shi Yangqing, Chang Caixia, Liu Xiaohong, et al. Calibration methods and progress for internal and external parameters of area-array camera [J]. Laser & Optoelectronics Progress, 2021, 58(24): 2400001. (in Chinese)
    [3] Zhang Xu, Zhang Chao, Zhan Yinhu, et al. A method of removing abnormal data from linear array total station astronomical measurement based on star position prediction [J]. Acts Geodaetica et Cartographica Sinica, 2023, 52(4): 571-578. (in Chinese)
    [4] Wang Zichen, Dai Ming, Li Gang, et al. Calibration and analysis of line-scan camera based on virtual stereo pattern [J]. Infrared and Laser Engineering, 2014, 43(3): 902-907. (in Chinese)
    [5] Wang Guohui, Qian Kemao. Review on line-scan camera calibration methods [J]. Acta Optica Sinica, 2022, 40(1): 0111011. (in Chinese)
    [6] Horaud R, Mohr R, Lorecki B. On single-scanline camera calibration [J]. IEEE Transaction on Robotics and Automation, 1993, 9: 71-75. doi:  10.1109/70.210796
    [7] Luna C A, Mazo M, Lazaro J L, et al. Calibration of line-scan cameras [J]. IEEE Transactions on Instrumentation Measurement, 2010, 59: 2185-2190. doi:  10.1109/TIM.2009.2031344
    [8] Li Dongdong, Wen Gongjian, Hui Bingwei, et al. Cross-ratio invariant based line scan camera geometric calibration with static linear data [J]. Optics and Laser in Engineering, 2014, 62: 119-125. doi:  10.1016/j.optlaseng.2014.03.004
    [9] Song Kechen, Hou Bin, Niu Menghui, et al. Flexible line-scan camera calibration method using a coded eight trigrams pattern [J]. Optics and Laser in Engineering, 2018, 110: 296-307. doi:  10.1016/j.optlaseng.2018.06.014
    [10] Liao Ruiying, Zhu Jigui, Yang Linghui, et al. Flexible calibration method for line-scan cameras using a stereo target with hollow stripes [J]. Optics and Laser in Engineering, 2019, 113: 6-13.
    [11] Draréni J, Roy S, Sturm P. Plane-based calibration for linear cameras[J]. International Journal of Computer Vision , 2011, 91: 146-156.
    [12] Hui Bingwei, Zhong Jinrong, et al. Derermination of line scan camera parameters via the direct linear transformation[J]. Optics Engineering , 2012, 51(11): 113201.
    [13] Hui Binwei, Wen Gongjian, Zhong Peng. A novel line scan camera calibrations technique with an auxiliary frame camera [J]. IEEE Transactions on Instrumentation Measurement, 2013, 62: 2567-2575. doi:  10.1109/TIM.2013.2256815
    [14] Sun Bo, Zhu Jigui, Yang Linghui, et al. Stereo line-scan sensor calibration for 3D shape measurement [J]. Applic. Optics, 2017, 56: 7905-7913. doi:  10.1364/AO.56.007905
    [15] Yao Ming, Zhao Zuyun, Xu Bugao. Geometric calibration of line-scan camera using a planar pattern [J]. Journal of Electronic Imaging, 2014, 23(1): 013028.
    [16] Niu Menghui, Song Kechen, Wen Xin, et al. The line scan camera calibration based on space rings group [J]. IEEE Access, 2018, 6: 23711-23721. doi:  10.1109/ACCESS.2018.2817629
    [17] Su D, Bender A, Sukkarieh S. Improved cross-ratio invariant-based intrinsic calibration of a hyperspectral line-scan camera [J]. Sensors, 2018, 18: 1885-1896. doi:  10.3390/s18061885
    [18] Wang Lan, Sun Bo, Long Changyu, et al. A novel method for calibrating intrinsic parameters of linear array cameras vased on collinear feature points[J]. Infrared and Laser Engineering , 2015, 44(6): 1878-1883. (in Chinese)
    [19] Tsai R. A versatile camera calibration technique for high-accuracy 3D machine vision metrology using off-the-shelf TV cameras and lenses [J]. IEEE Journal of Robotics Automation , 1987, 3: 323-344.
    [20] Fang Suping, Xia Xiaohua, Xiao Yan, et al. Lens distortion calibration method for linear array cameras[J]. Journal of Xi’an Jiaotong University , 2013, 47(1): 11-14.
    [21] Shi Chunlin, Zhang Chao, Chen Changyuna, et al. One-dimensional maximum entropy image segmentation algorithm based on the small field of view of measuring robot star map[J]. Acta Geodaetica et Cartographica Sinica , 2018, 47: 446-454.
    [22] Yamashita N, Fukushima M. On the rate of convergence of the Levenberg-Marquardt method [J]. Computing, 2001, 15: 239-249.
    [23] Fan Jinyan. A modified Levenberg-Marquardt algorithm for singular system of nonlinear equations [J]. Journal of Computational and Applied Mathematics, 2003, 21: 625-636.
    [24] Amini K, Rostami F, Caristi G. An efficient Levenberg-Marquardt method with a new LM parameter for systems of nonlinear equations [J]. Optimization, 2018, 5: 637-650.
  • [1] 熊琨, 何旭然, 王春喜, 李佳宾, 杨长浩.  用于高精度准直测量的鱼眼相机标定方法 . 红外与激光工程, 2024, 53(2): 20230549-1-20230549-11. doi: 10.3788/IRLA20230549
    [2] 纪运景, 杜思月, 宋旸, 李振华.  基于线结构光旋转扫描和光条纹修复的三维视觉测量技术研究 . 红外与激光工程, 2022, 51(2): 20210894-1-20210894-9. doi: 10.3788/IRLA20210894
    [3] 王鹤, 李泽明.  激光测距仪与相机信息融合过程中位姿标定方法 . 红外与激光工程, 2020, 49(4): 0413002-0413002-8. doi: 10.3788/IRLA202049.0413002
    [4] 熊琨, 王春喜, 吴跃, 王锴磊, 姜云翔.  基于三轴转台的多视场星敏感器标定方法 . 红外与激光工程, 2019, 48(4): 417002-0417002(6). doi: 10.3788/IRLA201948.0417002
    [5] 杨杰, 黄见, 苑克娥, 孙培育, 刘琳琳, 胡顺星.  小型化大气二氧化碳垂直廓线测量系统标定 . 红外与激光工程, 2019, 48(5): 517004-0517004(7). doi: 10.3788/IRLA201948.0517004
    [6] 任同群, 黄海亭, 王晓东, 刘钰.  微装配中变焦变倍视觉系统标定及自动聚焦 . 红外与激光工程, 2018, 47(11): 1117001-1117001(10). doi: 10.3788/IRLA201847.1117001
    [7] 董春梅, 任顺清, 陈希军, 王常虹.  激光陀螺捷联惯导系统的模观测标定方法 . 红外与激光工程, 2018, 47(9): 917007-0917007(9). doi: 10.3788/IRLA201847.0917007
    [8] 王小冷, 杨凌辉, 林嘉睿, 任永杰, 尹彤.  基于周向约束定位原理的周向接收器标定方法 . 红外与激光工程, 2016, 45(7): 717004-0717004(6). doi: 10.3788/IRLA201645.0717004
    [9] 刘源, 谢睿达, 赵琳, 郝勇.  基于机器学习的大视场星敏感器畸变在轨标定技术 . 红外与激光工程, 2016, 45(12): 1217004-1217004(9). doi: 10.3788/IRLA201645.1217004
    [10] 王建军, 黄晨, 李舰艇.  空间目标红外辐射测量系统标定技术 . 红外与激光工程, 2016, 45(4): 404002-0404002(5). doi: 10.3788/IRLA201645.0404002
    [11] 胡冬梅, 刘泉, 于林韬, 朱一峰.  LSSVM模型下的LCVR相位延迟特性标定方法 . 红外与激光工程, 2016, 45(5): 517004-0517004(4). doi: 10.3788/IRLA201645.0517004
    [12] 刘海庆, 杨凌辉, 任永杰, 邾继贵.  基于正交柱面成像相机的大尺寸三维坐标测量 . 红外与激光工程, 2016, 45(11): 1117002-1117002(6). doi: 10.3788/IRLA201645.1117002
    [13] 毕超, 刘京亮, 刘飞, 吕来鹏, 雷雯铭.  激光测头的光束空间矢量标定方法 . 红外与激光工程, 2015, 44(12): 3775-3781.
    [14] 王澜, 孙博, 隆昌宇, 邹剑, 邾继贵.  一种基于共线特征点的线阵相机内参标定方法 . 红外与激光工程, 2015, 44(6): 1878-1883.
    [15] 崔凯, 刘朝晖, 李治国, 高雄, 谢友金, 王振宇.  坐标变换在空间望远镜误差标定中的应用 . 红外与激光工程, 2014, 43(S1): 151-156.
    [16] 王建军, 黄晨, 高昕, 李舰艇.  红外辐射测量系统内外标定技术 . 红外与激光工程, 2014, 43(6): 1767-1771.
    [17] 李斌, 吴海英, 王文涛, 刘鹏.  一种简易红外标定系统的设计 . 红外与激光工程, 2014, 43(2): 458-463.
    [18] 孙高飞, 张国玉, 高玉军, 王凌云, 苏拾, 付芸, 王向东.  星敏感器地面标定设备的设计 . 红外与激光工程, 2013, 42(10): 2822-2827.
    [19] 杨国洪, 李军, 韦敏习, 侯立飞, 易涛, 刘慎业.  X射线分光晶体高阶衍射效率标定方法 . 红外与激光工程, 2013, 42(12): 3325-3329.
    [20] 张晓娟, 张国玉, 孙高飞, 刘石, 郑茹, 高玉军.  大视场LCOS拼接星模拟器标定方法 . 红外与激光工程, 2013, 42(11): 2996-3001.
  • 加载中
图(9) / 表(5)
计量
  • 文章访问数:  13
  • HTML全文浏览量:  4
  • PDF下载量:  3
  • 被引次数: 0
出版历程
  • 收稿日期:  2023-12-01
  • 修回日期:  2024-01-30
  • 刊出日期:  2024-04-25

基于2.5D标定扇的线阵相机内参数标定方法

doi: 10.3788/IRLA20230670
    作者简介:

    张旭,男,博士生,主要从事多传感器标定与数据融合方面的研究

基金项目:  国家重点研发计划项目(2023YFC3009400, 2023YFB2603702)
  • 中图分类号: TP391.41

摘要: 针对工业生产线上或集成设备中线阵相机内参数定期检校困难且成本高的问题,提出一种基于2.5D标定扇的线阵相机内参数标定方法。设计了一款特定的2.5D标定扇,构建包含镜头畸变的线阵相机内参数标定模型,顾及了标靶相对相机的两个姿态角,通过方程线性变换法解算模型参数初值,利用改进的Levenberg-Marquardt (L-M)算法加速优化相机参数的求解过程,给出了详细的解算步骤和数据处理过程,并通过仿真和实测数据验证了方法的可行性。结果表明:该线阵相机标定方法简单灵活,可得到大量分布规律的特征点对,参数标定精度不受相机移动精度和特定移动方向的限制;在扇骨面与靶面间夹角小于10°时,可以获得高精度、高一致性的相机内参数,特征点重投影误差的最大值优于0.28 pixel,平均RMSE为0.11 pixel,标准差仅0.01 pixel。

English Abstract

    • 相机标定是视觉测量的基本工作,需预先建立物方世界坐标与像方平面坐标之间的几何模型,利用特征点对求解模型参数[1]。常用视觉测量相机包括面阵相机和线阵相机,面阵相机标定的研究已经比较成熟[2],但线阵相机通常仅对一维空间成像,所以面阵相机的标定方法无法直接采用,且关于线阵相机内参数标定的研究文献相对较少[34]

      目前线阵相机内参数的标定方法一般从设计标靶、参数标定模型、解算方法入手[5],总体可分为静态标定和动态标定两类,标定目标包括2D和3D两种。静态标定方面,Horaud等人最早开展关于线阵相机标定的研究,根据线阵相机的静态成像模型提出了一种基于交比不变性的标定方法。该方法计算相对简单,但其参数标定精度受移动平台的位移精度影响,且未顾及镜头非线性畸变的影响[6]。为了避免使用高精度的移动平台,Luna的方法使用了3D标定物,但该方法对3D标定物的制造精度要求较高,特征线数量有限且提取困难[7]。Li等人借鉴了张正友的面阵相机标定方法,克服了上述问题,但需要额外的面阵相机进行辅助标定,且标定精度取决于面阵相机的标定精度和测量精度,这使得标定过程更加复杂[8]。Song等人设计了一种基于八卦图形的三维标定物,理论上提高了标定精度,但该标定物需要高额成本以保证制造精度,且特征提取过程较为复杂[9]。Liao等人设计了一款3D空心条纹标定物,虽解决了偏心误差的问题,但标定物制造难度大、特征点数量少,仅限实验环境下使用[10]。关于动态标定的研究,Draréni等人首先提出了一种基于二维靶标面的动态标定方法,该方法特征点提取方便且标定精度较高,但该方法计算复杂度较高,并要求标定物在特定的方向上运动[11]。Hui等人先后提出了使用2D标定物和3D标定物的标定方法,不仅考虑了镜头畸变的影响而且解除了相机只能在特定方向移动的限制,但需建立复杂的动态标定模型,难以满足现场环境下标定[1213]表1对不同标定方法进行了统计和对比[5]。总体上,3D标定物制作难度和成本比2D标定物更高,但标定操作的复杂度更低;动态标定比静态标定的方法复杂,但高复杂性的算法往往可以换取较低的重投影误差。此外,为获得高精度的标定结果,可增加特征点数量,且标定模型应顾及镜头畸变的影响[5]

      表 1  不同标定方法性能比较

      Table 1.  Performance comparison of different calibration methods

      MethodDimension of the
      calibration object
      Auxiliary equipmentAlgorithm difficultyRMSE of
      reprojection error
      Maximum reprojection
      error
      Contains lens
      distortion
      Horaud[6]2DPlatformLowNAN.a.No
      Luna[7]3DNoMediumNA0.71No
      Li[8]3DNoMedium0.3481.58Yes
      Sun[14]2DArea-array cameraHigh0.150NANo
      Song[9]3DNoHigh0.223N.a.No
      Yao[15]2DNoHigh0.460NANo
      Niu[16]3DNoMedium0.160-0.360NANo
      Su[17]3DNoMedium0.3771.76Yes
      Liao[10]3DNoMedium0.120-0.310NANo
      Hui[12]3DPlatformHigh0.310NAYes
      Hui[13]3DArea-array cameraHigh0.2200.36Yes
      Wang[18]2DLaserLow0.370.48No

      现有方法均需在实验环境下完成线阵相机内参数标定,对已投入运行的集成设备中线阵相机很难完成内参数的定期检校。针对以上问题和实际需要,提出一种基于2.5D标定扇的线阵相机标定方法。首先,受中国传统折扇的启发,设计一款2.5D的简易标定扇,解决特征点数量少和易丢失的问题;其次,在此基础上构建不受线阵相机运动方向严格限制的参数标定模型,可解算出相机与靶面不平行、不垂直的小角度误差,并顾及了镜头畸变参数;然后,利用特征点对建立标定方程组,将模型线性化方程的解作为参数初值,再基于改进Levenberg-Marquardt (L-M)算法加速优化标定参数的求解过程。

    • 2.5D是一种介于2D和3D之间的表现形式,在平面上模拟3D效果。2.5D标靶既有3D立体测量特性,又有2D制作简单、成本低的优点。现有2D或3D标定物常以黑白线条、多边形图案等为特征,特征数量较少且缺乏关联性,部分标靶制作难度大、成本高,不便于携带和存放。在文献[18]的基础上,为获得相同基准且不同角度的多组特征点,设计了一款2.5D的标定扇,将特征点均匀布设在每根扇骨上,利用扇骨间夹角获得不同位置的多组特征点。每根扇骨的尺寸如图1(a)所示;图1(b)为标定扇立体图与直径1 mm的特征圆,其中有13根扇骨,相邻扇骨间夹角为3°,每个扇骨面上12个特征圆;图1(c)中以第7根扇骨面为$ O - XZ $面建立右手坐标系,$ \theta $为扇骨面与$ O - XZ $面的夹角;将各面上的特征亮圆正投影到$ O - XZ $面得到标定扇靶面,如图1(d)所示;将其1∶1打印在250 g铜版纸上,并平铺于光学平面上作为标靶,如图1(e)所示。

      图  1  标定扇的设计

      Figure 1.  Design of the calibration fan

      上述设计的2.5D标定扇具有2D标靶低成本、高精度的优点,也兼具3D标靶获取不同位置特征数据的可行性。需说明的是,靶面上特征点数量虽多(可继续增加),但其呈规律性分布,因此便于与像点配对。

    • 线阵相机仅包含一维成像单元,沿像元排列方向满足中心投影特征。将2.5D标定扇近似平行的放置在相机前方,相机随导轨移动采集标靶图像,如图2。在此基础上,图3描绘了特征点在图像坐标系$ (Y) $、相机坐标系$ ({X_c},{Y_c}) $、世界坐标系$ ({X_w},{Y_w},{Z_w}) $中的成像原理。$ {Y_w} $轴与$ {O_c} - {X_c}{Y_c} $共面,沿标定扇的靶面向上,$ {Z_w} $轴平行于运动方向,$ {O_w} $为过单根扇骨上所有圆心的直线与扇骨底端的交点,$ ({D_x},{D_y}) $为$ {O_w} $在相机坐标系$ ({X_c},{Y_c}) $中的坐标。$ P $为立体扇骨面的特征圆,$ {O_w}P = Y $。$ P' $为$ P $在标定扇靶面上的投影点,$ \theta $为$ P $点所在扇骨面与靶面的夹角。标定时尽量保证靶面与运动方向垂直,且平行于成像平面,但实际测量中存在误差角。$ \alpha $为垂直误差角,$ \varphi $为平行误差角。

      图  2  采集标定扇的图像

      Figure 2.  Image acquisition of calibration fan

      图  3  线阵相机内参数标定模型

      Figure 3.  Internal parameter calibration model of line-scan camera

      依据上述几何原理分析,可得线阵相机成像模型为:

      $$ \frac{{y - {y_0}}}{f} = - \frac{{{D_y} + Y\cos \alpha \cos \theta \cos \varphi }}{{{D_x} + Y\cos \alpha \cos \theta \sin \varphi }} $$ (1)

      式中:$ f = {{{f_v}} \mathord{\left/ {\vphantom {{{f_v}} d}} \right. } d} $;$ {f_v} $为镜头焦距,$ d $为像素尺寸;$ y $为特征点在图像坐标系的坐标;$ Y $为世界坐标系中$ {Y_w} $方向的坐标。

      克服镜头畸变是提高模型标定精度的重要环节。镜头畸变是指由于镜头的设计、制造和装配而导致的成像点偏离其理论位置的现象,包括由镜头径向曲率不足引起的径向畸变、由光心不共线引起的偏心畸变和由相机装配误差引起的薄棱镜畸变。线阵相机可看作一种特殊的面阵相机,镜头畸变模型可由面阵相机的畸变模型推导。面阵相机的总畸变表示为[19]

      $$ \left\{ {\begin{array}{*{20}{c}} {{\Delta _x}(x,y) = {k_1}x({x^2} + {y^2}) + {p_1}(3{x^2} + {y^2}) + 2{p_2}xy + {s_1}({x^2} + {y^2}) + O[{{(x,y)}^4}]} \\ {{\Delta _y}(x,y) = {k_1}y({x^2} + {y^2}) + {p_2}(3{y^2} + {x^2}) + 2{p_1}xy + {s_2}({x^2} + {y^2}) + O[{{(x,y)}^4}]} \end{array}} \right. $$ (2)

      式中:$ (x,y) $为像平面坐标系中点的坐标;$ {\Delta _x} $和$ {\Delta _y} $分别为畸变的$ x $分量和$ y $分量。$ {k_1} $为径向畸变系数;$ {p_1},{p_2} $为偏心畸变系数;$ {s_1},{s_2} $为薄棱镜畸变系数;$ O\left[ {{{(x,y)}^4}} \right] $为畸变的高阶无穷小量。畸变的高阶分量有时会导致解的不稳定和过大的计算量,而一般不被考虑[19]。此外,由于线阵传感器是一维的成像元件,仅需考虑$ y $分量。因此,线阵相机的镜头畸变模型可简化为[20]

      $$ {\Delta _y} = a{(y - {y_0})^2} + b{(y - {y_0})^3} $$ (3)

      式中:$ a = 3{p_1} + {s_1} $、$ b = {k_1} $为待标定的畸变分量。综上,线阵相机的几何标定模型为:

      $$ \begin{gathered} y - {y_0} + a{(y - {y_0})^2} + b{(y - {y_0})^3} + \\ f\frac{{{D_y} + Y\cos \alpha \cos \theta \cos \varphi }}{{{D_x} + Y\cos \alpha \cos \theta \sin \varphi }} = 0 \\ \end{gathered} $$ (4)

      式中:$ f,{y_0},a,b $为线阵相机的内参数;$ \varphi ,\alpha ,{D_x},{D_y} $为外参数;$ \theta $为已知值;$ y,Y $为特征观测值。

      标定扇上特征点分布规律,便于获取解算相机标定模型参数所需的特征点对$ Y,y $。单根扇骨的特征点为一组观测值,不同扇骨间的特征点竖向排列,易于区分,依据标定扇设计原理可得特征圆相对于扇骨底端(世界坐标系原点)的距离$ Y $。特征点在图像中灰度明显,借鉴文献[21]中基于一维最大熵的亚像素边缘提取算法便可得到特征圆的亚像素边缘点,进而可拟合出圆心像素坐标$ y $。鉴于标定扇靶面上的特征点及其像点间具有相同的分布规律,即可建立特征点及其像点的对应关系,进而得到多组特征观测值$ \left\{ {Y,y} \right\} $。不难发现,不同组观测值间等效于在线阵相机的物方空间放置立体标靶;特征点数量多且分布规律,避免了特征丢失的问题,有效减少了数据配对出错率。值得说明的是,像点$ y $通过拟合圆心得到,不受运动平台移动精度的影响,这使得参数标定精度同样不受其限制。

    • 借鉴文献[18]中线阵相机内参数解算方法,将非线性未知参数用一组线性系数替代,从而解算出方程中的参数,在此将该方法的模型解作为线阵相机内参数的初值。参数初值的具体解算过程如下:将特征观测值$ \left\{ {Y,y} \right\} $代入标定方程,依据线性变换理论解算模型参数初值。首先,将相机成像模型中8个非线性未知参数利用一组独立的线性系数代替,则公式(4)可重写为:

      $$ \begin{gathered} {k_1}y + {k_2}{y^2} + {k_3}{y^3} + {k_4}yY + {k_5}{y^2}Y + {k_6}Y + {k_7} = {y^3}Y \\ \end{gathered} $$ (5)

      其次,利用单组特征方程解算线性系数。若单根扇骨上有$ m(m \geqslant 7) $个特征点,将每个特征点的$ y,Y $代入公式(5),联立$ m $个特征方程解得$ {k_1},{k_2},{k_3},{k_4},{k_5},{k_6},{k_7} $。线性系数与未知参数之间的关系式为:

      $$ \left[ {\begin{array}{*{20}{c}} {{k_1}} \\ {{k_2}} \\ {{k_3}} \\ {{k_4}} \\ {{k_5}} \\ {{k_6}} \\ {{k_7}} \end{array} } \right] = \left[ {\begin{array}{*{20}{c}} { - \dfrac{{(1 - 2a{y_0} + 3by_0^2){D_x}}}{{b\cos \alpha \cos \theta \sin \varphi }}} \\ { - \dfrac{{\left( {a - 3b{y_0}} \right){D_x}}}{{b\cos \alpha \cos \theta \sin \varphi }}} \\ { - \dfrac{{{D_x}}}{{\cos \alpha \cos \theta \sin \varphi }}} \\ { - \dfrac{{1 - 2a{y_0} + 3by_0^2}}{b}} \\ { - \dfrac{{a - 3b{y_0}}}{b}} \\ { - \dfrac{{ - {y_0} + ay_0^2 - by_0^3}}{b} - \dfrac{{f\cos \varphi }}{{b\sin \varphi }}} \\ { - \dfrac{{\left( {ay_0^2 - by_0^3 - {y_0}} \right){D_x} + f{D_y}}}{{b\cos \alpha \cos \theta \sin \varphi }}} \end{array}} \right] $$ (6)

      最后,利用多组特征方程的线性解计算模型参数初值。相机采集靶面特征时各扇骨的世界坐标系原点近似不变,仅改变了扇骨面与靶面间的旋转角,即$ f,{y_0},a,b,\alpha ,\varphi ,{D_x},{D_y} $不变,$ \theta $改变。每根扇骨可以解出一组$ [{k_{i1}},{k_{i2}},{k_{i3}},{k_{i4}},{k_{i5}},{k_{i6}},{k_{i7}}] $,${i = 1,2,\cdots,13}$。联立每根扇骨成像模型中的$ {k_{i3}} $和$ {k_{i6}} $,消去$ \varphi $可得公式(7),将$ {k_{i3}} $代入$ {k_{i7}} $可得公式(8):

      $$ \begin{gathered} D_x^2{[b{D_x}{k_{i6}} + ( - {y_0} + ay_0^2 - by_0^3)]^2} + D_x^2{f^2} = {({k_{i3}}f\cos \alpha \cos \theta )^2} \\ \end{gathered} $$ (7)
      $$ {k_{i7}}b{D_x} = {D_x}{k_{i3}}( - {y_0} + ay_0^2 - by_0^3) - f{D_y}{k_{i3}} $$ (8)

      将$ [{k_{i3}},{k_{i6}},{k_{i7}}] $,${i = 1,2,\cdots,13}$代入公式(7)和公式(8),建立26个方程,联系方程组,解出7个参数:$ f,{y_0},a,b,\alpha ,{D_x},{D_y} $。将$ f,{y_0},a,b,\alpha ,{D_x},{D_y} $代入每个扇骨的成像模型(公式(4)),可解出$ \varphi $。经上述过程可得标定模型各参数的初值。

    • 为进一步提高线阵相机标定模型参数的解算精度,避免解的过程中方程出现奇异问题,提出一种基于改进的L-M算法对参数初值进行优化。L-M算法结合了高斯-牛顿算法和梯度下降算法的优点,能够避免由于Hessian矩阵奇异而导致无法继续迭代的问题,常用于求解非线性方程。首先,将每根扇骨的特征点对$ \left\{ {y,Y} \right\} $和模型参数初值代入公式(4)中,并改写成$ F(x) = 0 $,一阶泰勒展开式为$ F(x + \Delta x) = F(x) + F'(x)\Delta x $,并组成标定方程组,从而将参数解算问题转化为计算函数$ G(x) $最小值的问题。

      $$ \begin{gathered} G(x) = \mathop {\min }\limits_{\Delta x} \left\{ {\frac{1}{2}{{(F({x_n}) + J\Delta x)}^{\mathrm{T}}} \cdot \left(F({x_n}) + J\Delta x + \frac{1}{2}\lambda \Delta {x^{\mathrm{T}}}\Delta x\right)} \right\} \\ \end{gathered} $$ (9)

      式中:$ x = \{ f,{y_0},a,b,\alpha ,{D_x},{D_y},\varphi \} $;$ J = F'(x) $。函数$ G(x) $关于$ \Delta x $的一阶导数为0,得:

      $$ \Delta x = - {({J^{\mathrm{T}}}J + \lambda I)^{ - 1}}{J^{\mathrm{T}}}F({x_n}) $$ (10)

      式中:$ \lambda $为信赖域半径;$ J $为雅克比矩阵,由特征点模型方程中每个参数的偏导数组成。

      $$ J = \left| {\begin{array}{*{20}{c}} {\dfrac{{\partial {F_1}(x)}}{{\partial {y_0}}}}&{\dfrac{{\partial {F_1}(x)}}{{\partial f}}}&{\dfrac{{\partial {F_1}(x)}}{{\partial a}}}&{\dfrac{{\partial {F_1}(x)}}{{\partial b}}} \\ {\dfrac{{\partial {F_2}(x)}}{{\partial {y_0}}}}&{\dfrac{{\partial {F_2}(x)}}{{\partial f}}}&{\dfrac{{\partial {F_2}(x)}}{{\partial a}}}&{\dfrac{{\partial {F_2}(x)}}{{\partial b}}} \\ \vdots & \vdots & \vdots & \vdots \\ {\dfrac{{\partial {F_{12}}(x)}}{{\partial {y_0}}}}&{\dfrac{{\partial {F_{12}}(x)}}{{\partial f}}}&{\dfrac{{\partial {F_{12}}(x)}}{{\partial a}}}&{\dfrac{{\partial {F_{12}}(x)}}{{\partial b}}} \end{array}} \right.\left. {\begin{array}{*{20}{c}} {\dfrac{{\partial {F_1}(x)}}{{\partial \alpha }}}&{\dfrac{{\partial {F_1}(x)}}{{\partial {D_x}}}}&{\dfrac{{\partial {F_1}(x)}}{{\partial {D_y}}}}&{\dfrac{{\partial {F_1}(x)}}{{\partial \varphi }}} \\ {\dfrac{{\partial {F_2}(x)}}{{\partial \alpha }}}&{\dfrac{{\partial {F_2}(x)}}{{\partial {D_x}}}}&{\dfrac{{\partial {F_2}(x)}}{{\partial {D_y}}}}&{\dfrac{{\partial {F_2}(x)}}{{\partial \varphi }}} \\ \vdots & \vdots & \vdots & \vdots \\ {\dfrac{{\partial {F_{12}}(x)}}{{\partial \alpha }}}&{\dfrac{{\partial {F_{12}}(x)}}{{\partial {D_x}}}}&{\dfrac{{\partial {F_{12}}(x)}}{{\partial {D_y}}}}&{\dfrac{{\partial {F_{12}}(x)}}{{\partial \varphi }}} \end{array}} \right| $$ (11)

      Yamashita证实$ {\lambda _n} = {\left\| {{F_n}} \right\|^2} $时,在一定误差范围内L-M算法具有二次收敛性[22];Fan提出$ {\lambda _n} = {\mu _n}\left\| {{F_n}} \right\| $,同样满足类似的收敛性,可逆转初始值远离解时收敛效果差的问题[23];Keyvan提出$ {\lambda _n} = \dfrac{{{\mu _n}\left\| {{F_n}} \right\|}}{{1 + \left\| {{F_n}} \right\|}} $,证明在解的一定误差范围内具有全局收敛和二次收敛性[24]。在此将控制步长设为:

      $$ {\lambda _n} = \frac{{{\mu _n}{{\left\| {{F_n}} \right\|}^2}}}{{1 + {{\left\| {{F_n}} \right\|}^2}}} $$ (12)

      $ {\mu _n} $的初值$ {\mu _0} $常取$ {J^{\mathrm{T}}}J $对角线元素的最大值。当结果远离解时,$ {\lambda _n} $接近标准L-M算法中使用的$ {\mu _n} $;当结果接近解时,$ {\lambda _n} $接近$ {\mu _n}{\left\| {{F_n}} \right\|^2} $,与文献[23]类似。$ \left\| {{F_n}} \right\| $取平方可加速迭代计算,提高计算效率。利用下降比率$ {r_n} $作为判断算法参数更新的依据,将每次计算的新点与之前步骤中最差的点进行比较,$ {r_n} $的定义为:

      $$ {r_n} = \frac{{\mathop {\max }\limits_{j = 1,\cdots,n} \left\{ {{{\left\| {{F_j}} \right\|}^2}} \right\} - {{\left\| {F({x_n} + \Delta {x_n})} \right\|}^2}}}{{{{\left\| {{F_n}} \right\|}^2} - {{\left\| {{F_n} + {J_n}\Delta {x_n}} \right\|}^2}}} $$ (13)

      通过公式(14)对$ {r_n} $判断,更新算法参数。

      $$\begin{split} & {x_{n + 1}} = \left\{ {\begin{array}{*{20}{c}} {{x_n} + \Delta {x_n}}&{{r_n} \geqslant {p_0}} \\ {{x_n}}&{{\rm{otherwise}}} \end{array}} \right.\\& {\mu _{n + 1}} = \left\{ {\begin{array}{*{20}{c}} {4{\mu _n}}&{{r_n} < {p_1}} \\ {{\mu _n}}&{{p_2} < {r_n} < {p_1}} \\ {\max \left\{ \dfrac{{{\mu _n}}}{4},m\right\} }&{{\rm{otherwise}}} \end{array}} \right. \end{split} $$ (14)

      式中:$ m,{p_0},{p_1},{p_2} $为常数,$ 0 < m \leqslant {\mu _n} $,$ 0 < {p_0} \leqslant {p_1} \leqslant {p_2} < 1 $。在迭代计算过程中设置了与梯度向量有关的终止条件和与$ x $有关的终止条件,如公式(15)所示,二者满足其一即终止计算。

      $$ \left\{ {\begin{array}{*{20}{l}} {\left\| {J_n^{\mathrm{T}}{F_n}} \right\| < {\varepsilon _1}} \\ {\left\| {{x_{n + 1}} - {x_n}} \right\| < {\varepsilon _2}(\left\| {{x_n}} \right\| + {\varepsilon _2})} \end{array}} \right. $$ (15)

      式中:$ {\varepsilon _1} $和$ {\varepsilon _2} $为迭代终止的条件参数。

      综上所述,线阵相机标定参数的解算和优化流程如图4所示。

      图  4  参数解算与优化流程

      Figure 4.  The process of parameter calculation and optimization

    • 标定物设计的合理性与观测值$ \left\{ {y,Y} \right\} $的获取精度直接影响参数标定结果的精度。$ Y $依据标定扇设计原理获得,无计算误差引入,而$ y $的误差来源于特征中心提取误差,是随机误差;同时,需探究标定扇骨面与靶面间旋转角$ \theta $的合理范围。为此,在Windows平台下利用编程软件仿真线阵相机标定数据,并分析相关因素对标定结果的影响。仿真实验中相机内外参数见表2

      表 2  仿真实验中相机的内参数和外参数

      Table 2.  The internal and external parameters of the camera in the simulation experiment

      Parameter Value Parameter Value
      y0/pixel 2048 $ f $/pixel 25/0.007
      $ a $ −1e−4 $ b $ 1e−10
      $ \alpha $/(°) 0 $ \varphi $/(°) 0
      $ {D_x} $/mm 1000 $ {D_y} $/mm 300

      为测试像点误差量级与扇骨面相对于靶面的旋转角$ \theta $对标定结果的影响,现设计仿真实验如下。首先,令$ \theta $=0,利用相机的成像模型和内外参数生成50个成像特征点坐标和对应的像点坐标;而后,改变$ \theta $的同时向像点坐标加入不同量级的随机误差,其中,旋转角$ \theta $以1°为增量从0°增加到20°,而对每个$ \theta $,像点坐标$ y $加入正态分布的误差$ N(0,{\sigma ^2}) $,方差$ \sigma $以0.01 pixel为增量从0 pixel增加到2 pixel,由此改变像点误差;对于每一组$ \theta $和$ \sigma $,进行1000次独立仿真实验,图5给出了基于文中标定方法解算相机内参数$ f $和$ {y_0} $的期望、方差以及绝对误差平均值,可以看出:

      图  5  相机内参数误差的期望、方差与绝对误差

      Figure 5.  Expectation, variance and absolute error of camera intrinsic parameter error

      1)随着$ y $误差等级的增加,$ f $的绝对误差逐渐变大,期望变化的幅度和方差均增大,即像点误差越大,$ f $越敏感;

      2)当$ y $的误差在0~0.6 pixel时,随着$ \theta $的增大,$ f $无明显变化;而$ y $误差继续增大时,增大$ \theta $,$ f $的绝对误差、期望、方差也随之增加;

      3)当$ \theta $在0°~10°时,增大$ y $的误差对$ {y_0} $的影响不明显;继续增大$ \theta $时,$ y $误差的增大会导致$ {y_0} $的绝对误差、期望、方差急剧增加,标定结果也越不稳定。

      值得说明的是,目前亚像素高精度定位算法的定位精度小于0.2 pixel,因此实验中特征提取误差对参数解算精度影响较小。此外,标定扇选用$ \theta \in [ - {10^ \circ },{10^ \circ }] $的扇骨面有利于得到精确、稳定的参数解,这是由于随着$ \theta $增大,相同像素误差$ \Delta y $所引起的数据误差$ \Delta Y $增大,从而引起特征点对$ y,Y $的精度下降,致使方程组解的稳定性降低,如图6所示。此外,在参数解算过程中,标定方程的高度奇异性也可能是标定误差被放大的原因。

      图  6  参数解不稳定的原因分析

      Figure 6.  Cause analysis of parameter solution instability

    • 实验选用武汉汉宁轨道公司研发的钢轨测量模块中的线阵相机完成内参数标定,相机有效像元4096个,像元尺寸7 μm,镜头焦距为25 mm,相机视场明亮,光照强度均匀。图7(a)中,相机模块于导轨上匀速运动(速度为6 mm/s),距离编码器以0.2 mm/次的频率向相机发送触发采集的信号,使单个像素的分辨率小于特征圆的半径,以确保采集标靶图像中的特征圆有效。相机沿导轨运动方向采集标定扇的二维图像如图7(b)所示,其中有13列像点,每列代表单根扇骨的特征圆。运用亚像素高精度定位算法处理图像,得到各特征点的像点坐标$ y $,根据提出的标定方法完成相机标定,得到每根扇骨的参数解算结果。利用第7组数据解算的线阵相机内参数见表3

      图  7  实验场景

      Figure 7.  Test scene

      表 3  第7根扇骨的参数标定结果

      Table 3.  Parameter calibration results of the 7th fan bone

      Parameter Value Parameter Value
      y0/pixel 2055.35 $ f $/pixel 3571.62
      $ a $ 3.01e−5 $ b $ −1.67e−8
      $ \alpha $/(°) 2.80 $ \varphi $/(°) 1.32
      $ {D_x} $/mm 1449.50 $ {D_y} $/mm 58.68

      扇骨面与靶面的旋转角$ \theta $在0°~±10°范围内时参数解算结果较稳定。图8所示为各扇骨解得的像主点坐标$ {y_0} $与焦距$ f $。其中,由第4~10根扇骨解算的相机内参数结果具有一致性,且各扇骨的$ \theta \in \left[ { - {{10}^ \circ },{{10}^ \circ }} \right] $。当扇骨面与靶面的夹角$ \left| \theta \right| > {10^ \circ } $时,参数解算结果不稳定,如第1、2、3、11、12、13扇骨,像主点和焦距的解算结果分别浮动约5 pixel和25 μm,且有随着$ \theta $增大而浮动范围增加的趋势,这与图5所得结论保持一致。

      图  8  13根扇骨解算$ {y_0} $与$ f $的结果。(a) $ {y_0} $的解算结果;(b) $ f $的解算结果

      Figure 8.  The results of 13 fan bones solving $ {y_0} $ and $ f $. (a) The result of $ {y_0} $; (b) The result of $ f $

      采用上述标定方法获得的稳定解可作为集成模块中线阵相机的准确内参数。为验证模型参数的正确性,将每根扇骨标定出的相机参数代入公式(4),将标定特征点重投影到像面,原像点与重投影点之差记为重投影误差$ {\rho _c} $。第4~10根扇骨共156个特征点,用单根扇骨上12个特征点解算相机参数,将其他特征点作为检验点,计算重投影误差以验证参数的标定精度。表4给出了每根扇骨所得重投影误差的最大值、最小值、均方根误差,图9统计了所有特征点重投影误差$ {\rho _c} $的分布占比,其中89%特征点的重投影误差小于0.20 pixel,最大误差优于0.28 pixel,单根扇骨上特征点重投影误差的RMSE为0.086~0.129 pixel。需要说明的是,采集10次标定扇图像完成相机内参数求解与精度验证,其中重投影误差的平均RMSE为0.11 pixel,标准差为0.01 pixel,参数标定结果精度高且一致性好,相较于文献[18](重投影精度为0.37 pixel),重投影精度提高了69.7%。

      表 4  重投影误差分析

      Table 4.  Analysis of reprojection error

      NumberMax/pixelMin/pixelRMSE/pixel
      40.1570.0020.086
      50.2410.0030.120
      60.1840.0260.104
      70.2620.0080.122
      80.2240.0270.124
      90.1940.0110.102
      100.2760.0020.129

      图  9  重投影误差分布占比

      Figure 9.  The proportion of re-projection error distribution

      另外,基于改进的L-M算法在不影响内参数标定精度的同时可有效提高计算效率。分别利用标准L-M算法和改进L-M算法优化第4~10根扇骨的参数初值,表5对比了两种优化方法的迭代次数和重投影误差的RMSE。相比于前者,改进的L-M算法将迭代次数减少约50%,而参数标定结果无明显差异。

      表 5  参数优化方法对比

      Table 5.  Comparison of parameter optimization methods

      Number Standard L-M Improved L-M
      Times RMSE/pixel Times RMSE/pixel
      4 18 0.086 10 0.086
      5 21 0.122 12 0.120
      6 15 0.105 7 0.104
      7 16 0.126 8 0.122
      8 18 0.125 11 0.124
      9 19 0.102 8 0.102
      10 22 0.131 14 0.129
    • 线阵相机常与其他传感器以模块化的形式安置于设备内部,对于这类相机内参数的定期检校成本高,困难程度大。鉴于此,提出一种基于2.5D标定扇的线阵相机内参标定方法。一方面,2.5D标定扇兼具3D标靶的立体测量效果和2D标靶低成本、高精度的优点,特征点数量多分布规律,避免了特征易丢失的问题,且特征点与像点易配对。另一方面,所构建的线阵相机标定模型顾及了镜头畸变误差和靶面与像面间两个方向的小角度姿态,相机移动方向与标定物的位置无严格限定。实验表明,采用$ \theta $<±10°的扇骨标定相机内参为最佳,标定结果精度高且一致性好,其中89%特征点的重投影误差小于0.20 pixel,最大误差优于0.28 pixel,平均RMSE为0.112 pixel。此外,相比于标准L-M算法,改进的L-M算法在不影响参数优化精度的同时将迭代次数减少一半,同时避免方程在解的过程中出现奇异。

      值得说明的是,该方法还适用于钢轨病害检测、车间质检、文物数字化等应用场景的线阵相机完成定期内参数检校,有利于传感器维持高精度测量状态,同时亦可为其他可见光、红外等线性光电探测器测量装置的内参数标定提供有益参考。然而,该标定方法也有局限性,要求线阵相机与标定靶面间应具有较小的姿态角以确保解算线性系数方程时$ {D_x} $、$ {D_y} $近似不变,才能获得准确的参数初值,否则会引起解的不稳定。

参考文献 (24)

目录

    /

    返回文章
    返回