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

%数据输入x_Yd=[129.0140.0108.588.0185.5195.5105.5157.5107.577.081.0162.0162.0117.5];y_Yd=[7.5141.528.0147.022.5137.585.5-6.5-81.03.056.5-66.584.0-38.5];z_Ft=[48686889988949];%由已知14点计算出14*14点（其中有两点重复）%计算QiQj的值Qi_QjQi_Qj=zeros(14,14);fori=1:14forj=1:14Qi_Qj(i,j)=sqrt((x_Yd(i)-x_Yd(j)).^2+(y_Yd(i)-y_Yd(j)).^2);endend[G_x,G_y]=meshgrid(x_Yd,y_Yd);%计算所求196个G点的坐标%求出每个G点对应的深度Z_gform=1:14fork=1:14if((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);fori=1:14G_Qi(i)=sqrt((G_y(m,k)-y_Yd(i)).^2+(G_x(m,k)-x_Yd(i)).^2);endforj=1:14G_Qj(j)=sqrt((G_x(m,k)-x_Yd(j)).^2+(G_y(m,k)-y_Yd(j)).^2);endG_Qij=zeros(14,14);[G_Qi,G_Qj]=meshgrid(G_Qi,G_Qj);G_Qij=min(G_Qi,G_Qj);%计算G点到Qi,Qj两点连线的距离G_Pijfori=1:14forj=1:14ifsqrt((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);endendend%计算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_pfori=1:14;forj=1:14;if(Qi_Qj(i,j)==0)|(z_Ft(i)==z_Ft(j))Z_p(i,j)=z_Ft(i);elseif((P_Qi(i,j)=P_Qj(i,j))&&(Qi_Qj(i,j)=P_Qi(i,j))&&(z_Ft(i)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);elseif((P_Qi(i,j)P_Qj(i,j))&&(Qi_Qj(i,j)P_Qj(i,j))&&(z_Ft(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)));elseif((P_Qj(i,j)=P_Qi(i,j))&&(Qi_Qj(i,j)=P_Qi(i,j))&&(z_Ft(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));elseif((P_Qi(i,j)=P_Qj(i,j))&&(Qi_Qj(i,j)=P_Qj(i,j))&&(z_Ft(i)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);elseif((P_Qi(i,j)=Qi_Qj(i,j))&&(P_Qj(i,j)=Qi_Qj(i,j))&&(z_Ft(i)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);elseif((P_Qi(i,j)=Qi_Qj(i,j))&&(P_Qj(i,j)=Qi_Qj(i,j))&&(z_Ft(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));elseZ_P(i,j)=1;endendendendendendendendend%计算G点深度Z_gZ_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);fori=1:14forj=1:14Z_g_3(i,j)=sum(sum(Z_g_1));endendZ_g(m,k)=sum(sum(Z_g_2./Z_g_3));endendend%程序一：插值并作海底曲面图x1=75:1:200;y1=-50:1:150;[x1,y1]=meshgrid(x1,y1);z1=-griddata(G_x,G_y,Z_g,x1,y1,'v4');meshc(x1,y1,z1)xlabel('x'),ylabel('y'),zlabel('z')figure;%程序二：插值并作出水深小于5的海域范围。x1=75:1:200;y1=-50:1:150;[x1,y1]=meshgrid(x1,y1);z1=-griddata(G_x,G_y,Z_g,x1,y1,'v4');%插值z1(z1=-5)=nan;%将水深大于5的置为nan，这样绘图就不会显示出来meshc(x1,y1,z1)xlabel('x'),ylabel('y'),zlabel('z')%产生的G点中有重复点，运行时有警告