2023A数学建模优秀论文

2023高教社杯全国大学生数学建模竞赛承 诺 书（。（（包括指导教师（包括网上查到的资料，必须按照规定的参考文献的表（。我们参赛选择的题号是（从A///D中选择项填写： A 我们的报名参赛队号为（8位数字组成的编号：1009072 所属学校（请填写完整的全名：东南大学 参赛队员(打印并签名)：1.吉张鹤轩 杨升 陈同广 指导教师或指导教师组负责人 (打印并签名)：（论文纸质版与电子版中的以上信息必须一致，只是电子版中无需签名。以）日期：2023年09月15日赛区评阅编号（由赛区组委会评阅前进行编号）：2023高教社杯全国大学生数学建模竞赛编号专用页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：阅人分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：嫦娥三号软着陆轨道设计与控制策略摘要1055.39kg26.71kgK均值聚25,1008.97（88,2068k8.9kg关键词：线规模型 序化传法 K均聚类 空线回归 体模型一．问题的提出背景介绍根据计划，嫦娥三号将在北京时间12月14号在月球表面实施软着陆。嫦娥三号如何实现软着陆以及能否成功成为外界关注焦点。目前，全球仅有美国、前苏联成功实施了13次无人月球表面软着陆。1210“”“”地“100问题重述1500N19.51W，44.12N-2641m。15km，100km6根据上述的基本要求，建立数学模型解决下面的问题：（1）确定着陆准备轨道近月点和远月点的位置，以及嫦娥三号相应速度的大小与方向。（2）确定嫦娥三号的着陆轨道和在6个阶段的最优控制策略。（3）对设计的着陆轨道和控制策略做相应的误差分析和敏感性分析。二．问题的分析问题一NASA1976[1]问题二由问题一已经得出近月点，即开始降落点位置，也知道每一阶段的状态，因此，降落轨SQP00最后根据各个阶段的最优方案，模拟出嫦娥三号着陆轨道即可。问题三三．模型假设由于月球扁率很小，可认为月球为球体，半径以平均半径为准，并且引力场分布均匀；由于侧面姿态调整喷射装置对燃料影响很小，为简化模型，认为飞行器变换姿态的由于着陆时间较短，所以诸如月球引力非球项、日月引力摄动等影响因素均可忽略由于预定着陆点海拔-2641m，因此在着陆过程中所使用的高度均不应是海拔高度，四．符号说明符号符号意义G引力常量M月球质量m嫦娥三号质量比冲gm月球表面重力加速度F推力J燃耗ˆ平坦度评价指标ˆ综合评价指标ˆ平均坡度评价指标五．模型建立与求解问题一模型建立Hf=100kmR=1737.013km5.1.1.1.图5.1.1.1近月轨道示意图（为了方便示意，本图不符合比例）a+c=Hf+R a-c=Hc+R联系万有引力定律与牛顿第二定律，可以列出嫦娥三号在近月点和远月点的运动学方程：Mm v2G0m f(ac)2 0fMm v2G0mcc(ac)2 0c0其中M72gG=21mt0上飞行器的质量，f与c为远月点和近月点的曲率半径。由椭圆的几何性质，在长轴两端点处的曲率半径分别为：c

22b=ac1eaa2f联立以上各式，带入参数，得：

=ac1eb近月点速度大小/s；远月点速度大小vf/s。下面开始研究主减速过程，该过程由近月点开始，切向速度v0vc，径向速度vr等于0，由切线正切制导率及相关报道可知，该阶段推力始终保持最大推力，且始终调整使推力与速度方向相反，即F7500N，Fv，最终应大致达到预定着陆点目标上空，且竖直方向速度为57m/s，高度共下降12000m。由于有端点限制条件，可以通过物理知识建立微分方程组，即采用微分方程模型进行分析求解。v2Fmxd r这样，问题就转化到了平面直角坐标系下，并且可以对速度和受力进行正交分解，如图。图5.1.1.2非惯性系下力学分析图示由此可得微分方程模型：dvy GM v2xasindt r2 rdvxacosvxvydt rdrvdt ym0mtm0mtvytanvxrβaFm=2940N/。至此，本问题已经完全抽象为数学语言，并且有初始参数：vy0，vx/s，r，=012000m模型求解NP猜测微元化后的差分方程依然可以有运动学规律得出：v v

(a

cos

vxnvyn)txn1

xn n nnrGM v2rvyn1

vyn(2rnr

ansinn

xn)tanm m

Fmnmtn1 nyn1 2xn1xnvxnt

ancosntyn1

ynvyn

t1a2

nsinn

t2其中为时微小减速大几秒为确描运曲一般取0.1s。带F，，0=0，0，0，vx0/s，vy00即可通过迭代法，求出在y=3000m时的x的值，该值即为飞行器转过的角度所对应的3000mθ于是我们通过改变θ的大小，试图输出大量运动轨迹簇，这一过程由计算机模拟生成。通过计算机仿真模拟，得到非惯性参考系下的运动轨迹簇。其中可以发现，在θ取时，高度3000m处竖直方向速度为，在误差允许的范围内可以接受，由此得出符合题意的主减速阶段飞行轨迹图，如图5.1.1.4.图5.1.1.4主减速阶段飞行轨迹图（非惯性参考系）图5.1.1.5水平与竖直速度变示意图 图夹角变化示意图5.1.1.5.385.21km。。由月球的球对称性，远月点纬度南北对调，经度东西对调，角度互补。1692.46m/s；远月100km1612.15m/s。问题二由于燃料损耗直接与推力有关，即Fm，其中=2940m/s为燃料比冲，也就有Fdtmdt另外，由动量定理mdvFdt，在近似处理下，也可将燃耗与速度变化联系在一起。过程一模型建立与求解5.1,3,也可r，极角为，角速度为，质量为m；又设v为嫦娥三号沿rF在IP；设月。图5.2.1.1二体模型力学分析示意图首先，从运动学角度进行分析。由动力学基本方程可得drtdv，所以drv①dt同理可得角变的关式为d②dt图5.2.1.2受力分析示意图再从动力学角度进行分析。如图，对嫦娥三号进行受力分析可知。F

GMmm= 月= 以及减发机径上的力Fsinm引 r2 r2由此，根据牛顿第二定律可以得到：m mv2③r2 rFsinmdv ④dt联立③，④可得联立③，④可得mFsinmv dv m ，用式v r化该可：r2 r dtdvFsinr2 ⑤dt m r2需要计入大小为第二定律可以得到Fcos2vmdvy。同样，利用vdt

r化简可以得到：F*cos2*v*dm ⑥.dt r由于飞船在运行过程利用燃料反冲制动，因此飞船得而质量是不断变化的。考虑比冲的定义为“火箭发动机单位质量推进剂产生的冲量”可得：dmF ⑦。dt ISP综上所述，联立式①，②，⑤，⑥，⑦即可得到完整的描述嫦娥三号与月球这一二体运动模型的方程组：drvdtdvFsinr2dt m r2ddtFcos2vdmdt rdmFdt ISP对于嫦娥三号燃料使用指标的衡量，我们可以反映到其动量上去，即将动量作为衡量燃tFtI料消耗的性能指标。由质量变化的推导，我们可以得到JItSP

dt。在轨道优化过程中,由于各状态变量的量级相差较大,寻优过程中可能会导致有效[2]

,mrefm0.rrefrrrefrref

mv2 m则：r,v,v

,I I

,F,F

refref

,m,rref

vref

SP

Fref

ref

rref

mrefr3refr3ref

,t

ttref

,tref

refvref

,.那么，嫦娥三号的动力学方程可改写为：drvdtdvFsin1r2dt m r2ddtFcosdmdt rdmFdt ISP又由第一问可以得到飞行器的初始条件和终端条件分别为:vr00

v0v0

1692.46m/s

rffv 57mfr01749.372km注：此处的终端条件中，径向速度和切向速度为理想情况下的结果，实际只能得到最优化的值。使用归一化条件可以得到初始条件和终端条件分别变为：vr00

vrf0v0

1692.46m/svref

vf

57vref

m/s0r1749.372kmrref00m/sFtF FtIJIt

dt= tfSP ISPFISP

都是定值，所以即要求时间最短的方案。[3]atiii0a0到a6七个变量后，根据所建立的二体模型方程即可确定a0a6七个变量，问题的约束条件为归一化后的f动力学方程以及初值条件和终端条件。由此需要求出燃料最优化方案，即性能指标J=FtfISP的最小值。对以上转化出来的问题，我们利用浮点数编码的遗传算法（FGA）进行求解。其步骤如下：72000N通过对约束条件的分析，可以选择交叉概率为0.6，变异概率为0.05。通过计算机仿真模拟，由遗传算法规划出的角度随时间控制函数参数向量的取值为：0

a67.5810

7.01105

4.30106

8.12108

6.031010

1.601012

1.401015运动轨迹见图5.2.1.2.图5.2.1.3运动轨迹图图5.2.1.4速度对比图 图5.2.1.5角变化图当下降12000m时速度v=57.12m/s，水平方向与垂直方向分别为vx48.63m/s，vy29.96m/s，符合题意。413.2sF恒为7500NJFt1055.39kg，剩余质量1344.61kg，水平位移为385.21km。1 过程二模型建立与求解向的速度减到0m/sF（的mm/s。图5.2.1.6受力分析示意图由此，首先进行运动学分析。在水平方向上，速度最终减小为有0m/s，所以有dvxa ①dt x在竖直方向上，由加速度的定义，有这样的关系：dvya ②dt y然后，在进行动力学分析。在水平方向上，根据牛顿第二定律可以得到：Fcosmaxgm

③2，其值由黄金代换公式GMgmR可2以确定。那么同理可以得到：

Fsin④至此，受力分析完毕。而另一方面，考虑到减速过程中燃料的消耗对质量产生影响，有：dmF ⑤dt ISP综合①，②，③，④，⑤式可以得到描述该阶段运行的方程：Fcosmdvxdtmgm

Fsinmdvy ⑥dtdmFdt ISP由主减速阶段的终值以及快速调整阶段的要求可以得到该方程初始条件和终端条件分别为：vx048.63m/s

vxf

0m/svy029.96m/s ⑦m01344.61kg

f90 ⑧oof90oo而力F的取值范围为1500NF7500N。F可以表示成一个多项式的形式：6t6Fati。另外，在快速调整阶段，燃料的性能指标同样为J

Fdt。ii0

tSPI化的参数为a0到a6以及JI同样使用浮点数编码的遗传算法来解决该非线性规划问题，由于具体步骤与主减速阶段基本一致，故此处不再赘述。在这样的算法下，求得的最优解为：vx0m/svy29.96m/s90o90oF3928.115Nt20.027sJ226.71kg过程三模型建立与求解2.4km100m100m2.4km1m/1m。3。图5.2.3.1高程数量分布曲线图95m120m80m80m120m80m120mK5时，可以较为清楚的区分出陨石坑与平地，如图。图5.2.3.2聚落色阶图。1ˆ2300×230025460×4608ˆS（以格数衡0格。ˆSˆˆS5，即ˆ5ˆSˆH3F2F2F2F2

的约束下完全独立，为了进一步简x y m则速度vxvt5.2.3.3知道，速度曲线2要比速度曲线1燃耗小，由此，可以建立这样的局部优化策略。图5.2.3.3速度曲线对比图通过该阶段的始末速度限制可知，初始水平速度为0，竖直速度vy0由上一阶段得出；段末悬停，速度为0。由此可见，要节省水平方向的燃耗，应最快的增大vx，接着匀速运动，最后对应最快的减小vx至0，从而完成该偏移过程。

F2F

F2=ma，同时，竖直方向匀减速yv2 HyFXvtxt3yy x a vy

J2mvxJv且联立前面各式可以发现vX3 3 x xH v2x3vvyx

xXa即vxXX500mXvx设偏移量评估指标为X为ABepXABpA200，B2.4，p0.0027X5.2.3.4.图5.2.3.4偏移量评价函数曲线图最终制定综合评估指标ˆˆˆ，作为一个方格对于降落点优先级的考量，指标评估对栅格化地图遍历并计算后得到如图5.2.3.5的结果。图5.2.3.5综合评估地图5.2.3.3127510X为15.26平方向最大速度vx1.87m/s，用时t3102.51sJ386.97kg。过程四模型建立与求解100m30m100m处，100m30m0m/s。0.1m/0.1m。100m×100m对于地图平坦度的分析所要考虑的因素不仅有安全半径大小ˆS，而且还要同时考虑到安全格内的平均坡度本次栅格化考虑到像素为0.1m，而飞行器本身体积较大，为了综合考虑，应将单元格大小取为5m×5m，即50×50.聚类分析结果如图5.2.4.1.图5.2.4.1聚类色阶图3最大安全半径分析计算，结果如图5.2.4.2.图5.2.4.2最大安全半径色阶图下面重点探讨如何得到栅格内的平均坡度。5.2.4.3.图5.2.4.3平均坡度色阶图下面制定坡度评价指标ˆ0101坡评指即ˆ1ˆ,将其为化数平度评指相耦maxmin合即可完成非线性加权，最终所得到的综合评价指标ˆˆˆ作为着陆点优先级的评估指5.2.4.4.图5.2.4.4综合评价指标色阶图5.2.4.2中（88m，56m）X38.47m。v2vtxXxx axvyt

2a

2vy HvygmFa2a2x4xym4目标函数为J42mvxmvymgmt要尽可能小。在该动态规划下可以得出最佳数值分配为：竖直缓速下降所用时间t25.65s，水平最大速度为vx1.44m/s，燃耗J420.68kg，段末质量m41210.25kg。过程五模型建立与求解30m4m4m4m0来说可以忽略不计，则由动量定理，设此阶段初始状态竖直速度分量为v5，质量为m5，则Fdtm5v5m5gmtmaFmaxgmFmaxH5

H4mgm

v2 5 5 2aGM1.63m/s2为月面重力加速度；由上一阶段R2得出v5FJ58.09kg，t510.42s。着陆轨道阶段最优控制策略燃耗主减速阶段5.2.1着陆轨道阶段最优控制策略燃耗主减速阶段5.2.11055.39kg快速调整阶段之后迅速变小，直到水平速度减到0.26.71kg粗避障阶段0，段末悬空。86.97kg精避障阶段水平方向制动方案与粗避障阶段方案基本一致；竖直方向0.20.68kg缓速下降阶段该阶段先匀速下降，在最后一段在全力减速至0.8.09kg问题三误差分析（1）27.3d500s0.0763°，即2313.14m（2）快速调整阶段水平位移大约600m，对应月面弧度为0.0198°，因此实际近月点纬（3）7500N1500N1s2000kgtvF0t可求出末速度v2.25m/s0.5s0.5s2m/s4000N（5%（4）（（01不计。1000kg10s05.51kgt，Fmgm其中F与m均为t函数解得mt1000t5%敏感性分析N1N甚至10N，1.3%。71%六.模型优缺点问题一的模型评价问题二的模型评价七.参考文献2023.09.ChuangCH,GoodsonTD,HansonJ.Fuel-optimal,low-and-medium-thrustorbittransfersinlargenumbersofburns[R].AIAA9423650,1994,158-166.NotonM.[M].1979...（）2003年438期.八.附录%————————————————————————————————————————%calc_proc.m计算第一问%————————————————————————————————————————T=0.1;t=0;x=0;y=0;x2=0;y2=0;Q=6.5044;sita=0+Q/180*pi;vx=1692;vy=0;m0=2400;m=2400;N=7500/m0;G=3844.6/m0;r=1749372;H=15000;ax=N-vy*vx/r;ay=G-vx^2/r;holdonxlabel'切向距离/m'ylabel'高度/m'history=[];i=0;while(y>-12000&&m>1000)x2=x+vx*T-0.5*ax*T;y2=y-vy*T-0.5*ay*T;vx=vx-ax*T;vy=vy+ay*T;r2=1749372+y2;G=3844.6/m0/(r^2)*(r2^2);m=m-2.55*T;t=t+T;N=7500/m;if(vx==0)ifvy<=0sita=pi/2;elseend

sita=-pi/2;elseend

sita=atan(vy/vx);if(sita)<=0sita=Q/180*pi;else

if(sita+Q/180*pi)<pi/2sita=sita+Q/180*pi;end

elseend

sita=pi/2;end

ax=N*cos(sita)-vy*vx/r2;ay=G-N*sin(sita)-vx^2/r2;plot([x,x2],[H+y,H+y2],'black','LineWidth',2);x=x2;y=y2;r=r2;history=[history;(i-1)*Tyvxvysita];i=i+1;%————————————————————————————————————————%calc_proc2.m模拟第二问最优解情况%————————————————————————————————————————T=0.1;t=0;x=0;y=0;x2=0;y2=0;Q=6.4;sita=fit(7);vx=1692;vy=0;m0=2400;m=2400;N=7500/m0;G=3844.6/m0;r=1749372;H=15000;ax=N-vy*vx/r;ay=G-vx^2/r;holdonxlabel'切向距离/m'ylabel'高度/m'history=[];i=0;while(y>-12000&&m>1000)x2=x+vx*T-0.5*ax*T;y2=y-vy*T-0.5*ay*T;vx=vx-ax*T;vy=vy+ay*T;r2=1749372+y2;G=3844.6/m0/(r^2)*(r2^2);m=m-2.55*T;t=t+T;N=7500/m;sita=fit(1)*t.^6+fit(2)*t.^5+fit(3)*t.^4+fit(4)*t.^3+fit(5)*t.^2+fit(6)*t+fit(7);ax=N*cos(sita)-vy*vx/r2;ay=G-N*sin(sita)-vx^2/r2;plot([x,x2],[H+y,H+y2],'black','LineWidth',2);x=x2;y=y2;r=r2;history=[history;(i-1)*Tyvxvysita];i=i+1;endif(m<=500)differ=200000;elseend

differ=t;if(abs(vx)>10)differ=200000;endif(abs(vy-57)>1)differ=300000;end%————————————————————————————————————————%calc1.m函数，用以使用遗传算法求解%————————————————————————————————————————function[differ]=calc1(Q)T=0.1;t=0;x=0;y=0;x2=0;y2=0;sita=Q*t;vx=1692;vy=0;m0=2400;m=2400;N=7500/m0;G=3844.6/m0;r=1749372;H=15000;ax=N-vy*vx/r;ay=G-vx^2/r;i=0;while(y>-12000&&m>1000)x2=x+vx*T-0.5*ax*T;y2=y-vy*T-0.5*ay*T;vx=vx-ax*T;vy=vy+ay*T;r2=1749372+y2;G=3844.6/m0/(r^2)*(r2^2);m=m-2.55*T;t=t+T;N=7500/m;if(vx==0)ifvy<=0sita=pi/2;else

elseend

sita=-pi/2;sita=atan(vy/vx);endif(sita)<=0sita=Q/180*pi;else

if(sita+Q/180*pi)<pi/2sita=sita+Q/180*pi;end

elseend

sita=pi/2;sita=Q*t;ax=N*cos(sita)-vy*vx/r2;ay=G-N*sin(sita)-vx^2/r2;x=x2;y=y2;r=r2;i=i+1;endif(m<=500)differ=200000;enddiffer=abs(sqrt(vy^2+vx^2)-57);%————————————————————————————————————————%calc2.m函数，用以使用遗传算法求解%————————————————————————————————————————functiondiffer=calc2(fit)T=0.1;t=0;x=0;y=0;x2=0;y2=0;Q=6.4;sita=fit(7);vx=1692;vy=0;m0=2400;m=2400;N_altered=7500;N=N_altered/m0;G=3844.6/m0;r=1749372;H=15000;ax=N-vy*vx/r;ay=G-vx^2/r;i=0;while(y>-12000&&m>1000)x2=x+vx*T-0.5*ax*T;y2=y-vy*T-0.5*ay*T;vx=vx-ax*T;vy=vy+ay*T;r2=1749372+y2;G=3844.6/m0/(r^2)*(r2^2);m=m-2.55/7500*N_altered*T;t=t+T;N=N_altered/m;sita=fit(1)*t.^6+fit(2)*t.^5+fit(3)*t.^4+fit(4)*t.^3+fit(5)*t.^2+fit(6)*t+fit(7);ax=N*cos(sita)-vy*vx/r2;end

ay=G-N*sin(sita)-vx^2/r2;x=x2;y=y2;r=r2;i=i+1;if(m<=500)differ=200000;elseend

differ=m0-m;if(abs(sqrt(vx^2+vy^2)-57)>1)differ=300000;end%————————————————————————————————————————%deal1.m函数，用以处理2400m处的俯拍相片%————————————————————————————————————————cd('./');A=imread('附件3距2400m处的数字高程图.tif');N=460;n=2300/N;data=zeros(N,N);[height,length]=size(data);fori=1:1:heightforj=1:1:lengthdata(i,j)=sum(sum(A(round(n*(i-1)+1):round(n*i)-(round(n*i)>2300),round(n*(j-1)+1):round(n*j)-(round(n*j)>2300))))/round(n*n);endenddata=uint8(round(data));%imshow(data);%第二个阶段，聚落分析平面data3=[];fori=1:1:Nforj=1:1:Ndata3=[data3;data(i,j)];endend[Idx,C,sumD,D]=kmeans(double(data3),5);C(1,3)=sum(C(:,1)<C(1,1));C(1,2)=sum(Idx==1);C(2,3)=sum(C(:,1)<C(2,1));C(2,2)=sum(Idx==2);C(3,3)=sum(C(:,1)<C(3,1));C(3,2)=sum(Idx==3);C(4,3)=sum(C(:,1)<C(4,1));C(4,2)=sum(Idx==4);C(5,3)=sum(C(:,1)<C(5,1));C(5,2)=sum(Idx==5);GND=sum(C(:,3).*(C(:,2)==max(C(:,2))));data4=zeros(N,N);fori=1:1:Nforj=1:1:Ndata4(i,j)=C(Idx((i-1)*N+j),3);endend%surface(data4,'EdgeColor','none');%第三个阶段选择合适区域data2=1-(data4==GND);score=zeros(N,N);fori=1:1:heightforj=1:1:lengthif(data2(i,j)==0)tmp=1;while(tmp<=i&&tmp<=N-i+1&&tmp<=j&&tmp<=N-j+1&&sum(sum(data2(i-tmp+1:i+tmp-1,j-tmp+1:j+tmp-1)))==0)tmp=tmp+1;end

end

endtmp=tmp-1;score(i,j)=tmp+50-exp(0.0027*sqrt((i-N/2)^2+(j-N/2)^2)*n);%surface(score,'EdgeColor','none');%————————————————————————————————————————%deal2.m函数，用以处理100m处的俯拍相片%————————————————————————————————————————cd('./');Aimread('附件4100m.tif');N=100;n=1000/N;data=zeros(N,N);[height,length]=size(data);fori=1:1:heightforj=1:1:lengthdata(i,j)=sum(sum(A(round(n*(i-1)+1):round(n*i)-(round(n*i)>2300),round(n*(j-1)+1):round(n*j)-(round(n*j)>2300))))/round(n*n);endenddata=uint8(round(data));%imshow(data);%第二阶段%分析斜度yuzumi=double(zeros(N,N));N2=50;n2=1000/N2;ana_width=n2+10;fori=1:1:N2forj=1:1:N2X=[];z=[];forx=max((i-1)*n2+1-(ana_width-n2),1):1:min(i*n2+(ana_width-n2),1000)fory=max((j-1)*n2+1-(ana_width-n2),1):1:min(j*n2+(ana_width-n2),1000)X=[X;1,x,y];z=[z;double(A(x,y))];endendend

end

b=regress(z,X);arc=[b(2),b(3),1];sita=acos(dot(arc,[001])/sqrt(dot(arc,arc)));yuzumi((i-1)*(N/N2)+1:i*(N/N2),(j-1)*(N/N2)+1:j*(N/N2))=sita;heion=(max(max(yuzumi))-yuzumi)/max(max(yuzumi));%surface(heion,'EdgeColor','none');%聚类分析data3=[];fori=1:1:Nforj=1:1:Ndata3=[data3;data(i,j)];endend[Idx,C,sumD,D]=kmeans(double(data3),4);C(1,3)=sum(C(:,1)<C(1,1));C(1,2)=sum(Idx==1);C(2,3)=sum(C(:,1)<C(2,1));C(2,2)=sum(Idx==2);C(3,3)=sum(C(:,1)<C(3,1));C(3,2)=sum(Idx==3);C(4,3)=sum(C(:,1)<C(4,1));C(4,2)=sum(Idx==4);GND=sum(C(:,3).*(C(:,2)==max(C(:,2))));data4=zeros(N,N);fori=1:1:Nforj=1:1:Ndata4(i,j)=C(Idx((i-1)*N+j),3);endend%surface(data4,'EdgeColor','none');%第三个阶段选择合适区域data2=1-(data4==GND);score=zeros(N,N);fori=1:1:heightforj=1:1:lengthif(data2(i,j)==0)tmp=1;while(tmp<=i&&tmp<=N-i+1&&tmp<=j&&tmp<=N-j+1&&sum(sum(data2(i-tmp+1:i+tmp-1,j-tmp+1:j+tmp-1)))==0)tmp=tmp+1;end

end

endtmp=tmp-1;score(i,j)=tmp;%.*heion(i,j);%surface(score,'EdgeColor','none');%————————————————————————————————————————%genpics.m各种画图指令的汇总%————————————————————————————————————————cd('./');clear;load('2400m处分析数据.mat');%聚落分析figure;cla;surface(data4,'EdgeColor','none');colorbar;saveas(gcf,'2400m处聚落分析.png');%合适区域cla;surf(score,'EdgeColor','none');colorbar;saveas(gcf,'2400m处落点评价.png');clear;load('100m处分析数据.mat');%聚落分析cla;surface(data4,'EdgeColor','none');colorbar;saveas(gcf,'100m处聚落分析.png');%合适区域cla;surf(score,'EdgeColor','none');colorbar;saveas(gcf,'100m处落点评价.png');clear;cla;calc_proc2;saveas(gcf,'第一阶段降落y-x轨迹.pn