1、1土地平整问题摘要随着社会的快速发展，对于山区的土地整理显得格外重要。为了更好地合理利用土地，选择合适的土地整理地点，达到节约成本的目的。我们会建立模型，根据几个重要的制约因素，使得土石方量最小，而且成本最低，在实际的操作过程中更接近于科学实施。针对问题一：利用软件 MATLAB 中的绘图功能，对数据进行处理，得到这片土地的三维图形与等高线图。针对问题二：采用二重积分的意义，即曲顶柱体的体积。利用分割、取值、求和、取极限的思想，将山地分割成无数个小曲顶柱体，然后对这些小曲顶柱体进行求和，并结合土石方量最小，挖填土时费用使用的最少限制条件，从而求得在什么海拔高度进行平整。对于平面位置的确定，利用
2、枚举法，总费用最少，求得在什么位置，什么海拔高度平整，已达到最优的标准。针对问题三：类比于问题二，采用相同的方法，在不同的海拔高度上平整，已满足二层的台阶状地块的要求，达到挖填土石方量最小的目标。最后，对所建立的模型和求解方法的优缺点给出了客观的评价，并指出了改进的方法。关键词：土石方量最小；二重积分；枚举法；2一、问题重述与分析1.1 问题的重述十堰市是一个山区城市，向山要地是十堰市发展的一个必然的选择，但是如何在一片山地之中选择合适的方位与开挖深度，从而使总的土石方量最小，就是一个有意义的命题了！某工厂为了在一片长度为 1500 米，宽度为 1000 米的山地之中，开挖出一个 800 米6
3、00 米平坦连续的长方形地块作为工厂的厂房地基，前期已经在本块土地上测量出长、宽每隔 30 米的网格的对应网格点的海拔高度。问题：用附件（1）中的数据画出工厂的这片土地的三维图形与等高线图；从什么地方，什么海拔高度平整一块 800 米600 米的连片土地能使总的土石方量最小？如果允许平整出来的土地为二层的台阶状地块，要求各地块的长、宽不少于 60 米，又将从什么地方、什么海拔高度分别开挖，能使总的土石方量最小？1.2 问题的分析1.2.1 问题一我们可以利用 MATLAB 软件直接调用 surf 函数和 contour 函数分别画出这片土地的三维图形和等高线图。 1.2.2 问题二考虑要平整一
4、块 800 米600 米的连片土地，且要使土石方量最小，我们需要确定这块 800 米600 米的连片土地的底面投影区域，以及海拔高度，使得其土石方量最小。让底面投影区域在长度为 1500 米，宽度为 900 米内枚举，此时计算出对应的土石方量最小的是，在这些所有的最小值中取得最小值的那块投影区域即为所求。我们采用此种方法求解。在平整土地的过程中，凸出的地方是要挖土的，凹下的地方是可以填土的，这样挖下来的土可以填充到凹下去的部分。通过计算体积的方式计算挖土和填土大小，求出最小值对应的海拔高度，即可确定整连片土地的具体高度。二、基本假设1、在平整土地的过程中，有些地方是要挖山的，但有些地方是要填土
5、的，假设填土的每立方米所需的费用为挖山的每立方米土石方所需费用的 1/3。2、假设除了挖土和填土以外，在平整土地的过程中其他作业（如登高）产生的费用都与 800 米600 米的连片土地所处的位置方向和海拔高度均无关。这样我们只需考虑挖土和填土的土石方量及费用，据此来考虑连片土地所处的位置方向和海拔高度。3、假设计算体积的过程中，分割的小曲顶柱体不能达到无穷小，而产生的误差忽略不计。实际计算中我们取一个很小的步长去划分，使其划分尽可能的3小，产生的误差忽略不计。三、符号约定四、模型的建立4.1. 问题一求解1利用 MATLAB7.0 软件的三维绘图功能，画出工厂这片土地的三维图形如图 4.1.1
6、 所示：符号 说明L平整连片土地块的长W平整连片土地块的宽D平整连片土地块的底面所在闭区域i平整连片土地块底面划分的子区域i子区域 的面积iDwV平整土地所需的挖土量f平整土地所需的填土量C平整土地挖土和填土所需要的总费用4图 4.1.1山地的三维图形2、这片土地的三维等高线图如图 4.1.2 所示，二维等高线图 4.1.3 所示：图 4.1.2山地的三维等高线图5图 4.1.3山地的二维等高线图4.2.问题二模型建立4.2.1. 平整块的海拔高度的确定通过土石方量费用最小的原则，确定平整土地海拔的开挖高度。平整土地快所对应的的底面区域记为 ，显然 是一个 800 米600 米的矩形区域。设D
7、顶部海拔高度是一个非负连续函数 ， 。假设 的位置已经确),(yxfD,定，下面我们来讨论在什么样的海拔高度下，土石方量费用最小。考虑底面区域 对应的山体是一个曲顶柱体，把区域 分成 个小区域n， ， ， ，为了简化运算，我们分成 个等大的区域，每个区域1D2nn都是边长为 的正方形。再分别以这些小区域的边界为准线，作母)(iid线平行于 轴的柱面，这些柱面将原曲顶柱体分成 个细的曲顶柱体。每个区z域 上任取一点 。当 很大时，所分的区域 很小，边长 的值很小，i ),(iniDd每个区域 对应的海拔高度 z 可以近似为一个相同的值 ，记所有小区iD ),(if域 的最小海拔高度为 ，最大海拔
8、高度为i 1h2h6（1）),(mini1fh（2）axi2ni设在平整块连片土地建在海拔高度为 处 。当小区域 所对应z21hiD的海拔高度 时，说明这块区域需要挖土，挖土量为 ，zfi),( iV（3）iiwzfVi ),(这里 表示 的面积。 时，说明这块区域需要填土，填土量为iiDi,fV（4）iiffzVi ),(总的挖土量 ，w（5）ki iikiwzfi11),(其中 为满足 条件的小区域 的个数。总的填土量 ， kzfi),(iDfV（6）li iilff fzVi11i ),(其中 为满足 条件的小区域 的个数。kzfi),(i当 ，挖土量大于填土量，这样挖出来的土足够填充凹
9、下的山地。总fwV的土石方量费用为：（7） ljjki ii fzzfC1j1f ),(3),(3 当 ，挖土量小于填土量，这样挖出来的土不足以填充凹下的山地。fwV需要从别的地方挖土来填充。总的土石方量费用为：（8） ljjfwfwf fzVC1j),(34)(313 因此7（9）fwljjki ii Vfzzf ，1j1 ),(3),( fwljj Vfz，1j),(34C如果让平整块的底面投影区域取遍所有可能的值，再根据上述思想求解此底面投影区域对应土石方量费用最小 的开挖海拔高度 。按照最小的 对应minChminC的底部和开挖海拔高度来开挖，所得平整的连片土地能使总的土石方量最小。4
10、.2.2. 平整块的底面区域的确定将平整块区域投影到底面（ 平面） ，确定开挖方向，即要确定平整块的xoy底面位置。显然底面是 800 米600 米的矩形，假设矩形的四个顶点分别为。山地是长度为 1500 米，宽度为 900 米的矩形区域，因此平整块DCBA、底面投影矩形 可以在 1500900 的区域范围内任意移动，四个顶点不能超出这个范围我们要取遍所有可能的位置，必须枚举所有的可能值。设矩形 的边长为 = =600， = =800，且各点坐标位置也已经ABWCL确定，设各顶点的坐标分别为 （ ， ） 、 （ ， ） 、 （ ， ） 、axyBbxyCcxy（ ， ） 。设 边的中点为 ，
11、边的中点为 ，那么唯一确定了线段DdxyEDF。EFO xyADBCEF由中点公式，中点 ， 的坐标分别可以表示为EF（ ， ） ， （ ， ） ，2bax2bayF2dcx2dcx8且 的边长为 = =800,因此只要矩形 的位置和大小确定，那么EFACABCD的大小和位置就唯一确定。反过来如果 的位置和大小确定了（ 分别EFFE,为短边的中点） ，那么矩形 的位置和大小也唯一确定。下面我们分别说明B矩形四个顶点坐标，以及四个边的直线方程计算方法。设 =800， 两点的坐标分别为 （ ， ） 、 （ ， ） ，为了EF、 exyfxfy使计算不重复，保证 的唯一性，我们不妨设 ，且当 时,
12、ffe,若 EF 所在直线的斜率存在，设为 。当 时，斜率存在，fexkefxefxy分以下四种情况讨论：(1)当 ,即 平行于 轴时，此时矩形如下图所示：efxEFyxyOADBCEF ）（ 0,Pyx此时，四条边对应的直线方程分别为： ， ： ， ： ， ： ，ABeyC2WxeCDfyA2Wxe四个顶点的坐标分别为：（ ， ）, （ ， ） ， （ ， ）, （ ，2xeeBeexefyDe）fy矩形 的充分必要条件是：),(0yxPACD9220Wxxeefey0(2)当 时，0k）（ 0,PyxxyOADBCEF此时，四条边对应的直线方程分别为: ,即： ABeyxy)(k10eky
13、x: ,即：CDff ff： 0122ekxyWyx: Bk任意两条直线的交点即为顶点， 、 、 、 四个顶点的坐标分别由以ABCD下方程组解出： 0ekyxfexWyk12100ekyxfexWk120ffkyx122exk0ffkyx122ekxWk矩形 的充分必要条件是：),(0yxPABCD0122ekxyWykxekff 0122ekxyykx(3)当 ，即 垂直于 轴时0kEF11xyO）（ 0,PyxADBCEF此时，四条边对应的直线方程分别为： ， ： ， ： ， ： ，ABexC2WyeCDfxA2eWy四个顶点的坐标分别为：（ ， ）, （ ， ） ， （ ， ）, （ ，
14、2eeBxeeefDxe）fy矩形 的充分必要条件是：),(0yxPACD220Wxxeefey0(4)当 ，0k12）（ 0,PyxxyoDABCEF此时，四条边对应的直线方程分别为: ,即： ABeyxy)(k10ekyx: ,即：CDff ff： 0122ekxyWyx: Ak任意两天直线的交点即为顶点， 、 、 、 四个顶点的坐标分别由以ABCD下方程组解出： 0ekyxfexWk120ekyxfexyk12130ffkyx122ekxWk0ffkyx122exk矩形 的充分必要条件是：),(0yxPABCD0122ekxyWykxekff 0122ekxyykx如果让 在平面内取遍所
15、有可能的值，再根据上述的思想求解可以确定平EF整块的底面投影区域的位置，再根据 4.2.1 的介绍求解海拔位置。为了让能取遍所有可能值，先在山地海底平面上任意取一点作为 点，再根据E点确定 点的位置。这样只需考虑 点取遍所有值的可能情况，根据 EF 的长E度，确定 点可能的有效位置。平面区域内任意取一点 ， ， ，以 为圆),(eyx90ex150ey心， = =800 为半径作圆，显然 F 可能取得圆周上的所有点。RL14O xy EnF1F3F2为了避免重复，令 ，且当 时, ,这样 的取值为上图所efyefyefxF示的圆弧上。这样线段 就可以确定下来。EF4.3 问题三模型的建立我们可
16、以继续使用问题二中的模型，但要对问题二中的模型结合问题进行修改。我们依然要确定平整块的海拔高度，而对平整块的底面区域的面积的求解我们依然采用问题二中的模型。由于平整出来的土地为二层的台阶状地块，所以我们可以想象成两个曲柱体，这两个曲顶柱体所对应的高度不同，所以我们需要求出两个海拔高度。记所有小区域 的最小海拔高度为 ，最大海拔高度为iD1h2h（10）),(mini1fh（11）axi2ni设在平整块连片土地建在海拔高度为 处，为了求解方便，我们不妨设21z、，若 。当小区域 所对应的海拔高度 时，21z21hziD12i),(zf说明这块区域需要挖土，挖土量为 。isumV（12）2i1i
17、),(),( iiw zfzfsumi 这里 表示 的面积。 时，说明这块区域需要填土，填土量iiD2i为 。fsumV（13）2i21i1 ),(),( iif fzfzVisum 总的挖土量 ，w15（14） ki iikiikiw zfzfVi 12111 ),(),( 其中 为满足 条件的小区域 的个数。总的填土量 ， 2i),(zfiDfV（15） li iiliilff fzfi 12211i ),(),( 其中 为满足 条件的小区域 的个数。k2i),(zfi当 ，挖土量大于填土量，这样挖出来的土足够填充凹下的山地。总fwV的土石方量费用为：（16） ljjljjki iikii
18、 fzfz zC1j221j 1211f ),(3),(3 , 当 ，挖土量小于填土量，这样挖出来的土不足以填充凹下的山地。fwV需要从别的地方挖土来填充。总的土石方量费用为： ljj ljjfwfwfz fzVVC1j22 1j1),(34 ),(34)(31 （17）因此（18）fwljjki ii Vfzzf ，1j111 ),(3),( fwljj Vfz，1j),(34C五、模型的求解5.1 问题二的求解根据 4.2.1 和 4.2.2 介绍的方法，编写 MATLAB 程序求解。值得说明的是，在计算过程中，我们根据已有的检测数据，用了插值的方法，多维插值拟合计算某一点的海拔高度，如线
19、性插值，三次样条插值等，计算思想在数值分析中有说明，在这里不做细说，直接调用 MATLAB 中的多维插值函数 interp2()求解即可。16根据算法的基本思想，作出程序流程图见附录 1，MATLAB 程序见附录 2.程序的运性结果为， ， ， ，A)765.0 2.1,(B)495.0 37.,(C).487.9,120(，海拔高度 ，即在以这样的 构成的矩形区域D).4362.1,79( hABD为底面，在海拔 米处 平整一块 800 米600 米的连片土地，所需要的土6石方量费用最低。5.2 问题三的求解程序的运性结果为， ， ， ，海拔1A)657.0 .,(1B)4.0 .,(1C)
20、.572.0,48(1D).436.0,2(高度 ， ， ， ，3.9h2A8. 5.36,B. .,2C.5,982D，海拔高度 ，即在以这样的 构成的矩形区域)24.,( 2h1BA为底面，在海拔 米处和在以 构成的矩形区域为底面，在海拔.722CA65.76 处平整一块土地所需要的土石方量费用最低。六、模型的评价与推广6.1 模型的优点本题中的模型使用了 matlab 软件进行了拟合，具有较高的可信性。同时模型中对于问题分了多种情况进行讨论，具有很高的参考价值。6.2 模型的缺点本题中的模型求解过程比较繁琐，运算量较大，运算效率相对较低。6.3 模型的推广本题中的模型可以应用于其它土地平
21、整问题中，只需将模型进行修改便可以推广。17参考文献1 陈纪修，於崇华，金路，数学分析 第二版，M，北京：高等教育出版社，1999. 2 李庆扬，王能超，易大义，数值分析 第五版，M，北京：清华大学出版社，2008. 3 韩明，张积林，李林，林杰等，数学建模案例，M，北京：同济大学出版社，2012. 4 侯进军，肖艳清，谭敏，高明柯，数学建模方法与应用，M，南京：东南大学出版社，2012. 5 宋世德，郭满才，数学实验，M，北京：高等教育出版社，2002. 6 姜启源等，数学模型，M，北京：高等教育出版社，2003. 7 盛聚等，概率论与数理统计，M，北京：高等教育出版社，1983. 8 薛定
22、宇，陈阳泉，高等应用数学问题的 MATLAB 求解，M，北京：清华大学出版社，2004。 9 何正风，MATLAB 在数学方面的应用，M，北京：清华大学出版社，2012.18%脚本文件 drawpicture.m%用附件中的数据画出工厂的这片土地的三维图形与等高线x=0:30:900;y=0:30:1500;X,Y=meshgrid(x,y);Z=xlsread(data,1,b3:az33);surf(X,Y,Z); %三维图xlabel(横坐标),ylabel(纵坐标),zlabel(海拔高度);title(山地三维图);figure,contour(X,Y,Z); %二维等高线title
23、(二维等高线图);xlabel(横坐标),ylabel(纵坐标),zlabel(海拔高度);figure,contour3(Z); %三维等高线title(三维等高线图);xlabel(横坐标),ylabel(纵坐标),zlabel(海拔高度);%函数文件 getypos.m%功能：根据已知 E 点，及 F 点的 x 坐标，取得满足 EF=800 且 yfye 所对应的 F 点的坐标function yf=getypos(xe,ye,xf)%输入:E 点的坐标(xe,ye),F 点的 x 坐标 xf%输出:F 的 y 坐标 yfyf=sqrt(800*800-(xe-xf)2) + ye;%函
24、数文件 isvalidpoint.m%功能：判断 P 是否是有效的点function ret = isvalidpoint( p )% p 是向量，长度为 2，分别存放 x,y 轴的坐标if (p(1) 0) k -1; %AB 和 AD 的直线方程系数b1 = pe(1)+k*pe(2),(-sqrt(1+k*k)*width/2-pe(2)+k*pe(1); %AB 和 AD 的直线方程对应的常数项pa = A1b1;A2 = 1 k; k (-1); %AB 和 BC 的直线方程系数b2 = (pe(1)+k*pe(2),(sqrt(1+k2)*width/2-pe(2)+k*pe(1)
25、; %AB 和 BC 的直线方程对应的常数项pb = A2b2;A3 = 1 k; k -1; %CD 和 BC 的直线方程系数b3 = pf(1)+k*pf(2),sqrt(1+k2)*width/2-pe(2)+k*pe(1); %CD 和 BC 的直线方程对应的常数项pc = A3b3;A4 = 1 k;k -1; %CD 和 AD 的直线方程系数b4 = pf(1)+k*pf(2),(-sqrt(1+k2)*width/2-pe(2)+k*pe(1); %CD 和 AD 的直线方程对应的常数项pd = A4b4;else if k pf(1) | (y0 (pf(2) + width/
26、2)ret = 0;endelse if x0 = pf(1) %EF 平行于 x 轴if (x0 (pf(1) - long/2) | (y0 pf(2) )ret = 0;endelsek = (pf(2) - pe(2)/(pf(1) - pf(1);if k 0ab = x0 + k*y0 - pe(1) - k * pe(2);cd = x0 + k*y0 - pf(1) - k * pf(2);ad = -k*x0 + y0 - sqrt(1 + k2)*width/2 - pe(2) + k * pe(1);bc = -k*x0 + y0 + sqrt(1 + k2)*width
27、/2 - pe(2) + k * pe(1);if (ab 0) | (cd 0)ret =0;endif (ab 0) | (cd 0)ret =0;endelse if k 0) | (cd 0)ret =0;endif (ab 0) | (cd 0)ret =0;endendendendendend21ad = -k*x0 + y0 - sqrt(1 + k2)*width/2 - pe(2) + k * pe(1);bc = -k*x0 + y0 + sqrt(1 + k2)*width/2 - pe(2) + k * pe(1);if (ab 0) | (cd 0)ret =0;en
28、dif (ab 0) | (cd 0)ret =0;endelse if k 0) | (cd 0)ret =0;endif (ab 0) | (cd 0)ret =0;endendendendendend%函数文件 getmaxminpoint.m%功能：取得矩形区域内的点对应的横纵坐标的最大值和最小值function xmin, xmax, ymin, ymax = getmaxminpoint( pa, pb, pc, pd )xs=pa(1), pb(1), pc(1), pd(1);ys=pa(2), pb(2), pc(2), pd(2);xmin = min(xs);xmax =
29、 max(xs);ymin = min(ys);ymax = max(ys);end22%脚本文件 main.m%计算并输出最小石方量对应的底面矩形坐标和海拔高度x=0:30:900;y=0:30:1500;long = 800;%平整块矩形区域的长width = 600;%平整块矩形区域的宽X,Y=meshgrid(x,y);Z=xlsread(data,1,b3:az33);%disp(Z);dp=30;%步长dr=5;%分区域用dh=0.2;%海拔高度步长direct=zeros(4,2);%开挖方向，存放底面投影的四个顶点seaheight=0;%开挖海拔高度minmincost=0;
30、%最小费用代价pexi=0;% E 的横坐标while pexi 0S = S + temp*dr*dr;% 挖土量elseT = T - temp*dr*dr;% 填土量endindex = index + 1;endif (S 0S = S + temp*dr*dr;% 挖土量elseT = T - temp*dr*dr;% 填土量endindex = index + 1;endif (S 0S = S + temp*dr*dr;% 挖土量elseT = T - temp*dr*dr;% 填土量endindex = index + 1;endif (S T)cost = 4*T/3;els
31、ecost = T/3 + S;endif (cost mincost) | (mincost = 0)mincost = cost;height = zh;endzh = zh+dh;end% 固定矩形区域 %if (minmincost mincost) | (minmincost = 0)minmincost = mincost;seaheight = height;direct = pa, pb, pc, pd;endpfxj = pfxj + dp;endpeyi = peyi + dp;endpexi = pexi + dp;endfprintf(the final direct);disp(direct);fprintf(the final seaheight);disp(seaheight);25%输出最后的计算结果fprintf(the final direct);disp(direct);fprintf(the final seaheight);disp(seaheight);
