嫦娥三号软着陆轨道设计与控制策略大学生数学建模竞赛B题获奖.doc

精选资料

嫦娥三号软着陆轨道设计与控制策略

摘要
本文首先在变推力发动机使加速度线性变化的条件下，给出了位置速度状态参数以及推力加速度、推力和秒流量的计算模型，建立了软着陆过程的运动方程，并根据质量计算公式得出了速度增量最小时燃耗最小的结论，此时满足着陆轨道最优的条件，并由上述结论求出该轨道下最小速度增量为1749m/s, 近月点到着陆点的月心角为，即确定了近月点的位置，根据近月点、远月点和月心在一条直线上也就确定了远月点的位置，然后根据开普勒定律的速度计算公式，确定了椭圆轨道中近月点和远月点的速度。近地点速度为1.691km/s，方向与月心和近地点连线方向垂直，即速度方向与月心和着陆点连线的夹角为82.2；远地点速度为1.612km/s，方向与月心和远月点的连线方向垂直，且与近月点速度方向相反。
其次，为了确定软着陆过程6个阶段的最优控制策略，本文根据不同阶段的运动特性建立了模型并给出了3种控制方案，并分别进行了误差分析和敏感性分析，分析了各方案中不同因素对方案的影响，得到可靠的误差范围，并进行了优化。
其中，方案一是针对主减速和快速调整段提出的。由于在主减速和快速调整段中，着陆器距离月面相对较高且着陆器走过的月面距离较长，将月球视为平面建立模型会带来较大的偏差。因此，本文将月球视为球体建立了三维动力学模型，表示出着陆器下降速度在坐标系三轴上的分量。然后给定初值进行迭代，从而求得协状态变量或中间变量，最终获得最优控制方案。然后采用蒙特卡洛打靶，假设各误差均符合正态分布，得出了着陆误差分布在1km范围内的结论且在绝大多数情况下着陆的水平速度不大于1m/s。
方案二是针对粗避障和精避障段提出的。为了避开障碍物，本文采用了基于最大类间方差法的故障检测法，通过这种方法，利用MATLAB对距离月面2400m和100m处的数字高程图进行分析，从而确定故障区域和安全区域。由于存在多个满足条件的区域可以保证着陆器安全着陆，本文又采用了基于螺旋搜索的着陆点选择方法，该方法可以在存在多个满足条件的安全着陆区域的情况下，兼顾能量消耗最少的原则，选择距离当前位置较近的区域实施着陆，但是该方法的误差范围较大，可能会对着陆区域的选择造成较大的偏差。
方案三是针对缓速下降阶段提出的。由于在缓速下降段中，着陆器距离月面很近，且着陆器几乎沿竖直方向下降。因此，本文在平面月球二维模型的基础上进行简化，建立了一维动力学模型，得出了该阶段不同过程加速度的计算方法，在考虑着陆的安全性的前提下，设计了缓速下降阶段的控制方案。该控制方法对初始位置和速度偏差的影响不敏感，因此也无法对初始偏差造成的着陆误差进行修正。
最后，对于上述模型进行了模型评价，给出了各模型和控制方案的优缺点。

关键词：软着陆轨道 运动方程 最优控制 三维动力学模型 基于最大类间方差法的故障检测法 基于螺旋搜索的着陆点选择方法 MATLAB 一维动力学模型 误差分析 敏感性分析

一、问题重述
嫦娥三号于2013年12月2日1时30分成功发射，它携带中国第一艘月球车，并实现中国首次月面软着陆。在整个“落月”过程中，嫦娥三号要完全依靠自主导航控制，完成降低高度、确定着陆点、实施软着陆等一系列关键动作。嫦娥三号在高速飞行的情况下，要保证准确地在月球预定区域内实现软着陆，关键问题是着陆轨道与控制策略的设计。其着陆轨道设计的基本要求：着陆准备轨道为近月点15km，远月点100km的椭圆形轨道，并从近月点开始着陆。整个软着陆过程共分为6个阶段，要满足各个阶段在关键点处所处的状态，并在此基础上确定其着陆轨道和各个阶段的最优控制策略，以进行精准控制，从而实现成功着陆。此外还需进行对设计的控制策略做相应的误差和敏感性分析。
二、问题分析

问题一要求确定着陆准备轨道近月点和远月点的位置，以及嫦娥三号相应速度的大小与方向。针对上述问题，我们需要建立模型得出近月点相对于着陆点的位置，然后根据近月点、远月点和月心在同一条直线上确定远月点的位置，再根据开普勒定律计算出速度。
问题二和问题三要求确定嫦娥三号的着陆轨道和6个阶段的最优控制策略并对控制策略做相应的误差分析和敏感性分析。因此，我们需要根据各个阶段的轨道特性和运动参数建立不同的模型，再制定相应阶段的控制方案，实现轨道优化，并进行误差和灵敏度分析，最后再分析各个模型的优缺点。

三、模型假设

1. 假设月球软着陆探测器由轨道器和着陆器组成。
2. 着陆器降落到月面是个单程任务，级它无须返回地球或在此飞起到达某一环月轨道。
3. 着陆地点由着陆器在环月轨道上的变轨点所决定。
4. 不考虑摄动影响且忽略月球自转。
5. 发动机为变推力发动机，且探测器加速度随时间呈线性变化。

四、符号说明

：软着陆加速度
：初始位置坐标
：初始速度
()：当地引力加速度
θ：初始点到着陆点的月心角
：月球半径
：月球引力常数
：单位时间燃料消耗的公斤数
：速度增量

五、模型建立与求解
5.1 模型一：软着陆过程的运动方程模型
5.1.1 模型说明
由于近月点、远月点和月心在一条直线上，因此只要确定了近月点的位置，即可确定远月点的位置。为了确定近月点的位置，本文建立了软着陆过程的运动方程模型，在变推力发动机使加速度线性变化的条件下，给出了位置速度状态参数以及推力加速度、推力和秒流量的计算模型，得出了速度增量最小时可实现最佳着陆的结论，并由此结论解得了此时近月点到着陆点的月心角，确定了近月点的位置，同时确定了远月点的位置。然后根据开普勒定律的速度计算公式，确定了椭圆轨道中近月点和远月点的速度。

5.1.2 模型的建立和求解
5.1.2.1 近月点和远月点位置的确定
在月心惯性坐标系中，采用变推力火箭发动机，使探测器加速度随时间呈线性变化。加速度可表示为：
；
； (1)

式中： 为常系数；t为时间。
软着陆轨道任一时刻的位置和速度分别为：
(2)
因式(1)中加速度为推力加速度与引力加速度之和，由此可得火箭发动机需施加的加速度：
(3)
为推力加速度，是时间的函数，且
(4)
r=；若探测器质量为m,发动机比冲为,则发动机推力： F= (5)
F的方向可由计算的推力加速度求得。单位时间燃料消耗的公斤数：
=- (6)
由式(6)可得质量计算公式为：
(7)
分别为探测器初始和最终的质量，。
综上所述：最小时燃耗最小，此时满足着陆轨道最优的条件，再根据着陆点的位置确定初始点即近月点的位置。因此，将求近月点位置的问题转化为求最小速度增量的问题。
对如图1所示的月球软着陆轨道示意图，设坐标系基准平面位于初始运动平面内，ox轴由月心指向初始点，oy轴沿运动方向垂直于ox轴，oz轴与ox、oy轴构成右手坐标系。着陆点位置的表示如下：
(8)

图1软着陆轨道示意图

式中：为月球半径，=1737km； 为初始点到着陆点的月心角。给定飞行时间区间[,]和着陆区域[,]，最小速度增量计算步骤如下：
1）计算飞行时间=。
2）计算着陆位置。代入(1)-(4)、(8)可求得与对应的速度增量。
3）计算着陆位置，。同样可求得。
4）当时，令
5）计算
6）计算速度增量、。
7）当时，令

8）当时，令

以上计算过程中，初始状态和着陆点的位置和速度分别为：
，
。初始状态为近月点状态。通过上述步骤，解得最小速度增量为1749m/s,对应的。
因此，近月点到着陆点的月心角为，且近月点高度为15km，即确定了近月点的位置。远月点在近月点与月心的连线上，且远月点高度为100 km，即确定了远月点的位置。

5.1.2.2 速度大小与方向的求解
已知月球和轨道常数如下：
近地点高度：h1=15km
远地点高度：h2=100km
月球平均半径：=1737.013km
椭圆轨道长半轴a=(2*+h2+h1)/2=1794.513km
椭圆轨道半焦距c=(h2-h1)/2=42.5km
椭圆轨道半短轴 b= =1794.0096km
偏心率e= =0.02368
远地点半径：r2=+h2=1837.013km
近地点半径：r1=+h1=1752.013km
月球引力常数：=4.8875*1012 m3/s2=4887.5km3/s2
根据开普勒定律[1]，卫星椭圆轨道上的运动速度为：
V=
所以，近地点速度：V2==1.691km/s，方向与月心和近地点连线方向垂直，即速度方向与月心和着陆点连线的夹角为82.2。
远地点速度：V1==1.612km/s，方向与月心和远月点的连线方向垂直，且与近月点速度方向相反。

5.2 最优控制策略的方案设计及误差分析和敏感性分析
5.2.1 嫦娥三号着陆轨道6个阶段的任务和参数分析
1)着陆准备轨道：近月点15km，远月点100km，将在近月点处以抛物线开始动力下降，初始速度1.691km/s。
2)主减速段：距参考月面高度从15到约3km。该段主要任务是减速制动，减小着陆器约1.7km/s的速度至约57m/s。
3)快速调整段：距月面高度从约3km到约2.4km，快速衔接主减速和后续的阶段。
4)粗避障段：距月面高度从约2.4km到约100m，通过光学成像敏感器检测大障碍，确定安全着陆区并避障，最终到达着陆区上方约100m。
5)精避障段：距月面高度从约100m到约30m，根据选择的安全着陆点，着陆器下降到着陆点上方30m处，水平速度接近于0 m/s。
6)缓速下降阶段：距月面高度从约30m到约4m。该段主要任务是消除水平速度，保证着陆器姿态垂直月面。
最后，在距月面约4m处，关闭发动机和推力器，着陆器自由下降到月面。

5.2.2方案一：主减速段和快速调整段的控制方案设计
5.2.2.1 方案说明
由于在主减速和快速调整段中，着陆器距离月面相对较高且着陆器走过的月面距离较长，将月球视为平面建立模型会带来较大的偏差。因此，本文将月球视为球体建立了三维动力学模型，表示出着陆器下降速度在坐标系三轴上的分量。根据文献[2]提出的方法设计了相应的控制方案，通过计算得到该方案中解析形式的推力角，从而控制主减速段和快速调整段轨道的下降方向。

5.2.2.2 模型二：三维动力学模型
首先，建立如下几个坐标系：
（1）参考惯性坐标系。原点O位于月球中心，轴位于环月轨道平面内且指向前进方向， 轴由月心指向初始着陆点，轴与，构成直角坐标系。
（2）下降轨道参考坐标系。原点O位于着陆器质心，轴由月心指向着陆器质心，轴位于当地水平面内且指向着陆器前进方向。
（3）着陆器坐标系。原点O位于着陆器质心，轴在制动推力矢量延长线上，延推力方向为正，，轴分别根据着陆器上仪器设备的安装而定，并与轴构成直角坐标系。
坐标系示意图以及着陆器位置与推力矢量关系如图1所示。（a）给出了各坐标系的示意和着陆器在坐标系中的位置，（b）给出了F在下降轨道参考坐标系中的位置。

图1 坐标系示意图以及着陆器位置与推力矢量关系

其中，α为在平面内的横向月心角；β为下降轨道平面内的纵向月心角；推力F与坐标系之间的2个推力方向角分别为推力方位角Ψ和推力仰角θ，推力方位角绕正轴旋转为正，推力仰角绕负轴旋转为正。
分别用u，v，w表示着陆器下降速度在坐标系三轴上的分量，于是有：w=，u=r，v=rsin。利用球坐标系与直角坐标系的关系可得到如下所示的下降轨道参考坐标系下的三维动力学模型：
(9)
对于(9)式表示的动力学模型，通常是给定初值进行迭代，从而求得协状态变量或中间变量，最终获得最优控制。该方法不利于在探测器上实现自主控制。文献[2]利用当前状态进行推力角控制量的单步优化控制。本文也采用这样的方法，具体计算见文献[2]。这里直接给出2个推力方向角的控制方程：
(10)
式中，下标f表示终端条件，r,u,v表示当前时刻的下降参数；表示当前时刻的径向加速度，为当前时刻的水平推力加速度。

5.2.2.3 方案一的误差分析和敏感性分析
此模型设涉及的误差源主要包括导航设备测量误差和执行机构误差2部分。还包括月球引力摄动和日、地引力摄动等环境干扰引起的误差。其中，测量误差包括地面测轨误差、惯性装置测量误差、多普勒测速雷达误差和推力误差。
本文给出了测量和推力综合误差情况下采用变推力着陆方案的着陆参数和着陆误差分布情况。以下关于误差的分析均采用蒙特卡洛打靶[3]，打靶次数为500次，假设各误差均符合正态分布。
着陆器在月面2个方向上打靶的着陆点散布情况如图2所示。

图2 变推力下的着陆位置误差分布情况
由图2可以看出，着陆误差分布在1km范围内。
着陆时在月面两个方向的速度误差的散布情况如图3所示。

图3 变推力下的着陆速度横向和纵向误差分布柱状图
图3表明在绝大多数情况下着陆的水平速度不大于1m/s，满足要求。

5.2.3方案二：粗避障和精避障段控制方案的设计
5.2.3.1 方案说明
粗避障和精避障段的主要任务是精确避障和下降。为了避开障碍物，本文采用了一种基于最大类间方差法[4]的故障检测方法，该方法可以自动地确定合理的阈值，对图像进行分割，对故障进行检测，避免了人工选择阈值的弊端，实现了自适应的故障检测。通过这种方法，对距离月面2400m和100m处的数字高程图进行分析，从而确定精确的故障区域和安全区域。由于存在多个满足条件的区域可以保证着陆器安全着陆，本文又采用了一种基于螺旋搜索[4]的着陆点选择方法，该方法可以在存在多个满足条件的安全着陆区域的情况下，兼顾能量消耗最少的原则，选择距离当前位置较近的区域实施着陆。

5.2.3.2 模型三：基于最大类间方差法的故障检测法
方差是图像中灰度分布均匀性的一种度量，方差值越大说明两类图像之间的差别越大，类间方差最大代表了分错的概率最小，分割的效果也就最好。阈值分割法的实质是按照某个准则函数求最佳阈值的过程。最大类间方差就是一评价函数为基础的算法。评价函数有很多种，其中类间最大方差和类内最小方差比较适合。由于类间最大方差只需要计算一阶统计数据，而类内最小方差需要计算二阶统计数据，因此类间最大方差更适合做评价函数，此方法就叫做最大类间方差法，其数学原理参见文献[4]。
利用最大类间方差法对距离月面2400m和100m处的数字高程图进行故障检测，检测结果如图4和图5所示(源程序见附录一和附录二)。

(a) 距月面2400m处的数字高程图
(b)

(b)灰度直方图

(c) 自适应阈值分割后的故障阴影图像
图4距月面2400m处的数字高程图故障阴影提取过程

(a) 距月面100m处的数字高程图

(b)灰度直方图

(c)自适应阈值分割后的故障阴影图像
图5距月面100m处的数字高程图故障阴影提取过程

图2和图3中，(a)是原始图像，(b)图是(a)的灰度直方图，反映了(a)的像素分布。通过最大类间方差法，不用人工设置阈值，程序自动确定阈值，然后对图像进行阈值分割即二值化，二值化后的图像结果如(c)所示。(c)中黑色阴影区域是故障区域，像素值为0；白色区域是安全着陆区，像素值为1。

5.2.3.3 模型四：基于螺旋搜索的着陆点选择方法
首先，确定着陆区域面积大小。一般是一个圆形区域或正方形区域。然后以预定着陆点存在的着陆方形区域开始，以螺旋形由内向外逐步搜索，如图6所示，并按一定的步骤对每个方形区域进行检测判定，如果该区域满足条件，则把该区域标定位安全着陆区，如果不满足条件，则按照规则对下一个区域进行检测判定，每次移动的距离可以适当调整。着陆区判定步骤参见文献[4]。当搜索完整个区域后将所有满足条件的区域作为备选着陆区，计算这个区域与预定着陆点的距离，将距离最近的一个区域作为第一着陆区，其它作为备份，当第一着陆区不再安全时，选择备选区实施着陆。

图6 螺旋搜索的着陆点选择
根据前文得到的二值图像，我们可以测算出图像中某一区域黑色阴影面积的面积即为该区域危险地势的面积，同样可以计算出该区域阴影面积占整个区域的比例，即该区域的危险地势的比例，该比例可以看做探测器能否在该地区安全着陆的指标和约束。将整个预定着陆区域图像进行分割，每个分割区为可能的着陆区，然后计算每个分割区的危险地势比例，进行比例，得出各个区域的安全系数排名，安全系数最高和较高的作为备选的着陆区。结合节省能量的原则，从备选着陆区里选出一个最合适的，其他作为冗余。
假设探测器着陆区域为一正方形，首先将检测到的数字高程图的灰度图像进行分割，分割的每块大小为探测器着陆时所需面积的值。Matlab对着陆图像的分割结果如图7和图8所示(源程序见附录二)。

图7 距月面2400m的障碍分割结果

图8 距月面100m的障碍分割结果

然后用Matlab(源程序见附录二)对每块区域进行阴影区域计算，得出各区域危险地势的比例，如表1和表2所示。
表1 距月面2400米着陆区域各分割区内危险地势比例数据表
编号
黑色像素点数
比例（%）
编号
黑色像素点数
比例
1
138983
65.6819
14
171898
81.2372
2
168007
79.3983
15
87953
41.5656
3
185367
87.6025
16
38479
18.1847
4
47376
22.3894
17
322
0.1521
5
0
0
18
110815
52.3700
6
73091
34.5420
19
124865
59.0099
7
73965
34.9551
20
12934
6.1124
8
97509
46.0817
21
0
0
9
41912
19.8071
22
128
0.0604
10
4460
2.1077
23
70788
33.4536
11
177
0.0836
24
15294
7.2277
12
23043
10.8898
25
18
0.0085
13
32997
15.5940

表2 距月面100米着陆区域各分割区内危险地势比例数据表
编号
黑色像素点数
比例（%）
编号
黑色像素点数
比例
1
1558
3.8950
14
15354
38.3850
2
3
0.0075
15
850
2.1250
3
22
0.0550
16
21736
54.3400
4
26
0.0650
17
29891
74.7275
5
14936
37.3400
18
11762
29.4050
6
7437
18.5925
19
1654
4.1350
7
20286
50.7150
20
416
1.0400
8
4179
10.4475
21
0
0
9
59
0.1475
22
488
1.2200
10
10439
26.0975
23
25
0.0625
11
26235
65.5875
24
181
0.4525
12
31130
77.8250
25
137
0.3425
13
25440
63.6000

兼顾能量消耗最少的原则，在满足着陆安全性的要求下，着陆区的选择应当尽可能接近当前着陆器位置。因此，距月面2400m时将区域11作为首选着陆点，区域5和区域17为备选着陆区；距月面100 m时将区域9作为首选着陆点，区域15为备选着陆区。

5.2.3.4 方案二的误差分析和敏感性分析
以距月面100m处着陆区域各分割区内危险地势比例数据表为例，假设探测器可以容忍的危险地势阈值极限值为3%，则区域2、3、4、9、15、20、21、23、24、25都是适合着陆的安全区。此时，虽然有效性较高，但是可靠性降低了。为了进一步提高着陆安全可靠性，设置阈值为1%，则区域2,3,4,9,21,23,24,25符合要求。再设置阈值为0.2%，则区域2,3,4,9,21符合要求。
因此，该方法的误差范围较大，要选择最佳着陆区域，应合理设置探测器可以容忍的危险地势阈值极限值。否则，会对着陆区域的选择造成较大的偏差，或者造成过多的燃料消耗，影响着陆的安全性。

5.2.4 方案三：缓速下降阶段的控制方案设计
5.2.4.1 方案说明
软着陆的垂直降落阶段，轨道高度在100m以下，探测器的垂直速度衰减到几米每秒的量级，探测器在姿控系统的控制下作盘旋运动，姿控系统的推力方向产生非常复杂的变化，这一阶段对探测器的轨道实时确定要求不高，主要考虑探测器的纵向速度对探测器着陆缓冲装置的影响，对姿控系统提出较高要求，要保证探测器姿态平稳，不会倾倒。探测器依靠自身携带的地形识别系统和测距测速敏感器，对降落过程实时制导。
由于在缓速下降段中，着陆器距离月面很近，且着陆器几乎沿竖直方向下降。因此，本文在平面月球二维模型[5]的基础上进行简化，建立了一维动力学模型，得出了该阶段不同过程加速度的计算方法，在考虑着陆的安全性的前提下，设计了缓速下降阶段的控制方案。该控制方法对初始位置和速度偏差的影响不敏感，因此也无法对初始偏差造成的着陆误差进行修正。

5.2.4.2模型五：一维垂直动力学模型
首先将月球视为平面建立与月球平面直角坐标系，如图9所示:

图9 月球平面直角坐标系

图中，原点O为下降轨道上制动发动机点火在月球表面上的阴影，为下降轨道参考系纵向平面，着陆器的下降轨迹位于此平面内。图8表示的是符合重力转弯软着陆的情况，即反推力F的方向与下降速度方向相反。对于这样的情况，沿坐标轴方向有如下的动力学方程：
(11)
式中，m为飞行器质量，在短时间内可是视为常量：为月球表面的重力加速度，始终垂直于月球表面且为常值；为飞行路径角，即为下降速度矢量与轴的夹角，从轴开始逆时针量起为正，为下降速度的模，。
由于缓速下降阶段着陆器几乎沿垂直方向下降，可在上述平面月球二维模型的基础上简化为一维垂直动力学模型，即将其中的飞行路径角设为90度，因此，(11)式可简化为：
(12)
其中，u为推力F的开关控制量。
一维垂直下降过程如图10所示：

图10 一维垂直下降过程示意图

5.2.4.3 方案设计以及方案的误差分析和灵敏度分析
缓速下降和自由落体阶段，着陆器依次经过悬停、匀加速、匀减速和关机降落几个过程。几个过程均符合牛顿定律，易得开关切换高度：
= (13)
其中，合加速度，。
考虑到着陆的安全性，在着陆段初始要进行段时间的悬停以对着陆区域进行成像勘察，且由于着陆时间很短，因此应保证着陆器平缓下降，尽量避免受制动发动机的开关冲击。
该控制方法对初始位置和速度偏差的影响不敏感，因此该方法也无法对初始偏差造成的着陆误差进行修正。

六、模型评价
模型一通过计算最小速度增量确定了初始点的位置，该计算方法简单，易于工程化，适用于变推力火箭发动机在月球上实现定点软着陆的任务，但不能具体求解出初始点的位置参数，只能求得初始点和着陆点的相对位置。
模型二将月球视为球体建立了三维动力学模型，并针对上述模型，给定初值进行迭代，从而求得中间变量，最终获得最优控制。但是该方法不利于在探测器上实现自主控制。
模型三和模型四分别进行了故障检测和着陆点的选择，初步避开了危险区域，得到了安全着陆的最佳区域，缺点是该方案有效性和可靠性不可兼顾，在设置探测器可以容忍的危险地势阈值极限值时可能会使安全性降低。
模型五在平面月球二维模型的基础上进行简化，设计了一维垂直动力学模型，并提出了缓速下降阶段控制方案的设计。优点是充分考虑了着陆的安全性，缺点是模型对初始位置和速度偏差的影响不敏感，因此该方法也无法对初始偏差造成的着陆误差进行修正。

七、参考文献
[1] 朱立东 ，《卫星通信导论》，出版社：电子工业出版社，2010年。
[2] Ueno S, Yamaguchi. 3-dimensional near-minimum fuel guidance law of a lunar landing module,AIAA Guidence, Navigation, and Control Conference and Exhibit, Portland, OR, 1999:AIAA-99-3983。
[3] 王鹏基，张熇，曲广吉. 月球软着陆飞行动力学和制导控制建模与仿真. 2009,39(3):521-527。
[4] 董士波 火星表面障碍识别与规避方法研究 哈尔滨工业大学硕士学位论文，2013。
[5] 褚桂柏，张熇. 月球探测器技术. 北京：中国科学技术出版社，2007. 152-154

附录一：运用Matlab编写的图像二值化处理源程序
clear all;clc;
I = imread('len_100.tif');
figure
imshow(I)
title('原始图像')
figure
imhist(I)
title('原始图像的直方图')
level = graythresh(I)
BW = im2bw(I,level);
figure
imshow(BW)
title('根据动态阈值做二值化')

附录二：运用Matlab编写的图像分割与统计源程序
clear all;
I=imread('len_100.tif');
%figure,imshow(a)

%将图像二值化
%I=a;
level = graythresh(I);
BW = im2bw(I,level);
figure
%subplot(231);
imshow(BW)
title('根据动态阈值做二值化')
a=I;

count=imhist(a);
[m,n]=size(a);
N=m*n;
L=256;
count=count/N;%%每一个像素的分布概率
count

for i=1:L
if count(i)~=0
st=i-1;
break;
end
end
st
for i=L:-1:1
if count(i)~=0
nd=i-1;
break;
end
end
nd
f=count(st+1:nd+1); %f是每个灰度出现的概率
size(f)
E=[];
for Th=st:nd-1 %%%设定初始分割阈值为Th
av1=0;
av2=0;
Pth=sum(count(1:Th+1));
%%%第一类的平均相对熵为
for i=0:Th
av1=av1-count(i+1)/Pth*log(count(i+1)/Pth+0.00001);
end
%%%第二类的平均相对熵为
for i=Th+1:L-1
av2=av2-count(i+1)/(1-Pth)*log(count(i+1)/(1-Pth)+0.00001);
end
E(Th-st+1)=av1+av2;
end
position=find(E==(max(E)));
th=st+position-1

for i=1:m
for j=1:n
if a(i,j)>th
a(i,j)=255;
else
a(i,j)=0;
end
end
end

%figure,imshow(a);
level = graythresh(a);
quan = im2bw(a,level);
figure
%subplot(232);
imshow(quan)
title('熵最大法分割图像')

lst=~quan;
fna1=(BW+lst);

% lst=(fna+quan);
figure
fna2=~fna1;
%subplot(233);
imshow(fna2);

%将图像进行分割
%figure
divd=142;
%fen=fna2(1:divd,divd*1+1:divd*(1+1));
%imshow(fen);
figure;
ci=1;
k=0;
tongji1=zeros(7,7);
tongji2=zeros(7,7);

for hang=0:6
for lie=0:6
fen=fna2((divd*hang+1):(divd*(hang+1)),(divd*lie+1):(divd*(lie+1)));
for i=1:142
for j=1:142
if fen(i,j)==0
k=k+1; %统计黑色像素点的个数
end
end
end
bili=(k/(142*142))*100;
tongji1(hang+1,lie+1)=k;
tongji2(hang+1,lie+1)=bili;
k=0;
subplot(7,7,ci);
ci=ci+1;
%figure;
imshow(fen);
end
end

THANKS !!!

致力为企业和个人提供合同协议，策划案计划书，学习课件等等
打造全网一站式需求

欢迎您的下载，资料仅供参考，如有侵权联系删 除！

可修改编辑