地形图的绘制

1、 海底地形图的绘制模型华东师范大学江* 常* 黄* 海底地形图的绘制模型摘要：该“海底地形图的绘制模型”首先通过各个角度的假设，建立一个较为理想的海洋测绘环境；考虑该地形图能够表示为在一个三维笛卡尔坐标内的函数，而“插值双三次样条函数”能够较好地对图形进行模拟，我们对函数的形式做出一些改进，以满足海底地形图的实际要求；为了使绘制的图形更贴近与现实，我们将待测的海洋区域以网格方式划分为若干小区域，通过对每个小区域的地形图的绘制，并将所有的小区域内的图形连接，最后生成整个待测海洋区域的海底地形图。一 问题重述海洋测绘船利用声纳绘制海底的地形图。测绘船上的声纳向海底发射声脉冲，随后接收从海底反射的脉

2、冲。发射的范围为与指向海底的铅垂线夹角从2°30°之间。船只以2米/秒的速度行进，声脉冲在海水中传播的速度约为1500米/秒。试建立绘制海底地形图的数学模型，并对绘制海底地形图的方法提出具体建议。二 模型假设1 造成海面的不平坦因素有很多，比如潮汐，风浪等，假设这些因素的影响在测量过程中间可以忽略。2 讨论区域的海底曲面是光滑的；更确切些，可以认为曲面的一阶、二阶导数是连续的。因为我们可以认为讨论区域为浅水海域，由于长期的海水水流作用，形成的是以砾石或沙为主要组成部分的海底，不存在珊瑚礁、水底峡谷、山脊等不可预料的突变地形。为了方便分析误差，假设海底临近的两点之间的与水平面

3、的夹角在15°±15°范围内，且反射回接收装置的声波线与铅垂线之间的夹角在16°±14°范围内。3 声波在传播途中受海水介质不均匀分布和海面、海底的影响和制约，会产生折射、散射、反射和干涉，会产生声线弯曲、信号起伏和畸变，造成传播途径的改变，以及出现声阴区，严重影响声纳的作用距离和测量精度1。假设海水介质非常均匀，造成声波在海水中传播方向不稳定的因素可以忽略。4 由于海水的浓度，温度，密度，压力等等因素起伏变化不定，声音在其中的传播速度肯定也在起伏变化，由于这种变化没有固定的规律可循，假设声波在海水中以一个稳定的速度传播。5 假设回波

4、接受装置接受到的回波是可辨别的，即能够确认是发射装置刚刚发射的声波的回波，不会出现辨别不出回波信号的问题。6 假设声波在到达海底前引起的回波都是不明显的，即假设收到的回波信号是遇到海底地面返回的。7 假设用于测量的声纳设备本身是精准的，至少在测量过程中不会因为声纳设备本身的原因而导致采集的数据出错。8 由于测量条件的限制 ，假设测量误差在小于10%的范围内都是允许的。 三 问题分析相关技术背景：1 声纳传播方式题目中要求建立利用声纳绘制海底地形图的数学模型，所以首先有必要了解一下声纳的工作原理。声波是目前已知的唯一能在海水中远程传播的波，声纳（SONAR）是利用声波在水下的传播特性，通过电声转

5、换装置和信号处理，完成对水下物体进行探测和定位识别的方法及所用设备的总称2。当我们通过声纳向水中发射声波脉冲时，声波能够以一定的速度在水中沿直线传播，并能从水底反射回波。为了得到某个目标位置的深度值，我们可以根据声波的这一特性，通过计算发射声波脉冲和回波到达的时间差来求得目标位置的深度。这是利用声纳测量海底地形的基本原理。为了更好的理解利用声纳测量海底地形的原理和过程，我们需要了解一下声波的传播过程。波在均匀各向同性介质中传播时，波面及波前的形状不变，波线（波的传播方向成为波线）也保持为直线，沿途不会改变波的传播方向。但是当遇到障碍物或是从一种介质传播到另一种介质中时，波面的形状和波线方向将发

6、生改变。通过惠更斯原理可以知道，在介质中，波传到的各点不论在同一波前或不同波前上，都可看作是发射子波的波源3。因此，测绘船发射的声脉冲是以一个球形波阵面的形式向海底传播的，而且在遇到海底障碍返回时，又是以球形波阵面的形式返回。2 曲面造型方法海底曲面可由插值双三次样条函数去模拟，前提是测出一些已知点的深度。关于插值双三次样条函数说明如下：定义：设在uw平面上的矩形区域R:a,b×c,d上给定一个矩形网络分割=u×w,其中u： a=u< u< u< <u=b,w: c=w< w< w< <w=d.凡在R上满足下述两条件的一个函数

7、x(u,w)称为双三次样条函数。在每个子矩形Rij：U，U ×W，W（i=1,2n;j=1,2,m）上，x(u,w)关于u和w都是三次多项式函数，即 (1.1) 在整个R上，函数x(u,w)的偏导数,(=0,1,2)是连续的，或简记为x(u,w)。 此外，在给定一数组（i=0,1,2n; j=0,1,2,m）之后，假如x(u,w)在节点处再满足插值条件，即： x(u,w)= （i=0,1,2n; j=0,1,2,m）(1.2)则称x(u,w)为插值双三次样条函数。给出以下记号：(1) 矩形R的四条边界上的所有节点处的一阶法向偏导数 (1.4)=(,) （j=0,1,2,m ; =0,

8、n）. =(,) （i=0,1,2,n; =0,m）(2) 矩形R的四个顶点处的二阶混合偏导数 （=0,n; =0,m）.(1.5)由存在唯一性定理可知：给定uw平面上的矩形R、它的一种矩形分割和(m+3)(n+3)个常数, , (i=0,1,n; j=0,1,m; =0,n; =0,m)之后，恰好存在一个插值双三次样条函数x(u,w),使它是以(1.2)为插值条件，以(1.4)和(1.5)为边界条件的。当我们应用上述定理到单个子矩形Rij时，便可知道：只要给定了Rij的四个定点上的函数值，两个方向的一阶偏导数以及二阶混合偏导数等16个数据，也就是给定了四阶方阵 C = 那么，Rij上的双三次

9、函数x(u,w)就被唯一地决定了。矩阵Cij中各元素的意义分别是=x(u); = x(), = x(), =.且易得出Rij上的双三次函数一定可唯一地写成x(u,w)=(u-u)A(h)C A(g)(w-w), (u,w)Rij, ( 1.7 )式中记号(u-u)代表关于变元(u-u)的1×4阶矩阵：(u-u)=(u-u)(u-u)(u-u),1,A(h)代表4×4阶矩阵：h=（i=1，2，.,n）g=(j=1, 2, ,m)符号T代表矩阵的转置。这样一来，双三次样条函数x(u,w)在子矩形Rij上的表示式（1.7）完全决定于矩阵Cij。Cij的左上角四个元素由插值条件（1

10、.2）直接给出，其余诸p,q,S共12个元素则需要按照边界条件(1.4)和(1.5)，并求解连续性方程获得。首先，固定w=w,则函数x(u, w)S（u，u）.由三次样条函数的m次连续性方程，我们有 (i=1,2,n-1), (1.9)式中= , =.在边界条件 (1.4) 式中，(1.9) 的边界切向 是已知的，由此解出 p (i=0,1,n) 。然后取j=0,1,m，我们便求得R上所有节点（u，w）处沿u方向的一阶偏导数(i=0,1,n;j=0,1,m). 对称地，当固定u= u时，则函数x(u,w) S（u，w）.于是得出相应的m连续性方程 (j=1,2, m-1), (1.10)式中记

11、= ，= 。其边界切向,和也是从(1.4)获得的，由此我们解出(j=0,1,m).这样，R上所有节点（u，w）处沿w方向的一阶偏导数（i=0,1,m）就全部确定了。为了计算所有节点处的二阶混合偏导数Sij，我们首先固定u= u,则函数 （u，w） S（w，w）。它在各节点w处的函数值（u，w）= 已从(1.9)求得。因此我们得到关于变元w的三次样条函数的m连续性方程 (j=1,2,m-1) (1.11)假如还给定了连续性方程（1.11）的边界切向Si0和Sim的话，那末就解得 (j=0,1,m).这样R上所有节点（u，w）处的二阶混合偏导数 (i=0,1,n;j=0,1,m)也就全部确定了。其

12、实，（1.11）的边界切向S和S是这样求得的：沿R的两条边界w=w与w=w,函数 (u, w)与(u, w) S（u;u）.它们在各节点ui处的函数值 (u, w)=q和( u, w)=q可从（1.10）获得。于是我们有关于变元u的m连续性方程 (i=1,2,n-1), (1.12)它们的边界条件由（1.5）式提供。由此解出S和S (i=0,1,n).至此，矩阵(1.6)中所有各元均已求得，因此R上的双三次样条函数x(u,m)也就被确定了，它们在各子矩形上的表达式就是（1.7）。余下的问题是：参数平面(u,w)上的定义域R及其分割应该如何确定呢？通常的方法是，适当固定某行j,而取u方向的节点间

13、距h,(i=1,2,.,n).再适当固定某列i，取w方向的节点间距g(j=1,2,.,m).若令定义域R=0,H×0,G,这里的H= G=.这样，我们就确定了定义域R及其分割。这种取法实际上是构造了两个方向上的累加弦长参数分割。累加弦长参数三次样条曲线具有良好的光顺性，因此现在构造的样条曲面的网络曲线将是光顺的。运用这种方法的前提，就是型值点x(u,w)之间的间隔比较均匀。倘若一侧间距很大，而另一侧间距很小，象扇形分布那样，我们将不会收到好的效果。前所讲的所谓适当固定某行j，是指这一行上各点x(i, j) 的间隔大致是一个平均值。因此我们作以下假设：讨论区域的海底曲面是光滑的；更确切

14、些，可以认为曲面的一阶、二阶导数是连续的。因为我们可以认为讨论区域为浅水海域，由于长期的海水水流作用，形成的是以砾石或沙为主要组成部分的海底，不存在珊瑚礁、水底峡谷、山脊等不可预料的突变地形。至此，所有问题都解决了；由测绘船测出规则点的深度值，则利用以上插值双三次样条函数的求解（规定边界上的一阶偏导数和二阶混合导数为0）可得出海底曲面；如果测的是随机点的深度值，则可由它来估算规则点的深度值，再求解；方法如下： 由随机点的深度值推测出分割节点上的水深值，我们采用在地球科学用得较多的IDW（inverse distance weighting）方法。IDW 方法是地球科学中常用的计算机插值技术，最

15、早可见于60年代的地质轮廓描述和矿产推测的有关文献中。IDW的基本思想是估计已知数据点对推测点的影响，假定推测点的值为周围已知点的值的线性组合。有前假定，每个已知点对推测点的影响反比于它们的距离的某种方式，这儿我们用距离的p次方。设已知点()的值为,则 X=,其中，为点(u,w)和点（）之间的距离，这儿K起着规范化的作用。对于IDW推测出的未知点的值，显然在已知点的值的范围之内，因此，只有靠近已知点的附近才会出现极值情况。上式中的p反映了未知点对已知点的偏倚情况，p越大，表示未知点的值受近距离点的值的影响越大（相对于远距离点）。例如，我们取p=3,依照上式就可以根据已知节点深度值计算出所有未知

16、格点的水深值（已知随机点的深度值用船去测），然后用前面介绍的双三次样条插值方法（规定边界上的一阶偏导数和二阶混合导数为0）就可求出海底拟合曲面。四 符号定义 声波在海水中的传播速度 测绘船的航速 声波从发射到接收所经过的时间 海底某一点（x,y）处的水深 在Tswave时间内测绘船走过的位移 在Tswave时间内声波走过的路程 反射声波走过的路程 海底临近的两点之间的倾斜角度 反射声波线与铅垂线之间的夹角 由于海底地形起伏引起的深度测量误差 由于轮船航行速度引起的深度测量误差 被测海域的平均深度 划分后得到的小区域的边长 一次声纳发射可以覆盖的区域的直径长度五模型建立概述：在现有的假设和技术背

17、景下，我们现建立如下测绘海底地形的方案：首 先将被测绘的海底区域以网格方式划分为若干小区域，然后分别测量网格上每个交叉点处的海底深度；通过双三次样条函数的计算方法，得出每个小区域的曲面图形；最后通过曲面整合的方式，得到整个海底地形图。误差分析：1 海底地形引起的误差由于海底地形是呈波动起伏状的，所以向海底传播的声脉冲的球形波阵面最前沿不会同时到达海底并反射，而是这个波阵面上的一点或几点首先遇到海底地面并反射回波，这个过程如图1.1所示。 图1.1从图中可以看出测绘船回波接受装置接受到的回波信号将不是位于轮船铅垂线正下方的A点反射回的回波信号，而是位于声脉冲覆盖范围内的另外一点B点的回波。而且B

18、点是声脉冲发射出来以后波阵面到达的第一个反射点，也就是整个声脉冲覆盖区域内与发射装置直线距离最近的海底表面的一点。这样测量的结果自然会与理论值有一个偏差，这是由于海底地形起伏引起的深度测量误差，我们可以利用全微分计算公式估算一下这个误差造成的影响是否在我们设定的误差范围之内。 图1.2根据图1.2所示，可以得到 (2.1)假设，是实际测量时和的偏差，则全增量的绝对值就是利用上述公式所产生的误差由于和都很小，我们可以用来近似地代替，这样就得到的误差为 (2.2)将=16°，=15°，= 14°，=15°代入式2.2，可以得到的绝对误差约为 由此可以得到的相

19、对误差约为 显然这一误差在我们的允许范围内，所以在实际测量过程中我们可以用的值来代替的值。 2 船速引起的误差接着我们需要考虑到船速对测量结果的影响。 图1.3如上图1.3所示，船在向正下方的A点发出声脉冲时处在测绘方向上的B点，此后经过时间船有了一段位移 (2.3)在C点，船收到了回波信号。声波所走过的路程 (2.4)所以实际测到的A点深度值为 (2.5) 由三角关系可得 (2.6) 所以由于船运动所引起的绝对误差为 （2.7） 由于未知，此绝对误差不好估算，我们可以通过估算相应的相对误差来看其影响。 (2.8) 由于= 2m/s, = 1500m/s，而海洋的平均深度在几十米到几千米的范围

20、内4，因此测量时声波在海水中传播到反射回来所花的时间很短，一般为几秒钟，这段时间里穿走过的路程很短，因此可以认为 (2.9) 将(2.9)代入(2.8)可得 = 这个误差远远小于所允许的误差范围，所以我们基本可以忽略船速对测量结果的影响。海底区域划分 最后我们需要考虑如何对一片海域划分，以得到所需的数据。很明显，划分时既不能将区域划分的太大，这样得到的数据太少，不足以反映出实际的海底地形特征；同时划分也不能太小，这样虽然会更真实地反映出实际的海底地形特征，但是会造成数据量过多，影响计算的效率。从这种考虑角度出发，我们认为划分的前提是利用较少的测量次数尽可能的覆盖整个区域，从而可以从一定程度上确

21、保采集到的数据能够比较完整的反映整个测量区域的海底地形分布。 考虑要测量的区域是一片如图1.3所示的正方形区域 ABCD，我们现在讨论如何将该区域进行划分，即等分的划分为(m+1)(n+1)个小正方形。 图1.4 如果能够确定每个划分后的小区域的边长，则划分就可以得到。如图1.5所示，若划分的小区域的大小满足如下条件：当测绘船分别在该小区域的对角线上的两个网格交叉点进行测量时，两处的声纳测量范围刚好相互外切。由于划分的小区域都处于声纳的测量范围，因此得到的数据能够较真实的反映该区域的海底地形，而且又不会产生划分过细的情况。 图1.5 由于声脉冲发射的范围为与指向海底的铅垂线夹角从2°

22、-30°之间，因此可以借助被测海域的平均深度，估算出一次声纳发射可以覆盖的范围，以此来确定。 (2.10) (2.11) 因此在划分小区域之前，还要进行一项工作，即估测北侧海域的平均深度，这可以通过使用上述方法随机测量几组数据来估测。 模型测试：模型已基本建立，现对模型进行模拟测试，以确定模型的可行性。A 假设：假设测绘的海域大小为4km * 4km = 16平方公里假设该海域平均深度为1200m（以南海的平均深度为例）按照已分析的海域划分方式（根据公式1.11），该海域应被划分为约1km *1km的小区域，共4 * 4 = 16块B 计算根据假设的海域平均深度1200m，在1200

23、 1200 * 30% 以内，随机产生5 * 5个交叉点数据：X= j i012340912.61020.61430.01148.5873.211318.7985.31315.01532.91472.02849.91138.11471.51458.6991.731312.5883.4960.31165.41268.441266.31527.41171.31124.1875.0根据公式（1.9）（1.10）（1.11）（1.12），所计算出的矩阵p q s为：p = j i0123400000010137.110793.2732310.0179020-195.9429-248.5929-309.

24、7714030340.9607-163.0018-173.43210400000q = j i0123400000010-111.9054436.52148.6196020405.1179244.3286-420.9321030-326.9411251.1643168.28390400000S = j i012340000001032.9034148.2061-107.0064020-162.0528-97.5673210.8365030-12.0243-440.9082232.47840400000以海域最左上角的小区域为例，生成的函数为z(u,w) = (-2451598904247826

25、143/1231453023109120*(u-1)3+3444389734408632031/1231453023109120*(u-1)2-6183/5)*(w-1)3+(629309/280*(u-1)3-886377/280*(u-1)2+6183/5)*(w-1)2+(-171867/140*(u-1)3+47721/28*(u-1)2-5643/5)*(w-1)-57115/56*(u-1)3+409167/280*(u-1)2-108该小区域的海底地形图为：地形图尚未模拟出来通过将4 * 4个小区域的海底地形图的连接，生成的整个待测海域的地形图为：地形图尚未模拟出来六模型分析及改

26、进双三次样条拟合是很有用的一种方法，它具有局部性和连续光滑性；但是IDW算出的结点水深处于已知数据点的凸包之中，它还不是“保山脊”的，即相似已知点与不相似已知点具有相同的影响。另一缺点是模型无法给出误差的解析估计。直接用差分法国用IDW方法等求出足够精细分割的节点上的水深值，采用计算机绘图的“以值代曲”原理用直线联结相邻叠点来描述曲面，计算量会比较大，但能更多地反映已知数据点的影响。参考文献：1 声纳技术 2 声纳技术 3 惠更斯原理 /jp2005/25/dzjc/p02/ch07/sec06/index.htm4 挑战地球深渊海沟5 苏步青，刘鼎元 计算

27、几何 上海科学技术出版社 1981附录：测试用例中各个子矩形上的曲面表达式：(坐标轴为u和w) (-2451598904247826143/1231453023109120*(u-1)3+3444389734408632031/1231453023109120*(u-1)2-6183/5)*(w-1)3+(629309/280*(u-1)3-886377/280*(u-1)2+6183/5)*(w-1)2+(-171867/140*(u-1)3+47721/28*(u-1)2-5643/5)*(w-1)-57115/56*(u-1)3+409167/280*(u-1)2-108, (-5014

28、601/3920*(u-1)3+7801413/3920*(u-1)2-11244/5)*(w-2)3+(310463/560*(u-1)3-464127/560*(u-1)2+11244/5)*(w-2)2+(-32039174199553129/175921860444160*(u-1)3+326231932640540383/1231453023109120*(u-1)2-9197/5)*(w-2)-1849/16*(u-1)3+15621/80*(u-1)2-2047/5, (3740951597090580629/1970324836974592*(u-1)3-64493214151

29、865571299/19703248369745920*(u-1)2-1171/2)*(w-3)3+(-56641/28*(u-1)3+476853/140*(u-1)2+1171/2)*(w-3)2+(2022723/3920*(u-1)3-755211/784*(u-1)2-867)*(w-3)+437951/560*(u-1)3-143523/112*(u-1)2+563/2, (374533/140*(u-1)3-518397/140*(u-1)2-1613/5)*(w-4)3+(-374533/140*(u-1)3+518397/140*(u-1)2+1613/5)*(w-4)2+(

30、17852838784066013929/9851624184872960*(u-1)3-51558524178336618467/19703248369745920*(u-1)2-5979/10)*(w-4)+206869/280*(u-1)3-266901/280*(u-1)2+2753/10 (2223142849447798043/1231453023109120*(u-2)3-1419772112766089031/615726511554560*(u-2)2+40201050963649683/123145302310912-466017243926214367/123145302

31、3109120*u)*(w-1)3+(-831179/280*(u-2)3+553443/140*(u-2)2-70583/140+115173/280*u)*(w-1)2+(46661/28*(u-2)3-30333/14*(u-2)2-3621/35-38391/140*u)*(w-1)+364569/280*(u-2)3-250113/140*(u-2)2+85067/140-38391/280*u, (4355587/3920*(u-2)3-3155247/1960*(u-2)2-714653/392+559023/3920*u)*(w-2)3+(-378613/560*(u-2)3+

32、233631/280*(u-2)2+549697/280+627/112*u)*(w-2)2+(-727774835689052901/1231453023109120*(u-2)3+593772327273121977/615726511554560*(u-2)2-1061229695822088481/615726511554560-20358792909534943/1231453023109120*u)*(w-2)+8311/80*(u-2)3-6057/40*(u-2)2-3339/8+3507/80*u, (-828019085719304573/1970324836974592*

33、(u-2)3+27133476527957160451/19703248369745920*(u-2)2-5104189319418553713/19703248369745920-1047367524438357733/1231453023109120*u)*(w-3)3+(38729/28*(u-2)3-186381/70*(u-2)2+16859/35+104091/140*u)*(w-3)2+(-1297473/3920*(u-2)3+322401/392*(u-2)2-218409/392-1483941/3920*u)*(w-3)-345613/560*(u-2)3+298119/

34、280*(u-2)2+12073/56-121377/560*u, (-72355/28*(u-2)3+113973/35*(u-2)2-181319/70+17361/28*u)*(w-4)3+(72355/28*(u-2)3-113973/35*(u-2)2+181319/70-17361/28*u)*(w-4)2+(-7974101270744457073/9851624184872960*(u-2)3+20947641702275428867/19703248369745920*(u-2)2-35633387505921369457/19703248369745920+49999804

35、346535583/246290604621824*u)*(w-4)-227291/280*(u-2)3+127083/140*(u-2)2-78279/140+17361/56*u (-24861685079709507961/9851624184872960*(u-3)3+859928044083609667/281474976710656*(u-3)2-796352789571249301/307863255777280+262161426676411819/615726511554560*u)*(w-1)3+(1067823/280*(u-3)3-39321/8*(u-3)2+2434

36、59/70-20574/35*u)*(w-1)2+(-62553/28*(u-3)3+56283/20*(u-3)2-182137/70+13716/35*u)*(w-1)-442293/280*(u-3)3+84039/40*(u-3)2-30661/35+6858/35*u, (-4216851/3920*(u-3)3+29607/16*(u-3)2-2609689/980+251199/980*u)*(w-2)3+(58641/112*(u-3)3-95511/80*(u-3)2+447997/140-9909/28*u)*(w-2)2+(7121225808387579527/9851

37、624184872960*(u-3)3-1163050205721919153/1407374883553280*(u-3)2-559495994715743381/307863255777280+85703004557897131/615726511554560*u)*(w-2)+3489/80*(u-3)3+12819/80*(u-3)2-9827/20+1053/20*u, (217849144663418143/1231453023109120*(u-3)3-22994498988539915/61572651155456*(u-3)2-4658127818853749383/1231

38、453023109120+791781255832591249/1231453023109120*u)*(w-3)3+(-127683/140*(u-3)3+29739/20*(u-3)2+191043/70-30249/70*u)*(w-3)2+(-579119/3920*(u-3)3-3489/80*(u-3)2-1981089/980+53583/196*u)*(w-3)+284261/560*(u-3)3-62943/80*(u-3)2-23889/140+1713/28*u, (261323/140*(u-3)3-53271/20*(u-3)2+126051/70-21684/35*

39、u)*(w-4)3+(-261323/140*(u-3)3+53271/20*(u-3)2-126051/70+21684/35*u)*(w-4)2+(116958271844871635/246290604621824*(u-3)3-38995691707564043/61572651155456*(u-3)2-799123480336700039/1231453023109120-121833742012064879/1231453023109120*u)*(w-4)+183847/280*(u-3)3-36669/40*(u-3)2+19547/14-10842/35*u (23249945540632397879/9851624184872960*(u-4)3-14888581823860983887/4925812092436480*(u-4)2+9328218952510060381/2462906046218