数学建模题 年降雨量计算

组号183

B题、中国水坝对区域降水的影响1.摘要：

本文通过建立数学模型研究了中国水坝对区域降水影响问题。对于气象空间站分布不均匀，使得中国大陆平均降雨量不能直接计算，并且很难得到某地区非常准确的降雨量数字，我们采用根据距离加权来计算某一点的降雨量，根据距离它最近的m个点来计算该点的降雨量。在建立模型求解中，我们着重解决了以下问题：1、用matlab编程处理所给xls信息；2、借助c++实现我们做的模型，并进行稳定性测试。3、将算法移植到matlab上，解出精确度为1度的地图上的点的降雨量信息。4、借助matlab将中国地图大致范围求出。5、分析某地区的降雨量变化

声明：由于原始数据坐标问题，导致画出图像与真实情形相差太大，故借助matlab将错误数据更正。

2.问题重述

根据附件中的材料，研究中国水坝对区域降水的影响。

建立相应的数学模型，并解决的如下问题：

1.估计1951年——2008年中国大陆的年平均降水量；

2.估计1951年——2008年某一地区的年降水量，即给出某一地区

的经度和纬度，用所建模型计算出该地区的年降水量。按照你的

方法，估计水坝地区的降水量（1951年——2008年）。

3.研究中国水坝对区域降水的影响。（注：影响可能是多方面的。

可能会增加某地区的降水，也可能会减少另一地区的降水，还

可能会对某一地区的降水无影响。请大家从多个层面考虑这个问

题。）

3.基本假设

a)假设经过修改的数据真实可靠。

b)假设大坝是平均分布在全国各地的。

c)假设大坝没有因年代久远或水量过大而影响蓄水量，并且一直完好如初。

4.符号说明:

m为距离任意点(x,y)最近的点的个数

未知点(x,y)的降雨量

为已知点的年平均降雨量

为第i个已知点第j年的降雨量

为m个最近点中第i个点与任意点(x,y)的距离

为第i个计算出来的点的降雨量，

n为计算过的点的个数。

5.术语说明：

已知点预测：在验证求未知的是否准确的时候，假设一个离已知点很近的点为未知点，求出它的降雨量，与刚取的已知点比较，看差距大小。

下文提到的c++程序只有一个，就是附录3中给的

6.模型的建立与求解

6.1模型的建立：

由题目中附件3可以看出，气象站在全国并不是平均分布的，所以不能用加起来求平均值的方法，我们利用距离位权法建立了数学模型，以求出任意一点的平均降雨量。

设任意一点(x,y)降雨量为R(x,y)则：

其中：

为距离任意点(x,y)最近的点的个数

为已知点的年平均降雨量

r j为第i个已知点第j年的降雨量

w i为m个最近点中第i个点与任意点(x,y)的距离

m点的取值和R(x,y)的精确度有关，若m很大，则会包括所有城市，虽然进行已知点验证时很精确，但不符合实际情况，若m很小，则精确度会下降，关于m的取值，将会在下边的可靠性分析中讨论。

6.2模型可靠性分析：

根据利用c++编出来的程序，可以验证，当m>60时，进行已知点验证，与原降雨量差距很小，但是不符合实际，因为某地区降雨量不会和很远距离的降雨量有太大相关性。根据c++程序验证，取m=15。

本模型对于气象站分布较密集的地方精确度较高，但对于西部地区气象站分布不均且数量有限情况下，可靠性会下降。从c++程序来

看（去掉70行处注释符），当m=15时进行已知点预测的差别大的主

要在编号140以后的地区。

6.3问题求解：

6.3.1问题一的解:

借助matlab将数据网格化大致算出中国降雨量可能会覆盖到的地方如附件2。如图1

纬度

图1

经度

得到了中国大致的限制方程：

-0.72*x0+94.72-y0<0

其中x0,y0为当时要构造的点的坐标。

年平均降雨量R总为：

其中：

为第i个计算出来的点的降雨量，n为计算过的点的个数。

由此，年平均降雨量R总求出

图2为求出的全国降雨量的分布（精确到1度）程序在附录3中

图2全国降雨量分布

6.3.2问题二的解:

同模型建立过程。

6.3.3问题三的解：

我们选取全国1个地区作为我们的分析对象：东北区（1-34）。采用所建立的模型，借助一元线性回归来分析降雨量变化。

图2东北地区平均降雨量分布

图3东北地区降雨量逐年分布水平

利用matlab算出每年东北地区平均降雨量的一次拟合曲线，再不考虑人为因素时得到初步结论：水坝的修建会减少降雨量。

然后利用matlab程序计算东北地区年降水量和大坝修建的相关系数（附录4）

求得相关系数为-0.114，可以看出东北地区的降雨量和水坝的修建基本无关。由下图

上图可以看出，大坝在1960年附近和2000年附近时候有大幅上升，但东北地区的降水量波动不是很大。故得到最终结论：东北地区降雨量

和全国大坝修建情况无关。

7.参考资料

8.附件

Matlab调试环境：2010a

C++调试环境：vs2008

附录1：以下程序用来初始化：

clear;

%%%%%%%初始化%%%%%%%%%

xls = xlsread('2009A2.xls');

dam = xlsread('2009A1.xls');

for i = 2: 161

x(i-1) = xls(i,3);

y(i-1) = xls(i,4);

end

for i = 2:161

for j = 5:62

zz(i-1,j-4) = xls(i,j);

end

end

for i = 1:160

z1(i) = zz(i,1);

z2(i) = zz(i,58);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%大坝容量年增加量%%%%%%%%%%%

damYearX = 1:2009;

damYearY = zeros(2009,1);

damX = zeros(4607,1);

damY = zeros(4607,1);

water = zeros(1,58);

year = 1951:2008;

for i = 1:4607

damX(i) = dam(i,1);

damY(i) = dam(i,2);

end

for i = 1:4607

damYearY(damX(i,1),1) = damYearY(damX(i,1),1)+damY(i,1);

end

for i = 1:2008

if i > 1950

water(i-1950) = damYearY(i);

end

end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

附录2，画图构造降雨量覆盖图

xtmp = linspace(min(x),max(x),80);

ytmp = linspace(min(y),max(y),80);

[X,Y] = meshgrid(xtmp,ytmp);

Z1 = griddata(x,y,z1,X,Y);

Z2 = griddata(x,y,z2,X,Y);

%mesh(X,Y,Z1);

mesh(X,Y,Z2);

附录3

#include

#include

#include

#include

#define MAX 200

using namespace std;

double sum(double*a,int m)

{

double ans = 0;

for(int i=0;i

ans+=a[i];

return ans;

}

int main()

{

ifstream cin("aa.txt");//aa.txt文件有MATLAB生成程序附录(5)，将aa.txt放到工程文件夹下

double x[MAX],y[MAX];

double x0,y0;

int i,j,m;

double save[MAX],rainPerSite[MAX];

int min[MAX];

for(i=0;i

double tmp[MAX];

//scanf("%d%lf%lf",&m,&x0,&y0);

for(i=0;i<160;i++)

cin>>x[i];

for(i=0;i<160;i++)

cin>>y[i];

for(i=0;i<160;i++)

cin>>rainPerSite[i];

for(m = 2;m<150;m++)

{

int count=0;

cout <<"m= " << m << " ";

for(int k=0;k<160;k++)

{

x0 = x[k]-0.01;

y0 = y[k]-0.01;

for(i=0;i<160;i++)

save[i]=sqrt(fabs((x[i]-x0)*(x[i]-x0))+fabs((y[i]-y0)*(y[i]-y0)));

for(i=0;i

{

for(j=0;j<160;j++)

if(save[min[i]]>save[j])

{

min[i] = j;

tmp[i] = save[min[i]];

}

save[min[i]] = 100000;

}

for(i=0;i<3;i++){//列出了最近的3个点的信息，去掉注释符可看见

//printf("%.2lf ",rainPerSite[min[i]]);

//printf("%.2lf ",tmp[i]);

//printf("%d ",min[i]);

}

//cout <

double ans = 0;

for(i=0;i

ans += (1-tmp[i]/sum(tmp,m))*rainPerSite[min[i]];

//cout <

if((1-fabs((rainPerSite[0]*.7+rainPerSite[1]*0.3)-ans)/ans)>0.1)count++;

//cout<<(1-fabs((rainPerSite[0]*.7+rainPerSite[1]*0.3)-ans)/ans)<

//if(m==15)cout<<"站点: "<

}

cout <<"与原点相差大于10%的点的个数" <

}

return 0;

}

附录4

function re=res(x,y,n)

%RES 求的是最小二乘法的相关系数r，其中%x，y为数据，n是数据的个数

lxx=sum(x.^2)/n-(sum(x)/n)^2;

lyy=sum(y.^2)/n-(sum(y)/n)^2;

lxy=sum(x.*y)/n-sum(x)/n*sum(y)/n;

re=lxy/(sqrt(lxx*lyy));

end

附录5

以matlab程序下用来写入记事本供c++用z90 = rot90(zz);

a = sum(z90);

a= a./58;

fp = fopen('aa.txt','wt');

for i =1 : 160

fprintf(fp, '%.2f ', x(i));

end

fprintf(fp,'\n');

for i = 1:160

fprintf(fp,'%.2f ',y(i));

end

fprintf(fp,'\n');

for i = 1:160

fprintf(fp,'%.2f ',a(i));

end

fclose(fp);