本发明属于海洋测绘技术领域，具体涉及一种机载lidar测深海气界面折射改正方法。
背景技术：
机载lidar测深系统具有测量精度高、测点密度高，工作效率高、机动性强、测量连续性等特点，特别适合于浅水区、海岛礁附近区域等复杂地形的快速探测，可以实现海岸线水上水下地形的无缝拼接。对于满足海洋、测绘、水利、交通、外交、海军等多方面的迫切需求，具有重要意义。
由于波浪、潮汐等外界因素的影响，海面不断地发生随机起伏，改变了海气界面的几何形状。当机载lidar测深系统发射的蓝绿激光穿过海气界面时，瞬时海面倾斜会造成折射光线路径发生偏离，严重影响海底地形测量精度。目前大部分折射改正主要依靠扩大激光脉冲发散角覆盖多个波浪来减小折射影响的方式，该改正方法是基于激光扩束以平滑瞬时海面实现的，然而较大的光斑自身就降低了测深精度，已不适用于当前主流的小光斑机载lidar测深系统。
因此，有必要提出一种机载lidar测深海气界面折射改正方法，以提高机载lidar测深精度。
技术实现要素：
针对现有技术中存在的上述技术问题，本发明提出了一种机载lidar测深海气界面折射改正方法，设计合理，解决了瞬时海面倾斜造成的折射光线路径偏移问题，有效提高了机载lidar测深精度。
为了实现上述目的，本发明采用如下技术方案：
一种机载lidar测深海气界面折射改正方法，包括以下步骤：
步骤1：计算瞬时海面激光点坐标(x,y,z)；
步骤2：基于最小二乘和波浪谱理论，利用步骤1中计算的瞬时海面激光点坐标，构建瞬时三维海面模型；具体包括如下步骤：
步骤2.1：瞬时三维海面模型函数的确定，其能够看成是连续的傅里叶变换频谱的叠加，且为连续可导曲面；所述瞬时三维海面模型函数为：
其中，(x,y,z)为瞬时海面三维坐标，ζi、wi、t、αi和εi分别代表波浪的振幅、角频率、时间、方向角和初相；g＝9.8m/s2；m为余弦波叠加次数(m＝30)；
步骤2.2：针对步骤2.1中的瞬时三维海面模型函数，确定其初始参数，根据jonswap谱，将公式(1)中的角频率ωi设为从1.1到4.0，且角频率间隔为0.1，i＝1,2,3,...,30；同时，计算出30个相应的初始振幅αi；30个方向角和30个初相设置为零；
步骤2.3：利用levenberg-marquardt算法进行瞬时三维海面拟合，获得瞬时三维海面模型；
步骤3：根据步骤2中构建的瞬时三维海面模型，根据公式(2)计算瞬时海面斜率；通过计算公式(1)的偏微分，得到公式(2)；
其中，
步骤4：以步骤2中构建的瞬时三维海面模型作为光线追踪平台，追踪穿过海气界面的每束激光，计算折射光线单位向量；
步骤5：根据步骤4中获得的折射光线单位向量，对海底激光点的点位偏差和测深误差进行改正。
优选地，在步骤1中，具体包括如下步骤：
步骤1.1：获取机载lidar测深数据，读取测深激光数据中的全波形数据，并对波形数据进行去噪和拟合，确定出测深波形；
步骤1.2：结合发射激光天顶角、载体瞬时位置/姿态信息，计算瞬时海面激光点坐标(x,y,z)。
优选地，在步骤4中，海气界面光线追踪的具体方法如下：
步骤4.1：基于步骤3中获得的海面斜率
其中，(nx,ny,nz)分别为单位法向量n在x、y、z方向上的三个分量；
步骤4.2：基于步骤1获得的瞬时海面三维坐标(x,y,z)，根据公式(4)，计算入射光线的单位向量ilaser：
步骤4.3：结合步骤4.1获得的折射界面单位法向量n和步骤4.2获得的入射光线单位向量ilaser，通过余弦定理，根据公式(5)，计算激光入射角ilaser：
其中，入射角ilaser为锐角；
步骤4.4：计算激光折射角rlaser，折射角rlaser为锐角；有两种计算方法：
第一求解方法，根据snell法则，激光折射角rlaser表示如公式(6)所示：
第二求解方法，由于激光折射角rlaser是由折射光线单位向量rlaser和法向量n交汇形成的锐角，因此利用公式(7)计算折射角rlaser：
步骤4.5：计算折射光线单位向量rlaser(rx,ry,rz)；由于折射光线单位向量rlaser有三个未知量，进行解算需要构建三个联立方程：
方程一：将步骤4.4中的公式(6)和(7)带入sin2rlaser+cos2rlaser＝1，并兼顾公式(5)，得到目标方程，如公式(8)所示：
方程二：考虑到rlaser为单位向量，向量rlaser表示如公式(9)所示：
方程三：因为向量n、ilaser、rlaser在同一平面内，所以三个向量满足共面方程，如公式(10)所示：
ilaser×n·rlaser＝0(10)；
通过构建三个联立方程(8)、(9)和(10)，计算折射光线单位向量rlaser(rx,ry,rz)。
优选地，在步骤5中，改正海底激光点的点位偏差和测深误差的具体方法如下：
步骤5.1：根据公式(11)，对海底激光点的点位m(xm,ym,zm)进行折射改正：
其中，na为空气折射率，na＝1，nw为海水折射率，nw＝1.33，δt为激光从海面到海底的旅行时；
步骤5.2：根据公式(12)，对海底激光点水深值d进行折射改正：
本发明所带来的有益技术效果：
本发明提出了一种机载lidar测深海气界面折射改正方法，与现有技术相比，本发明方法考虑到瞬时海面倾斜造成的折射光线路径偏离变化，并通过构建的瞬时海面模型作为光线追踪平台，追踪每一束激光脉冲，并对海底激光点的平面位置和水深提供有效的折射改正，解决了瞬时海面倾斜造成的折射光线路径偏移问题，有效提高了机载lidar测深精度。
附图说明
图1为本发明一种机载lidar测深海气界面折射改正方法的流程图。
图2为机载lidar测深过程中瞬时海面倾斜造成折射光线路径偏移的示意图。
图3a为本发明中机载lidar测深系统获得的瞬时海面激光点云图。
图3b为本发明中构建的瞬时三维海面模型图。
图4为本发明中海面斜率的示意图。
图5为本发明中机载lidar测深光线追踪的示意图。
具体实施方式
下面结合附图以及具体实施方式对本发明作进一步详细说明：
本发明提供了一种机载lidar测深海气界面折射改正方法，其流程如图1所示。机载lidar测深过程中瞬时海面倾斜造成的折射光线路径偏移，如图2所示。改正方法具体包括以下步骤：
步骤1：获取机载lidar测深数据，计算出瞬时海面激光点坐标；
具体的，采用机载lidar测深系统获取测深激光数据。其中有效的测深激光点是指去除陆地上的激光点，只留下水面上的激光点。
进一步的实施例中，步骤1具体包括如下步骤：
步骤1.1：获取机载lidar测深数据，读取测深激光数据中的全波形数据，并对波形数据进行去噪和拟合，提取出测深波形；
步骤1.2：结合发射激光天底角、载体(飞机)瞬时位置/姿态信息计算瞬时海面激光点坐标。
具体实施时，由于海气界面激光点的折射改正是在地心地固坐标系进行的，因此，需要进行一系列的坐标转换，最终得到海面激光点在地心地固坐标系下的坐标，如图3a所示，为获得的100m×100m格网区域的瞬时海面激光点云。
步骤2：基于最小二乘和波浪谱理论，利用步骤1中计算的瞬时海面激光点云数据构建瞬时三维海面模型。
具体的，以海洋学中的海浪谱理论为研究基础，瞬时海面模型可以看成是由许多个不同波高、不同周期、不同相位、不同运动方向的余弦叠加，通过拟合532nm蓝绿色激光获取的瞬时海面激光点云，可以计算每个余弦函数的参数。从而构建瞬时三维海面模型，准确地反映机载lidar测深过程中的实际海气界面变化情况，为折射改正提供光线追踪平台。
进一步的实施例中，步骤2具体包括如下步骤：
步骤2.1：瞬时三维海面模型函数的确定，其可以看成是连续的傅里叶变换频谱的叠加，且为连续可导曲面；所述瞬时三维海面模型函数，如公式(1)所示：
其中，(x,y,z)为瞬时海面三维坐标，ζi、wi、t、αi和εi分别代表波浪的振幅、角频率、时间、方向角和初相；g＝9.8m/s2；m为余弦波叠加次数(m＝30)；
步骤2.2：针对步骤2.1中的瞬时三维海面模型函数，确定其初始参数。根据jonswap谱，公式(1)中角频率ωi设为从1.1到4.0，且角频率间隔为0.1，i＝1,2,3,...,30；同时，计算出30个相应的初始振幅αi；30个方向角和30个初相设置为零；
步骤2.3：利用levenberg-marquardt算法进行瞬时三维海面拟合，获得瞬时三维海面模型。
具体实施时，由于海面是瞬时变化的，因此选取的瞬时海面激光点云应为同一个航带，这样可以避免不同航带测量时间内瞬时海面的变化。如图3b所示，为利用图3a中的瞬时海面激光点云构建的瞬时三维海面模型。
步骤3：根据步骤2中构建的瞬时三维海面模型计算瞬时海面斜率。
具体实施时，计算瞬时海面斜率需要获得瞬时海面倾角。如图4所示，海面倾角
进一步的实施例中，步骤3中瞬时海面斜率的计算基于公式(2)：
其中，
步骤4：以步骤2中构建的瞬时三维海面模型作为光线追踪平台，追踪穿过海气界面的每束激光，计算折射光线单位向量。
具体实施时，需要确定机载lidar测深光线追踪路线，如图5所示。n为反射面c的单位法向量，ilaser为海面上入射光线的单位向量，rlaser为水体中折射光线的单位向量。在以上3个向量中，n和ilaser分别可以通过海面斜率和激光点p的地理坐标获得。rlaser(rx,ry,rz)则需要解算3个联立方程组获得。
进一步的实施例中，步骤4具体包括如下步骤：
步骤4.1：基于步骤3中获得的海面斜率
其中，(nx,ny,nz)分别为单位法向量n在x、y、z方向上的三个分量；
步骤4.2：基于步骤1获得的瞬时海面三维坐标(x,y,z)，根据公式(4)，计算入射光线的单位向量ilaser：
步骤4.3：结合步骤4.1获得的折射界面单位法向量n和步骤4.2获得的入射光线单位向量ilaser，通过余弦定理，根据公式(5)，计算激光入射角ilaser(为锐角，小于90°)：
步骤4.4：计算激光折射角rlaser(为锐角，小于90°)，有两种计算方法：
第一求解方法，根据snell法则，激光折射角rlaser表示如公式(6)所示：
第二求解方法，顾及激光折射角rlaser是由折射光线单位向量rlaser和法向量n交汇形成的锐角，也可以利用公式(7)计算折射角rlaser：
步骤4.5：计算折射光线单位向量rlaser(rx,ry,rz)。由于折射光线单位向量rlaser有三个未知量，进行解算需要构建三个联立方程：
方程一：将步骤4.4中的公式(6)和(7)带入sin2rlaser+cos2rlaser＝1，并兼顾公式(5)，得到目标方程，如公式(8)所示：
方程二：考虑到rlaser为单位向量，向量rlaser表示如公式(9)所示：
方程三：因为向量n、ilaser、rlaser在同一平面内，所以三个向量满足共面方程，如公式(10)所示：
ilaser×n·rlaser＝0(10)
通过构建三个联立方程(8)、(9)和(10)，计算折射光线单位向量rlaser(rx,ry,rz)。
步骤5：根据步骤4中获得的折射光线单位向量，对海底激光点的点位偏差和测深误差进行改正。
进一步的实施例中，步骤5具体包括步骤：
步骤5.1：根据公式(11)，对海底激光点的点位m(xm,ym,zm)进行折射改正：
其中，na为空气折射率(na＝1)，nw为海水折射率(nw＝1.33)，δt为激光从海面到海底的旅行时；
步骤5.2：根据公式(12)，对海底激光点水深值d进行折射改正：
综上所述，本发明提供了一种机载lidar测深海气界面折射改正方法。考虑到瞬时海面倾斜造成的折射光线路径偏离变化，方法包括：通过构建的瞬时海面模型作为光线追踪平台，追踪每一束激光脉冲，并对海底激光点的平面位置和水深提供有效的折射改正。本发明解决了瞬时海面倾斜造成的折射光线路径偏移问题，有效提高了机载lidar测深精度。
当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。