水道测量的设计与解法Word格式.docx
《水道测量的设计与解法Word格式.docx》由会员分享，可在线阅读，更多相关《水道测量的设计与解法Word格式.docx（11页珍藏版）》请在冰豆网上搜索。
。
但是一个加权平均并不能完全反映曲面的走势。
例如，再简单的一维情况下，当互相间隔4米的三个点A、B、C，若A点深8米，B点深4米，那么按照上面的算法就应该得到C点深度为4.8米，而这个4.8米显然与我们日常的认识：
C点深度应该小于4米，是相悖的。
从这个例子可以看出远点A对C点的影响大于B点而致使C点计算得到的深度大于4米。
所以需要在加入一项，使得近点的影响更大，于是引入趋向项以增强近点对于所要求的点的影响。
然而在二维平面上，G点并不在Q1、Q2的连线上，于是由G点向Q1、Q2的连线作垂线，设垂足为P，则Q1、Q2对G的影响可以在P上体现出来。
于是，再次引入关于P点的反距离加权影响，最终得表达式：
根据以上表达式可推测出：
1、每个已知点对未知点的影响表现在两者间距离的平方分之一；
2、任意一已知点对于位置点都有影响，即，未知点的水深随着任一已知点的变化而变化；
3、任意两点对于未知点的影响取决于其中的三个距离；
4、若要求一个已知点A对某一未知点的影响，那必须把除这个点之外的所有已知点与A两两组合求出对未知点的总影响。
（不足之处：
我们认为Z的左边与等式右边的式子一定是成正比的关系，但不一定是相等，即比例系数=1。
应当存在比例系数K，最后应当利用已知点，代入公式求出K。
）
以上所叙述的影响由以下几点原因产生：
1、两个已知点与未知点的距离；
2、未知点与两已知点连线的垂直距离；
3、两已知点之间的距离。
有以上的出的Z的坐标结合之前得到的X、Y坐标，经过双三次样条插值法便可拟合出曲面。
（2）极小曲率差分迭代法
在这一方法中，论文直接对区域作细密分割，只要求出每个节点上的水深值，就能直接描绘海底曲面。
以下皆为求出海水深度的讨论：
将分割好的区域中不在网格上的点作微小平移，论文认为，这样做不会影响最后结果。
将处理过的坐标代入微分方程：
这时有边界条件：
1、分割矩形区域四条边上作用力为零；
2、切线零混合矩；
3、u（xn,yn）=wn。
离散化海底曲面总平方曲率得：
当涉及二维时，可以用简单曲率逼近表达式：
最后让uij出的导数为零，再代入调整后的数据，即可解出uij，最后将坐标画在坐标纸上便可得到海底大致曲面。
（3）IDW解法
将X、Y的区域作等距离分割，由假设三，我们可以设每个点对于未知点的影响表现为距离的p次方，则
；
其中
，
，Di即为调整后的格点与未知点间的距离。
在上式中，指数p反映了未知点对于已知点的离散情况，p越大，说明两者离得越近，影响越大。
在得出了Z的坐标之后，再通过双三次样条插值法便可得出所要求图形。
程序与图形：
反距离平方格式与双三次样条拟合法程序：
#include<
stdio.h>
math.h>
floatdis（floatx1,floaty1,floatx2,floaty2）
{
floata=0;
a=（float）sqrt（（x2-x1）*（x2-x1）+（y2-y1）*（y2-y1））;
return（a）;
}
floatXP（floatx1,floatx2,floatx3,floaty1,floaty2,floaty3）
floatxu,xd,x;
xu=x3*（x1-x2）*（x1-x2）+（y1-y2）*（x2*（y1-y3）-x1*（y2-y3））;
xd=（x2-x1）*（x2-x1）+（y2-y1）*（y2-y1）;
x=xu/xd;
return（x）;
floatYP（floatx1,floatx2,floatx3,floaty1,floaty2,floaty3）
floatyu,yd,y;
yu=y3*（y2-y1）*（y2-y1）+（x2-x1）*（（x2*y1-x1*y2）+x3*（y2-y1））;
yd=（x2-x1）*（x2-x1）+（y2-y1）*（y2-y1）;
y=yu/yd;
return（y）;
voidmain（）
floatQ[3][14]={
{77,81,88,105.5,107.5,108.5,117.5,129,140,157.5,162,162,185.5,195},
{3,56.5,147,85.5,-81,28,-38.5,7.5,141.5,-6.5,84,-66.5,22.5,137.5},
{8,8,8,8,9,6,9,4,8,9,4,9,6,8}};
floatX[15]={75,77,81,88,105.5,107.5,108.5,117.5,129,140,157.5,162,185.5,195,200};
floatY[16]={-100,-81,-66.5,-38.5,-6.5,3,7.5,22.5,28,56.5,84,85.5,137.5,141.5,147,150};
floatG[3][240];
floatGP=0,PQ1=0,PQ2=0,Q1Q2=0,GQ1=0,GQ2=0,ZP=0,ZGU=0,ZGD=0,zgu=0,zgd=0,xp=0,yp=0;
intm=0,n=0;
for（inti=0;
i<
15;
i++）
{
for（intj=0;
j<
16;
j++）
{
G[1][m++]=Y[j];
G[0][n++]=X[i];
}
}
for（intk=0;
k<
240;
k++）
G[2][k]=0;
zgu=0;
zgd=0;
for（inti=0;
14;
for（intj=0;
{
GP=Q[0][j]*（Q[1][i]-G[1][k]）-Q[0][i]*（Q[1][j]-G[1][k]）-G[0][k]*（Q[1][i]-Q[1][j]）;
if（i!
=j）
{
xp=XP（Q[0][i],Q[0][j],G[0][k],Q[1][i],Q[1][j],G[1][k]）;
yp=YP（Q[0][i],Q[0][j],G[0][k],Q[1][i],Q[1][j],G[1][k]）;
}
GQ1=dis（Q[0][i],Q[1][i],G[0][k],G[1][k]）;
GQ2=dis（Q[0][j],Q[1][j],G[0][k],G[1][k]）;
if（GQ1>
GQ2）GQ1=GQ2;
Q1Q2=dis（Q[0][i],Q[1][i],Q[0][j],Q[1][j]）;
PQ1=0;
PQ2=0;
PQ1=dis（xp,yp,Q[0][i],Q[1][i]）;
PQ2=dis（xp,yp,Q[0][j],Q[1][j]）;
if（i==j）ZP=Q[2][j];
elseif（PQ1==PQ2&
&
PQ1!
=0）ZP=（Q[2][j]+Q[2][i]）/2;
elseZP=（Q[2][j]*PQ1-Q[2][i]*PQ2）/（PQ1-PQ2）;
if（ZP<
0）ZP=0;
ZGU=ZP*（float）pow（（GP*GP+GQ1*GQ1+Q1Q2*Q1Q2）,-1.0）;
zgu=ZGU+zgu;
ZGD=（float）pow（（GP*GP+GQ1*GQ1+Q1Q2*Q1Q2）,-1.0）;
zgd=zgd+ZGD;
}
G[2][k]=zgu/zgd;
FILE*fp;
fp=fopen（"
file2.txt"
"
w+"
）;
intl=0;
for（l=0;
l<
l++）
if（l==0）fprintf（fp,"
{%.3f"
G[2][l]）;
elseif（l<
239）fprintf（fp,"
%.3f"
elsefprintf（fp,"
%.3f}"
fclose（fp）;
反距离平方格式与双三次样条拟合法图像
散点图：
画等高线程序：
clc,clear
x1=[129.0140.0108.588.0185.5195.5105.5157.5107.577.081.0162.0117.0162.0];
y1=[7.5141.528.0147.022.5137.585.5-6.5-81.03.056.584.0-38.5-66.5];
z1=[48686889988499];
plot（x1,y1）;
x2=75:
1:
200;
y2=-50:
150;
z2=griddata（x1,y1,z1,x2,y2'
'
v4'
[C,h]=contour（x2,y2,z2,10）;
clabel（C,h）
colormapcool
等高线图样
IDW程序（VB）：
OptionExplicit
OptionBase1
PrivateSubCommand1_Click（）
Dima（）AsVariant,b（）AsVariant,c（）AsVariant,iAsInteger
a（）=Array（77,81,88,105.5,107.5,108.5,117.5,129,140,157.5,162,162,185.5,195）
b（）=Array（3,56.5,147,85.5,-81,28,-38.5,7.5,141.5,-6.5,84,-66.5,22.5,137.5）
c（）=Array（8,8,8,8,9,6,9,4,8,9,4,9,6,8）
Dimsum1AsSingle,sum2AsSingle,sum3AsSingle
DimxAsInteger,yAsInteger,zAsSingle,jAsInteger
sum1=0
sum2=0
Forx=75To200Step5
Fory=-100To150Step10
Fori=1To14Step1
sum1=sum1+1/（（x-a（i））^2+（y-b（i））^2）^（3/2）
Nexti
Forj=1To14
sum2=sum2+c（j）/（sum1*（（x-a（j））^2+（y-b（j））^2）^（3/2））
Nextj
Open"
C:
\数学建模\text1"
ForAppendAs#1
Print#1,Int（sum2*100+0.5）/100
Close1
List1.AddItemInt（sum2*100+0.5）/100
Nexty
Nextx
EndSub
画图程序：
clc
clear
clf
x=75:
5:
y=-100:
10:
Z=[];
%z内数据有以上VB程序获得。
见text1
z1=reshape（z,26,26）;
xi=75:
0.5:
yi=-100:
zi=interp2（x,y,z1,xi'
yi,'
spline'
mesh（xi,yi,zi）
holdon
X=（75:
200）;
Y=（-100:
150）;
Z=-5*ones（251,126）;
mesh（X,Y,Z）
Z
由IDW法所的图像
升级会员 