海底地形测量图的插值模型_海底地形

海底地形测量图的插值模型

摘要：

随着全球经济一体化和信息技术的发展，企业之间的合作日益加强，跨地区甚至跨国合作制造的趋势日益明显。国际上越来越多的制造企业不断地将大量常规业务 “外包”（outsourcing）出去给发展中国家，而只保留最为核心的业务（如市场、关键系统设计和系统集成、总装配以及销售）。在这些合作生产的过程中，大量的物资和信息在更为广阔的地域间转移、储存和交换，国际物流活动将日益频繁，港口作为国际物流活动主要的载体，在国际贸易与国际经济合作中愈来愈发挥着极其重要的作用。但是 海底的地形是十分复杂的，它不仅分布有巍峨的海底山脉、平缓的海底平原，而且还有许多陡峭的海底深沟。为使轮船进出港口安全，就需要了解港口航道的海底地形。所以我们就要对海底地形有足够的了解，预测出哪些区域是船只的禁入点，避免船只在危险区域搁浅。为了研究此问题，我们利用低潮时测得14组数据，并进行了6个基本假设，在此基础上我们便利用插入与拟合的思想来利用光滑曲面来模拟海底曲面。但考虑到保凸性及光滑性的要求，我们采用双三次样条函数来模拟海底曲面。首先利用测量的14个随机点的深度，以随机点的坐标将待测矩形区域划分成1414的网格。然后我们通过某种加权平均来逼近未知网格接点上的深度，采用距离的倒数作权重反映出距离越小影响越大。求出所有14×14个网格点上的深度后，调用IMST中的双三次样条子程序，通过插值得到海底曲面，然后再加细网格，划分成50×50的网格，计算这加细网格接点上的深度;最后找出两个危险区域分别在深度为4Ft的两个点的周围，并借助于Matlab中的绘图程序，绘出海底的二维、三维网格图及等高线图。以不同的颜色将水深小于船只的吃水深度标示出来，作为船只的避免进入区域，并作出水深小于船只吃水深度的海域范围,并绘出等高线。在找几个近似大小的海底拓扑地图，每个随即取14个数据点，把我们的模型应用上去，通过比较实际地图及由模型画出地图比较接近。这说明该模型具有十分有效的实际作用。该模型充分利用了已知点的数据，给出了求未知网格点上深度的近似方法，用保凸性较好的双三次样条拟合了海底曲面，得到了较为满意的结 果。但在实际计算中，三次样条可能会导致数值上的不稳定，遇到这种情况，可以用加细网格点的办法来解决，也可以用稳定性较好的B—样条来拟合。

一、问题的提出

随着全球经济一体化和信息技术的发展，企业之间的合作日益加强，跨地区甚至跨国合作制造的趋势日益明显。国际上越来越多的制造企业不断地将大量常规业务 “外包”（outsourcing）出去给发展中国家，而只保留最为核心的业务（如市场、关键系统设计和系统集成、总装配以及销售）。在这些合作生产的过程中，大量的物资和信息在更为广阔的地域间转移、储存和交换，国际物流活动将日益频繁，港口作为国际物流活动主要的载体，在国际贸易与国际经济合作中愈来愈发挥着极其重要的作用。

但是 海底的地形是十分复杂的，它不仅分布有巍峨的海底山脉、平缓的海底平原，而且还有许多陡峭的海底深沟。为使轮船进出港口安全，就需要了解港口航道的海底地形。

表1给出了水面直坐标（x，y）处的水深z，这是在低潮测得的。如果船的吃水深度为5Ft，试问在矩形区域(75，200)×(50，50)中行船应避免进入哪些区域？其中水平方向的坐标x, y以Yd(=0.914m)为单位，水深方向以Ft(=30.48cm)为单位。

二、模型假设

1. 海底是光滑的，且无暗礁;

2. 每个给定的数据都影响着其它未知点的深度，且距离越近，影响越大; 3. 任何两个数据点之间深度的变化都影响着其它未知点的深度;

4. 对于给定两个数据的变化对于某一未知点深度的影响，取决于三个距离：两个数据点连线与该未知点的垂直距离、该未知点与离它最近的那个数据点之间的距离及两个给定数据点之间的距离;

5. 两个数据点之间深度的变化对某一未知点深度的影响与这两个数据点连线线性传播; 6. 每一个给定的数据点对某一未知点深度的影响它们之间的距离平方成反比。 三、分析与建摸

根据假设条件海底是光滑的，无暗礁，因此很自然地想到利用光滑曲面来拟合海底曲面。例如可以用二维拉格朗日(Lagrange)插值或双三次样条函数来逼近。考虑到保凸性及光滑性的要求，我们采用双三次样条函数来拟合.

为了用双三次样条函数插值，必须知道xoy平面内所

有网格上的深度，而所给数据的14个随机点并不构成任何网格。所以 第一步 生成网格：

我们采取最简单的方法是过个14 数据点分别作平行于x，y 轴的直线，划分成不规则1414的网格。

第二步 确定那些未知点数据网格上的深度： 应该说所有数据对未知网格的深度都有影响，只是越靠近的数据点影响越大。由于我们对海底所知甚少，所以只能通过某种加权平均来逼近未知网格接点上的深度，采用距离的倒数作权重反映出距离越小影响越大。但仅用加权平均来逼近未知点的深度，它不能反映数据点深度的变化趋势 。

设Q1，Q2点的深度分别为4Ft和8Ft，G是Q1Q2连线上未知深度的点(见图1)，

Q1Q2Q2G4Ft，求G点的深度zg。

根据光滑假设，由点Q2经点Q1到G点的深度应渐渐变浅，因此，未知点G的深度应小于4Ft。

4

Q2

4

G

下面用三种外推公式加以分析：

Z(Q1,Q1,G)Zg

GQ1

2



Z(Q2,Q2,G)GQ21GQ2

2

2

1GQ1

2

4.8Ft

其中z(Q1, Q1, G)和z(Q2, Q2, G)分别表示Q1, Q2点的深度4Ft和8Ft。显然，zg不符合小于4Ft的期望。

4

Q2

4

G

2. 线性外推公式

Zg0Ft

这个外推值又太小，也不符合实际情况。

4

Q2

4

G

4

Q2

4

G

2. 线性外推公式

Zg0Ft

这个外推值又太小，也不符合实际情况。

4

Q2

4

G

2. 线性外推公式

Z(Q,Q,G)Z(Q,Q,G)Z(Q,Q,G)



GQGQGQQQ3.4Ft

ZgGQGQGQQQ12112

其中z(Q1, Q2, G)为线性外推值。这种外推比较合理。

4

Q4

G

现将上述一维情况变成二维情况，即未知点G不在已知点Q1，Q2的连线上。根据假设4，Q1，Q2点对G点深度的影响取决于三个距离：GP,GQ,GQ，其中P是G到Q1Q2 延长线的垂足。利用上面的分析结果修改权因子，得到如下的加权平均：

Z(Q,Q,G)Z(Q,Q,G)Z(Q,Q,G)

2GQGQGPGQQQZ

g111



2GQGQGPGQQQ12112

P

Q1

QQ2

Q2

GP

1

G

zg

为了考虑所有给定点的影响，将上述加权平均推广到所有点对。设G是某一未知深度的网格点，Qi，Qj 是已知深度的网格点，记

GQmin{GQ,GQ}，Pij是G到

ijij

Qi，Qj 连线的垂足。采用下面的加权平均来逼近G点的深度。

Z(Q,Q,P)NNijij



i1j1GPGQQQ

ijijijZ

gNN1



i1j1GPGQQQ

ijijij

这里z(Qi, Qj, Pij)是由Qi，Qj点的深度线性外推的Pij点的深度。

四、模型求解

求出所有14×14个网格点上的深度后，调用IMST中的双三次样条子程序，通过插值得到海底曲面，然后再加细网格，划分成50×50的网格，计算这加细网格接点上的深度;最后找出两个危险区域分别在深度为4Ft的两个点的周围，并借助于Matlab中的绘图程序，绘出海底的二维

三维网格图及等高线图。 五、模型检验

检测模型比较简单，只要找几个近似大小的海底拓扑地图，每个随即取14个数据点，把我们的模型应用上去，

通过比较实际地图及由模型画出地图即可。 六、模型的评价

该模型充分利用了已知点的数据，给出了求未知网格点上深度的近似方法，用保凸性较好的双三次样条拟合了海底曲面，得到了较为满意的结果。但在实际计算中，三次样条可能会导致数值上的不稳定，遇到这种情况，可以用加细网格点的办法来解决，也可以用稳定性较好的B—样条来拟合。

七MATLAB实现代码

%数据输入

x_Yd=[129.0 140.0 108.5 88.0 185.5 195.5 105.5 157.5 107.5 77.0 81.0 162.0 162.0 117.5]; y_Yd=[ 7.5 141.5 28.0 147.0 22.5 137.5 85.5 -6.5 -81.0 3.0 56.5 -66.5 84.0 -38.5];

z_Ft=[ 4 8 6 8 6 8 8 9 9 8 8 9 4 9];

%由已知14点计算出14*14点（其中有两点重复）

%计算QiQj的值Qi_Qj Qi_Qj=zeros(14,14); for i=1:14

for j=1:14

Qi_Qj(i,j)=sqrt((x_Yd(i)-x_Yd(j)).^2+(y_Yd(i)-y_Yd(j)).^2); end end

[G_x,G_y]=meshgrid(x_Yd,y_Yd); %计算所求196个G点的坐标

%求出每个G点对应的深度Z_g for m=1:14 for k=1:14

if ((m==k)|((m==12)&&(k==13))|((m==13)&&(k==12))) %排除Qi,Qj为同一点和重复点的情况

Z_g(m,k)=z_Ft(m); continue; else

%计算GQij=min{GQi,GQj} G_Qi=zeros(1,14); G_Qj=zeros(1,14); for i=1:14

G_Qi(i)=sqrt((G_y(m,k)-y_Yd(i)).^2+(G_x(m,k)-x_Yd(i)).^2); end

for j=1:14

G_Qj(j)=sqrt((G_x(m,k)-x_Yd(j)).^2+(G_y(m,k)-y_Yd(j)).^2); end

G_Qij=zeros(14,14);

[G_Qi,G_Qj]=meshgrid(G_Qi,G_Qj); G_Qij=min(G_Qi,G_Qj);

%计算G点到Qi,Qj两点连线的距离G_Pij for i=1:14

for j=1:14

if sqrt((y_Yd(i)-y_Yd(j)).^2+(x_Yd(i)-x_Yd(j)).^2)

G_Pij(i,j)=abs((G_x(m,k)-x_Yd(i)).*(y_Yd(i)-y_Yd(j))-(G_y(m,k)-y_Yd(i)).*(x_Yd(i)-x_Yd(j)))./sqrt((y_Yd(i)-y_Yd(j)).^2+(x_Yd(i)-x_Yd(j)).^2); end end end

%计算PQij=min{PQi,PQj} P_Qi=sqrt(G_Qi.^2-G_Pij.^2); P_Qj=sqrt(G_Qj.^2-G_Pij.^2); P_Qij=min(P_Qi,P_Qj); % Z_g_wang=zeros(14,14);

%根据平面几何关系计算出由Qi,Qj点的深度线性外推的Pij的深度Z_p for i=1:14; for j=1:14;

if (Qi_Qj(i,j)==0)|(z_Ft(i)==z_Ft(j)) Z_p(i,j)=z_Ft(i); else ((P_Qi(i,j)z_Ft(j)))

Z_p(i,j)=z_Ft(j)-(abs(z_Ft(i)-z_Ft(j))).*P_Qij(i,j)./Qi_Qj(i,j); else ((P_Qi(i,j)

z_Ft(i)))

Z_p(i,j)=abs(z_Ft(i)-(abs(z_Ft(j)-z_Ft(i)).*P_Qij(i,j)./Qi_Qj(i,j))); else ((P_Qj(i,j)z_Ft(i)))

Z_p(i,j)=z_Ft(j)+(abs(z_Ft(j)-z_Ft(i)).*P_Qij(i,j)./Qi_Qj(i,j)); else ((P_Qi(i,j)z_Ft(j)))

Z_p(i,j)=z_Ft(i)+(abs(z_Ft(i)-z_Ft(j))).*P_Qij(i,j)./Qi_Qj(i,j); else ((P_Qi(i,j)z_Ft(j)))

Z_p(i,j)=z_Ft(j)+(abs(z_Ft(i)-z_Ft(j))).*P_Qj(i,j)./Qi_Qj(i,j); else if if if if if if

((P_Qi(i,j)z_Ft(i)))

Z_p(i,j)=z_Ft(i)+(abs(z_Ft(j)-z_Ft(i)).*P_Qi(i,j)./Qi_Qj(i,j));

else

Z_P(i,j)=1; end

end

end

end

end

end

end

end

end

%计算G点深度Z_g

Z_g_1=ones(14,14)./(G_Pij.^2+G_Qij.^2+Qi_Qj.^2); Z_g_2=Z_p./(G_Pij.^2+G_Qij.^2+Qi_Qj.^2);

for i=1:14

for j=1:14

Z_g_3(i,j)=sum(sum(Z_g_1));

end

end

Z_g(m,k)=sum(sum(Z_g_2./Z_g_3));

end

end

end

%程序一：插值并作海底曲面图

x1=75:1:200;

y1=-50:1:150;

[x1,y1]=meshgrid(x1,y1);

z1=griddata(x_Yd,y_Yd,z_Ft,x1,y1,"v4");

meshc(x1,y1,z1);

%程序二：插值并作出水深小于5的海域范围。

x1=75:1:200;

y1=-50:1:150;

[x1,y1]=meshgrid(x1,y1);

z1=griddata(x_Yd,y_Yd,z_Ft,x1,y1,"v4"); %插值

z1(z1>=5)=nan; %将水深大于5的置为nan，这样绘图就不会显示出来 meshc(x1,y1,z1)

图表 1插值海底曲面图

图表 2插值并作出水深小于5的海域范围