ABSTRACT

The traditional bathymetric survey involves a large number of discrete, irregularly distributed water depth data. In order to solve the problem of complex calculation and slow speed of traditional fitting method, the idea of array algebra is applied to the fitting of the sounding grid. The least squares solution of the array algebra is established by taking the biquadratic polynomial fitting as an example, and the least square solution of the quadratic matrix equation with the help of the normal equation is given. Moreover, the paper gives the basic steps to achieve the improved quasi legal, efficiency and accuracy of the experimental data to determine unknown point water depth on the inverse distance weighted method, global polynomial interpolation and least square fitting method of array algebra is compared. The results show that the efficiency of array algebra fitting is better than that of the other two methods, and can be applied to the data processing of hydrology survey.

Keywords:Array Algebra, Hydrology Survey, Biquadratic, Least Square Fitting

基于阵列代数的测深数据拟合方法

张瑞辰，刘雁春，边少锋，李厚朴

海军工程大学导航工程系，湖北 武汉

收稿日期：2017年9月8日；录用日期：2017年9月24日；发布日期：2017年9月30日

摘 要

传统水深测量涉及到大量离散的、分布不规则的水深数据，针对传统拟合方法计算复杂、速度慢的问题，该文将阵列代数思想应用于测深网格拟合。以双二次多项式拟合为例，建立了阵列代数最小二乘解模型，给出了基于正规方程的双二次矩阵方程的最小二乘解法；给出了改进拟合法实现的基本步骤，最后通过实验数据对反距离加权法、全局多项式法及阵列代数最小二乘拟合方法求得未知点水深的效率及精度进行比较。结果表明，在水深测量中引入阵列代数拟合结果效率要高于另外两种方法，可应用于水深测量数据网格化处理中。

关键词 :阵列代数，水深测量，双二次，最小二乘拟合

Copyright © 2017 by authors and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

随着一带一路建设的推进，高效而精确地测量海底地形就显得尤为重要。这就需要将有限的等深线数据，在科学的插值方法下，得到尽可能准确的连续水深数据 [1] 。水深测量数据具有数据大、离散化、网格化、范围广、时空域等特点。传统的水深插值方法如反距离加权插值、全局多项式插值法等，运用这些方法可以较好的内插出水深值，得到较好的海底地形模型。但由于时代的限制，反距离加权插值算法虽简便易行，但计算量大，速度慢；全局多项式插值的适用范围相对有限，易出现病态方程组，对精度影响较大，计算复杂 [2] [3] [4] 。而阵列代数可以有效地将网格状大规模的水深数据块分解为小数据块进行高效率的处理和分析，大大提高了效率；并与最小二乘拟合思想相结合，借助正规方程求出最小二乘解，在一定程度上提高了精度，在水深测量内插计算中具有较好的优势 [5] [6] 。

2. 基于阵列代数的双二次多项式

传统双二次多项式可表示为：

$F\left(x,y\right)=a_{00}+a_{10}x+a_{20}{x}^2+a_{01}y+a_{11}xy+a_{21}{x}^{2}y+a_{02}{y}^2+a_{12}x{y}^2+a_{22}{x}^2{y}^2$ (1)

将(1)通过阵列代数的思想变换为矩阵形式：

$F\left(x,y\right)={\displaystyle \underset{i=0}{\overset{2}{\sum }}{x}^i}{\displaystyle \underset{j=0}{\overset{2}{\sum }}{a}_{ij}{y}^j}=\left(1,x,x^2\right)\left(\begin{array}{ccc}{a}_{00}& a_{01}& a_{02}\\ a_{10}& a_{11}& a_{12}\\ a_{20}& a_{21}& a_{22}\end{array}\right)\left(\begin{array}{c}1\\ y\\ y^2\end{array}\right)$ (2)

线状测深模式是海洋测深的基本作业模式 [2] 。船载测深仪在海面沿计划测线边进行航行测量，获得多条测线上的水深数据，在测区按一定间距加测垂直交于计划测线的检查测线，以便获得交叉点附近的差值数据来检核测深成果。这样，测线和检查线便形成网状，如图1所示。

假设取九条主测线，垂直于主测线选择九条检查测线，如图所示。可以得到 $9\times 9$ 个测点 $\left(x_i,y_i\right)$ 及对应的水深值 $F\left(x_i,y_i\right)$ ， $\left(i=1,\cdots ,9\right)$ 。则得到：

$F_{9\times 9}=A_{9\times 3}{X}_{3\times 3}{B}_{3\times 9}$ (3)

对于形如 $F=AXB$ 形式的矩阵方程。设观测阵列方程为式(4)，即求其最小二乘解最为关键。

$V=AXB-F$ (4)

图1. 抽象到平面后网状测线

3. 基于正规方程的最小二乘解法

本文用 $R_{m\times n}$ 表示所有 $m\times n$ 阶实矩阵全体。 $A\otimes B$ 表示矩阵 $A$ 与 $B$ 的Kronecker乘积，即若 $A=\left(a_{ij}\right)$ ， $B=\left(b_{ij}\right)$ ，则 $A\otimes B=\left(a_{ij}B\right)$ 。 $vecA$ 表示矩阵 $A$ 的按行拉直，即若 $A={\left[a_1^{\text{T}},a_2^{\text{T}},\cdots ,a_m^{\text{T}}\right]}^{\text{T}}$ ，其中 $a_i$ 为 $A$ 的行向量，而 $vecA={\left(a_1{a}_2\cdots a_m\right)}^{\text{T}}$ 。

设 $A\in R_{m\times n},B\in R_{n\times m}$ ，考虑线性矩阵方程：

$F_{m\times m}=A_{m\times n}{X}_{n\times n}{B}_{n\times m}$ (5)

运用Kronecker乘积，方程(5)等价于

$\left[A\otimes B^{\text{T}}\right]\cdot vecX=vecF$ (6)

从而得到方程(6)的正规方程：

$\left[A^{\text{T}}\otimes B\right]\left[A\otimes B^{\text{T}}\right]vecX=A^{\text{T}}\otimes B\cdot vecF$ (7)

由(7)进一步化简得到：

$\left[A^{\text{T}}A\otimes B{B}^{\text{T}}\right]vecX=\left(A^{\text{T}}\otimes B\right)\cdot vecF$ (8)

将(8)还原可得：

$A^{\text{T}}AXB{B}^{\text{T}}=A^{\text{T}}FB$ (9)

由于 $A$ 与 $B$ 均为范德蒙矩阵，所以可知 $B{B}^{\text{T}}$ 与 $A^{\text{T}}A$ 皆为可逆矩阵，于是由方程(9)得到其最小二乘解为：

${\hat{X}}^{\text{T}}={\left(B{B}^{\text{T}}\right)}^{-1}B{F}^{\text{T}}A{\left(A^{\text{T}}A\right)}^{-1}$ (10)

由图2，可以看出，由于 ${\left(B{B}^{\text{T}}\right)}_{3\times 3}$ 为 $3\times 3$ 的矩阵， ${\left(A^{\text{T}}A\right)}_{3\times 3}$ 也为 $3\times 3$ 的矩阵，因此，计算速度优于普通的接线性方程组。

4. 实验对比分析

为验证插值效率和所构造的双二次插值最小二乘内插模型的性能，取双曲抛物面函数：

$f={\left(\frac{x}{0.5}\right)}^2+{\left(\frac{y}{0.6}\right)}^2$ (11)

写成矩阵形式为：

图2. 阵列代数将大矩阵分解为小矩阵块的过程

$f=\left(\begin{array}{ccc}1& x& x^2\end{array}\right)\left(\begin{array}{ccc}0& 0& 4\\ 0& 0& 0\\ 2.77777\cdots & 0& 0\end{array}\right)\left(\begin{array}{c}1\\ y\\ y^2\end{array}\right)$ (12)

取 $x\in \left[-100,100\right]$ 并等分为8份，得到结点 $x_i\left(i=0,1,\cdots ,8\right)$ ；取 $y\in \left[-100,100\right]$ 并等分为8份，得到结点 $y_i\left(i=0,1,\cdots ,8\right)$ 。可得到样本点 $\left(x_i,y_j\right)\left(i,j=0,1,\cdots ,8\right)$ 。采用基于阵列代数的双二次最小二乘拟合模型进行测试，用时0.010秒，拟合结果如下：

$Check-f=\left(\begin{array}{ccc}1& x& x^2\end{array}\right)\left(\begin{array}{ccc}6.4\times {\text{e}}^{-10}& -2.2\times {\text{e}}^{-11}& 4.0\\ -4.4\times {\text{e}}^{-12}& 5.8\times {\text{e}}^{-13}& -4.4\times {\text{e}}^{-15}\\ 2.8& -9.6\times {\text{e}}^{-15}& 7.3\times {\text{e}}^{-17}\end{array}\right)\left(\begin{array}{c}1\\ y\\ y^2\end{array}\right)$ (13)

取 $x\in \left[-80,80\right]$ 并等分为49份，得到结点 $x_i\left(i=0,1,\cdots ,49\right)$ ；取 $y\in \left[-100,100\right]$ 并等分为49份，得到结点 $y_i\left(i=0,1,\cdots ,49\right)$ 。可得到检查点 $\left(x_i,y_j\right)\left(i,j=0,1,\cdots ，49\right)$ ，求其拟合精度的公式为：

$\sigma =\sqrt{\frac{{\delta }_1^2+{\delta }_2^2+\cdots +{\delta }_n^2}{n}}=\sqrt{\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\delta }_i^2}}{n}}=\sqrt{\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(f_i-check\_f_i\right)}^2}}{n}}\left(n=2500\right)$ (14)

其中， $\sigma$ 为标准差， $f_i$ 为真值， $check\_f_i$ 为拟合值， ${\delta }_i$ 为真值 $f_i$ 与拟合值 $check\_f_i$ 的差值， $n$ 为试样总数 $\left(n=2500\right)$ ， $i$ 为各次拟合值。

由式(14)可得精度结果：

$\sigma =3.3659\times {10}^{-17}$ (15)

如图3所示，图中曲面为基于阵列代数的双二次插值最小二乘拟合模型及原函数图像的叠加结果图，可以误差很小，肉眼不可见。

用Matlab编程，使用反距离加权插值法(IDW)，同样的采样点，当 $p=4$ 时得到内插标准差较小，标准差为：

$\sigma =1.49\times {10}^{-2}$ (16)

同时，将内插出的点与原函数点绘制成图，进行比较。可以看出，当权重为4时，内插出的图像有很多小的凹陷，不平整，如图3(b)所示。用Matlab编程得到的反距离加权插值法的插值效果一般，该方法总共需要花费时长为0.209 s，远远超出基于阵列代数的双二次插值最小二乘内插模型的时间。

采用传统不完全双二次函数进行拟合。不完全双二次函数如下式(17)所示。

$f\left(x,y\right)=a_0+a_{1}x+a_{2}y+a_{3}xy+a_4{x}^2+a_5{y}^2$ (17)

图3. (a) 本文方法图像与原函数图像对比图(b) IDW( $p=4$ )与原函数图像对比图

表1. 三种方法用时与精度比较

用Matlab编程用时为0.072秒，得到的拟合结果为：

$\begin{array}{l}f\left(x,y\right)=a_0+a_{1}x+a_{2}y+a_{3}xy+a_4{x}^2+a_5{y}^2\\ a_0=a_1=a_2=a_3=0,a_4=4,a_5=2.7778\end{array}$ (18)

即：

$f\left(x,y\right)=4x^2+2.7778y^2$ (19)

标准差为：

$\sigma =7.105\times {10}^{-17}$ (20)

比较三种方法的用时和精度，可以得到表1。

由表1可以看出，相比较传统双二次多项式最小二乘拟合，基于阵列代数的双二次多项式最小二乘拟合算法能够有效地分解大型稀疏矩阵，不需要复杂的奇异值分解等分解形式，即可求出拟合值，同时具有较好的精度和效率。相比于常用水深内插方法反距离加权插值法，在用时和精度方面，本论文方法都有较大优势。

5. 结束语

双二次插值最小二乘内插模型属于整体内插方法，容易受区域内极值的影响，适用于较平缓地区；同时，阵列代数的应用要求所给的数据必须是格网数据。基于阵列代数的双二次插值最小二乘内插模型，将水深测量的大数据块分割成若干小数据块，同时将双二次多项式与最小二乘拟合相结合得出基于正规方程的最小二乘解，一定程度上提高了计算精度，同时有效得提高了计算机的运算效率，减小了计算机存储量，具有非常重要的实用价值。

基金项目

国家自然科学基金项目(41576105, 41631072)。

文章引用

张瑞辰,刘雁春,边少锋,李厚朴. 基于阵列代数的测深数据拟合方法
Fitting Method of Hydrographic Data Based on Array Algebra[J]. 应用数学进展, 2017, 06(06): 795-800. http://dx.doi.org/10.12677/AAM.2017.66096

参考文献 (References)