前言
本文仅为学习matlab过程中的一些记录,奈何才疏学浅,如有错误希望各位大佬不吝赐教。
一、解决的问题
在研究河相关系公式需要率定河相关系系数时,需要先得到自然河宽B和平均水深H等基本参数,其中平均水深H又可以由河道断面面积A/自然河宽B得到,因此为了根据已有的河道各断面的地形数据获得不同水位下的自然河宽、断面面积及断面最大水深编写了以下的代码。
二、所需数据
输入的基础数据包括两个:
1、河道各个断面的地形数据。它反映了某一断面河床横断面形态,由一个个地形点构成,地形点包括两个维度的参数分别是到左岸的起点距和高程;
2、设计水位数据。它反映了每个断面所研究的若干个设计水位;
代码中的数据下载地址如下:
链接:https://pan.baidu.com/s/114scyBlUO-XzzPz2Pr6APQ
提取码:bum1
三、解决思路
1、河道的断面面积可以利用matlab的积分函数trapz()来计算。
2、问题的难点在于如何得到不同水位下某一断面的自然河面宽度,采用的思路是:取左岸低于水位的第一个地形点和右岸低于水位的第一个地形点,取两个地形点之间的起点距差值近似计算自然河面宽度B(见图)。
可以看出这是一种近似方法不是十分精确,但是考虑到地形点数据足够密集,河道宽度较大,因此结果的误差在可以接受的范围之内。
3、断面最大水深则可以利用matlab中max()函数;
四、matlab代码
代码如下:
%2021/8/12
%1、数据导入
%2、求解各断面不同水位下自然河宽及断面面积
%3、数据导出
clc;clear
%% 数据导入
excel_path='D:\个人文件夹\matlab学习\test.xlsx';
M=xlsread(excel_path,1); %导入地形数据
N=xlsread(excel_path,2); %导入设计水位数据(表头不读入)
[MX,MY]=size(M); %获得矩阵M的行列数,行数MX,列数MY
[NX,NY]=size(N); %获得矩阵N的行列数,行数NX,列数NY
%% 求各断面在不同水位下自然河宽及断面面积
m=0.5*MY; %m为总的断面数
CS=[]; %定义空矩阵以存储所有断面自然河宽B、断面面积A、断面最大水深Hm数据
LC=[]; %定义空矩阵以存储航道里程数据
%各断面循环
for k=1:m
x=2*k-1;
lc=M(1,x); %获得各断面的航道里程
LC=[LC;lc]; %存储所有断面的航道里程数据
C=[]; %定义空矩阵以存储某一断面B、A、Hm数据
%各设计水位循环,求一个断面的B、A、Hm
for j=1:NY
h=N(k,j); %各设计水位依次赋值给h,(第一行为表头)
%寻找左岸水面与地形交点
ls=[]; %定义空矩阵以存储左岸水面与地形交点数据
for i=2:MX %自左到右地形点循环,从2开始是去除了表头
z=M(i,x+1); %各地形点高程依次赋值给z
if (z-h)<=0 %以水位以下的第一个地形点作为水面与地形交点
l=i;
break
end
end
ls=[ls l]; %存储左岸水面与地形交点数据
%寻找右岸水面与地形交点
rs=[]; %定义空矩阵以存储右岸水面与地形交点数据
for i=MX:-1:2 %自右到左地形点循环
z=M(i,x+1);
if (z-h)<=0
r=i;
break
end
end
rs=[rs r]; %存储右岸水面与地形交点数据
%获得水深数据
H=[]; %存储水深数据
for i=l:r
ho=h-M(i,x+1); %求解水深,即为水位-地形高程
if ho<0 %去除河道中凸出水面的浅滩等的影响
ho=0;
end
H=[H,ho]; %存储水深数据,行向量形式
hm=max(H); %获取最大断面水深hm
end
b=M(r,x)-M(l,x); %求解自然河宽,即为右岸交点-左岸交点
X=M(l:r,x); %X为起点距数据
H=H';
a=trapz(X,H); %求解断面面积,采用梯形积分
c={b a hm ' '}; %整合自然河宽、断面面积及最大水深数据
C=[C c]; %存储某一断面B、A、Hm数据
end
CS=[CS;C]; %存储所有断面B、A、Hm数据
end
%% 数据导出
str={'航道里程','枯水B','枯水A','枯水Hm','备注','中水B','中水A','中水Hm','备注','洪水B','洪水A','洪水Hm','备注'}; %表头
CST=[str;[num2cell(LC),CS]];
xlswrite(excel_path,CST,3);