数学建模卫星和飞船的跟踪测控-毕业论文.doc

卫星和飞船的跟踪测控高教社杯全国大学生数学建模竞赛承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B/C/D中选择一项填写）： 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）： 参赛队员 (打印并签名) ：1. 2. 3. 指导教师或指导教师组负责人 (打印并签名)： 日期： 年 月 日赛区评阅编号（由赛区组委会评阅前进行编号）：高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：1.问题的提出根据问题背景和所提供的信息，本文致力于解决以下三个问题：(1)在假设所有测控站都与卫星或飞船的运行轨道共面的情况下，计算至少应该建立多少个测控站才能对其进行全程跟踪测控；(2)若一个卫星或飞船的运行轨道与地球赤道平面有固定的夹角，并且卫星（飞船）在离地面高度为H的球面S上运行，综合考虑地球自转时该卫星(飞船)在运行过程中相继两圈的经度的差异等因素，计算此时要实现全程跟踪测控的话应该建立的测控站的最少数目；(3)收集我国一个卫星或飞船的运行资料和发射时测控站点的分布信息，并分析这些测控站点对该卫星所能测控的范围，本文选择神舟七号宇宙为研究对象。2 模型假设1.假设地球是规则的球形，半径为6375千米；2.卫星沿闭合的圆形轨道或者椭圆形轨道运行；3.每次监控卫星监控站都能以最大的监控范围内正常工作；4.监控船的位置是灵活的，因此可以动态的监控卫星的运行；5.忽略各个地面监控站的海拔差异，都认为是分布在距离地心为6375千米的地表；6.卫星运行轨迹在地表的投影区域展开近似认为是矩形。3 符号说明 地球的半径 卫星（飞船）运行轨道离地面高度 最少观测站数目（未取整） 最少观测站数目（取整） 即图2中AOB，圆心和观测点的连线与圆心和卫星轨道连线的夹角 即图2中ABO 卫星运行轨道与赤道的夹角 万有引力常量，取值 地球质量，取值 卫星绕地球运行周期 地球的自转周期，即24小时 问题二中的最少测控站数目 地球自转在高度为的线速度 卫星的线速度 测控站测控圆锥在地表投影圆的半径 测控覆盖率4 问题分析在卫星或飞船的发射以及运行过程中对卫星的测控是非常重要的，而其核心问题是测控站点的布设问题。这正是本文要解决的问题。问题一是在所有测控站都与卫星运行轨道共面的假设下进行考虑的，在这种理想的情况下，要想使用最少的测控点，那它们应该是这样分布的：测控点均匀分布；两个相邻测控点辐射出的圆锥的母线恰好相交于卫星轨道平面。在截面中，圆锥母线、地球半径、卫星运行半径三者组成三角形，并且有一个角为93度，于是利用正弦定理、余弦定理即可以求解三角形，测控站数目即为与单个夹角的比值，当然应该取整。问题二是问题一的深入，考虑到了两点因素：一是卫星轨道与赤道呈固定夹角；二是地球自转对飞船相继两圈的经度位置造成的差异。这就决定了卫星的运行轨迹在地表的投影是卫星的运动和地球的自转运动的合运动。实际上将地表展开成近似矩形，卫星的运动轨迹是若干个曲线相交成的网状，其容纳在一个矩形带中。要实现对卫星的全覆盖转化为对这个带的覆盖。然后要探讨观测点的位置问题，讨论发现将测控点设置在赤道上是最节省测控点的。最后改进问题一的算法和模型，求解出最少的探测点数目。另外的一种思路是，将两个运动合成，类比于一个新的卫星绕一个不动的地球在转，求出新卫星的参数，进而求解出要设置的观测点数目。问题三是理论问题一、二的实际应用。本文选取神舟七号宇宙飞船为研究对象。首先要收集起神七的发射、运行、测控等各个方面的详细资料，然后将问题分解为两个部分，一是发射过程的测控问题，二是运行过程的测控问题。对于测控覆盖率的求解，我们将神七留在地表的投影带计算出来，然后将11个测控点加入到投影带图形中，再计算出每个测控站的测控圆的面积，它与投影带面积之比即为覆盖率。对于发射过程，问题的难点在于对详细发射过程、进入预定轨道欠的轨迹等资料的了解，但是这样的资料属于宇航局核心资料是难以查询的，因此我们首先根据比较成熟的计算发射过程轨道的知识和微分方程模拟发射过程，进而以此为依据计算发射过程的测控覆盖率。5 模型的建立与求解5.1 问题一的模型建立在问题一的模型建立与求解的过程中，我们把卫星（飞船）的运行轨道分为圆轨道和椭圆轨道两种情况讨论。5.1.1 情况一：卫星（或飞船）运行轨道为圆轨道假设地球半径为R，卫星（或飞船）飞行轨道离地面的距离为H。实际上，每个测控站的监控范围都是一个以测控站为顶点的曲顶圆锥体，又因为题中假设所有测控站都与卫星的运行轨道共面，所以我们取卫星运行轨道的截面(如图2)作为研究的对象，那么每个测控站的可监控范围的辐射图形为扇形。如果各个观测点均匀分布、且相邻两个观测点的辐射网络在卫星运行轨道切面处恰好交汇，此时所有的测控站辐射范围恰好覆盖卫星的轨道平面，那么此时的观测点数目是最少的，设为N，如图2所示。图2 地球表面与卫星运行轨道截面图从上图可知，在三角形AOB中，由正弦定理可得： （5.1.1）解得： （5.1.2）从而有 （5.1.3）则可得如下关于的表达式： （5.1.4）而N值要取不小于以上计算数值的整数。从（5.1.4）可知，是一个关于H的单调减函数，而又因为N只能取不小于的整数，故。其实际意义在于，卫星离地面越高，地面上的单个测控点辐射的区域会越大，所以需要的测控点就越少，但无论卫星离地面多么远，都至少要用3个测控点才能完全覆盖，具体的计算数值是5.1.4式。图3绘制了观测站点数目与卫星离地高度的函数关系图像，并确定了在N取得整数点时的临界高度值，在图中均已标注。又依据公式5.1.4计算出了在测控点分别为3，4，5，29时，飞船离地面的临界高度H（见表1）。表1 观测站点数目与卫星离地临界高度的关系表 （高度单位：km）观测站的数目34567891011离地最小高度7557.93102.31795.51201.6873.5670.4534.7439.0368.6观测站的数目121314151617181920离地最小高度315.1273.5240.2213.3191.1172.5156.9143.5131.9观测站的数目212223242526272829离地最小高度121.9113.1105.498.592.486.982.077.573.5由上表，我们可用matlab软件画出观测站点数目与卫星离地临界高度的关系图：图3观测站点数目与卫星离地高度的函数图像5.1.2 情况二：卫星（或飞船）运行轨道为椭圆轨道卫星（或飞船）运行轨道为椭圆轨道时，其测控站数目计算结果与卫星（或飞船）运行轨道是圆轨道情形的计算结果是一致的，因为只要在近地点有站点可以探测到，那么到了远地点则一定能探测到，这是因为离地越远，测控站点的探测范围越大，因此测控站点更容易捕捉到离地远的卫星，所以在计算最小站点数目的时候，只要将H赋以近地点距离之值即可。5.2 问题二的模型建立对问题二，我们首先根据条件得到卫星（或飞船）的运行轨迹投影带，把卫星（或飞船）的运行轨道与地球赤道平面的固定夹角看作是一个固定的变量，然后分一般卫星和地球同步卫星两种情况讨论，分别得到了地面观测站数目与夹角和卫星飞行高度的模型，讨论了夹角的变化对地面观测站的影响以及给出了在什么地点建立观测站能达到观测站数目最少，并给出了几类常见卫星的高度与地面测控站数目图。5.2.1 模型一 图4 卫星运行与地球表面模拟图如图4所示，假设卫星的运行轨道与赤道平面的夹角为，卫星在离地面高度为H的球面上运行。由万有引力定律，可以得到卫星运行的线速度V和运行周期T。其中M是地球质量，；G是万有引力常量，。由于夹角是固定不变的，因此在地球不发生自转的情况下，卫星的运行轨迹在地球表面的投影是一个固定的圆。但是由于在卫星绕地球运行的同时地球自身在发生自转，因此投影到地球表面的轨迹是地球的自转和卫星的旋转这两个运动的合成。图5为卫星在地球表面上的实际投影。所以实际上，在卫星飞行一周后，地球已经自转了时间T（卫星的周期），因此两条轨迹的投影会相隔一段距离S： （5.2.1）其中R是地球平均半径。我们更深入的研究可知，如果地球是球形的，且质量分布均匀，卫星绕地球按照圆轨道飞行，则在地球自转一周的时间里，卫星可以飞行24/T圈，投影到地球表面就会形成24/T条圆形轨迹。相邻两条轨迹与赤道的交点之间的距离都为S。如图5中-、-、-、-、-间的距离。相邻两圈的经度差异为 （5.2.2）图5 卫星的运行轨迹投影带关于图5的解释说明：这是地球表面的一个展开图，横坐标轴是赤道，端点分别是东经180度和西经180度；曲线是卫星运行轨迹在地表的投影；每一条完整的曲线代表卫星绕地球运行一周的投影轨迹；曲线的数目为24/T。由于卫星运行轨道平面与赤道平面夹角为，则卫星的运动轨迹在地面的投影必定分布在北纬和南纬之间。由图6可知，弧长与半径的关系为：根据两个扇形相似性，可以得到：图6 地球与卫星相对位置所以卫星在天球面S上的运行轨迹处在宽度为S1的带状环形区域内，并且 （5.2.3）测控站所选取的位置不同，会使得监控区域的重叠和相对位置发生变化，进而影响到所需测控站的数目。下面分别考虑把测控站建在边缘区域和赤道的两种特殊情况。情形一：当观测站恰好位于这个带状环形区域的边缘时候，即处在北纬度或南纬度时，可以计算出当观测站完全覆盖高度为H的卫星在这一点可能出现的地方时，的最大角度（见图7）。图7 情形一：观测站建在宽带边缘由余弦定理可得：解上面的方程组可得：即将上式带入方程组消去，可得求解得：将回代至方程组中，可以求得： (5.2.4)当大于上式的计算结果时，建在边缘的观测站无法监控到卫星运行时可能出现的全部区域；当小于上式的计算结果时，建在边缘的观测站可以完全监控到卫星在此点上空运行时可能到达的任何区域。情形二：当观测站建在赤道上时。如下图，ED是赤道平面，C是观测站位置，ACD和BCD都等于87度，GED和FED 为。CE=R，CD=H，弧GDF位于球面S上。图8 情形二：观测站建在赤道上这种情形的情况是直观的，因为是卫星轨道和赤道平面的交角，所以；而且观测站的观测角度有87度，因此在观测站上空的圆锥体范围内，无论卫星轨迹如何移动，都会被监控站监控到。将两种情形作对比，若将测控站建在投影带的边缘，不妨假设测控站建在北纬度处，那么此时要达到完全覆盖的目的，就一定要覆盖到到南纬度处（投影带的下边缘），那么其跨度达到了2；而若将测控站建在赤道处，那么要实现完全覆盖的话，就要覆盖到南北纬度处，而此时卫星运动轨迹在地球上投影具有的对称性，观测站可以最大限度的监控卫星在圆环带区域上空可能出现的位置，那么其跨度只有。所以，很明显的，如果将测控站建在投影带的中线即赤道处是最理想的，即所需要测控站数目最少的情形。现在我们继续深入探讨，考虑当卫星在离地面高度为H的球面S运行时，需要多少个站点才能完全监控住卫星可能出现的区域。由于卫星的运动轨道与地球赤道所成的角度是固定的。前面已经根据运动的合成分析过了，卫星的轨迹在地球表面的投影是一个圆环带。卫星运行轨道面与赤道面的斜交角度为，我们将测控站都建立在赤道上。假设最少需要个观测站才能恰好完全覆盖卫星所到达的区域，这个时候应该恰好是当观测站1结束对卫星的监控时，观测站2恰好收到该卫星的信号。利用这个原理我们可以计算得到和地球半径R、轨道高度H的关系式。 图9上图中，由余弦定理可得又由正弦定理可得 联立以上三个方程，解方程组得 (5.2.5)从而我们得到了一个关于半径、高度与地面测控站数目的数学模型。现在我们取地球半径为，对一些有代表性的卫星，利用它们的数据画出上面的函数关系如下：图10 观测点最少数目与卫星离地面高度关系5.2.2 模型二因为卫星在绕着地球转，而地球同时也在自转，所以最终形态是两个运动的合成。我们考虑将两个运动合成，因为地球的自转对卫星的投影轨迹产生了影响（出现了纬度差异），我们将地球的自转合成到卫星的运行中，即将地球假定为不再自转，而此时卫星的运动与真实的运动便不再相同了，等价于另外一颗不同参数（运行高度、周期等）的卫星在绕着不能自转的地球运行。假设新卫星的运行高度为，线速度为。情形一：卫星呈一定锐角倾斜角自西向东运行图11 同向旋转时运动的合成此时，卫星的旋转与地球的自转是同向的，那么新卫星的线速度为： （5.2.6）由万有引力定律可知，而 因此，又因为 故得 （5.2.7）又根据问题一的结论5.1.4式，有： 将式5.2.7代入到式5.1.4中，可以得到： （5.2.8）当然，实际的N值要取不小于的整数值。情形二：卫星呈一定锐角倾斜角自东向西运行图12 反向旋转时运动的合成此时，卫星的旋转与地球的自转是反同向的，那么新卫星的线速度为： （5.2.9）于是得 从而得 (5.2.10)同理，将5.2.10代入5.1.4，得到： （5.2.11）5.3 问题三的模型建立及求解5.3.1神七相关的资料及数据：神舟七号载人飞船是中国神舟号飞船系列之一，于2008年9月25日21点10分04秒988毫秒从中国酒泉卫星发射中心载人航天发射场用长征二号F火箭发射升空。飞船于2008年9月28日17点37分成功着陆于中国内蒙古四子王旗主着陆场。神舟七号飞船共计飞行2天20小时27分钟。神舟七号飞船实现了中国人在太空中的第一次行走，具有划时代的意义。神舟七号的运行轨道与赤道的夹角为42，其近地点距地面距离为200km，远地点距地面距离为343km。神舟七号发射和运行过程中设立了11个固定站点、开动了5艘远望船舰对其进行跟踪测控，这11个站点的经纬度等详细信息如表3.1，为研究方便，本文已经对其进行了编号。表2 各个监控站的地理位置站 点位置1北京站2喀什站3和田站4东风站5青岛站6渭南站经度11623E7559E79E10110E12022E10930E纬度3954N3928N3707N42N3603N3414N站点位置7厦门站8纳米比亚站9卡拉奇站10马林迪站11圣地亚哥站远洋一、二、三、五、六号经度11804E1829E6702E405E7027W没有固定位置，分布在各大洋，相继移动追踪卫星信号纬度2426N2257S2451N317S3326S5.3.2飞船在发射过程中的测控从神舟七号的发射日志中我们了解到如下关于神舟七号飞船的发射资料：21时09分许：火箭点火21时10分：神舟七号飞船升空点火第120秒,，火箭抛掉逃逸塔点火第159秒 ，火箭一二级分离成功点火第200秒，整流罩分离点火第500秒，二级火箭关机点火第583秒时，飞船与火箭成功分离飞船在上升阶段有三个站，第一个是东风站，就是发射场；第二个是渭南站，在西安附近；第三个是青岛站。这三个测控站负责飞船在上升段的测量，因为在上升阶段火箭一直在中国境内（更确切的说，一直在发射场附近），所以三个测控站实现了100%的测控覆盖率。其次是入轨阶段，有两条测量船“远望一号”和“远望二号”进行实时跟踪测控。飞船入轨的时候有很多动作，如捕获地球、建立正常运行姿态、太阳帆板要展开、判断轨道是否正确等，因为远望号是可以调整位置的，所以覆盖率也达到了100%。青岛站在入轨后1分钟还可以看，和“远望一号”测量船可以接上。这样，飞船入轨以后5到6分钟的情况地面都可以完全监测到。入轨二十分钟以后，“远望二号” 船再进一步跟踪判断飞船入轨运行情况。飞船入轨以后的测控情况参见下面模型的详细分析。5.3.3 飞船在运行过程中的测控因为神州七号宇宙飞船的运行轨迹与赤道的夹角为42，所以由问题二的结论可知，飞船的轨迹在地表的投影是北纬42和南纬42之间的宽带，将这个宽带展开，近似地看作一个矩形，则矩形的长为赤道长度，宽度为南北纬42之间的弧长，即然后将这11个固定站点定位到这个长为、宽为的宽带当中去。其具体方法如下：以展开地球的南纬42纬度线为横坐标轴、以经线为纵坐标轴建立坐标系，横纵坐标均以距离为度量，单位为千米。将地球近似看做球形，则每一度经度、每一度纬度的跨越距离均为。基于此，将上述11个坐标转化为如下表3的坐标。表3 转化后的站点的坐标序号123456纵坐标9005.48957.78699.39236.385828382.3横坐标9814645069278267107019954序号7891011纵坐标7350.54257.19427304.72094.7横坐标6688440022899118191871至此便完成了11个测控点到地面的投影。然后计算每个测控点的测控半径，实际上计算的是每个测控点所辐射的圆锥在地表投影的圆的半径，由图13可得方程组：图13根据上面得方程组可得，解以上方程，可得 ， （舍去）在中，由余弦定理可得进一步，得（弧度），换算成角度是11.4206度。最后在扇形AOB中，求弦长AB，即我们要求的观测站观测范围在地面投影的半径：也就是说每个站点所辐射的圆锥投影到地面的圆的半径为1268.6km，因此可以在宽带中以每个站点为圆心，以1268.6km为半径作出每一个圆，如下图所示，红色线分别为北纬42和南纬42对应的线。 图14 测控站点测控区域地面投影图由上图可知，要计算测控站点的覆盖率，可以转化为近似的计算测控站点所辐射的圆的面积与宽带总面积的比例。宽带的面积为 （5.3.1）圆所覆盖面积计算如下：从图14可以看出，有两个完整的圆包含在宽带中并与其他圆不相交，面积为 （5.3.2）有一个圆有一部分落在了宽带内部，其方程为其与x轴交点为（23749，0）、（22049，0），并且它落入宽带中的部分圆的面积为： （5.3.3）余下的8个圆相互交叠，需要计算其覆盖的面积。首先得到每一个圆的方程： （5.3.4）以上8个圆的交叠中共有6个关键点的交点，它们决定了覆盖区域的轮廓。又因为有两个交点非常接近，为了便于计算，我们将其视为同一个点。另外还有两个重要的点是其与北纬42线的交点，将其覆盖面放大如下图：图15 交叠圆放大、分割图计算方程组（5.3.4），可以得到六个边界点坐标分别为A（5604.5，8011.06）、B（7859.71，7838.92）、C（8728.24，8054.0175）、D（10551，7262.2）、E（11970，8565.0）；与北纬42线的两个交点坐标为G（5211.55，9236）、F（11788.4，9236）。将以上重叠区域分割为8个部分，分别为GKA、KABH、AB、BHIC、ICDJ、JDEF、DE、EF。因此重叠部分面积计算为： （5.3.5）其中代表线段AB下方弓形的面积，代表线段DE下方弓形的面积，而弧GA、弧BC、弧CD、弧EF的曲线方程见5.3.4式。代入实际数据得：所以测控覆盖率为： （5.3.6）另外，以上计算的只是11个固定的测控点的测控覆盖率，而神舟七号在发射和运行过程中还有5艘远望船舰。因为这5艘船是可以随时变动地点的，所以我们近似认为其测控辐射的圆都是完全落在宽带中的，则有，那么此时全部16个测控点的覆盖率为： （5.3.7）由以上的建模和分析，我们发现，在发射过程中，由于上升阶段火箭偏离发射中心的距离并不大而且附近布设的测控点数目又多，所以基本可以达到100%的测控率；而在运行阶段，测控的覆盖率维持在13%至25%的范围，但是由于在本次神七发射过程中“天链一号”中继卫星的同时发射，使得覆盖率远远上升，达到了60%以上，可以较好的完成测控任务。6 模型的评价与改进6.1模型的优点1.本文依据严谨的算法建立了科学恰当的数学模型，成功解决了所提出的三个问题；2.在解决问题二时，本文从两个不同的角度出发，一个是从一个投影的角度考虑，一个是从运动合成的角度考虑，建立了两个完全不同的模型，但均取得了良好的计算结果；3.问题三的求解中，本文对发射过程和运行过程分别讨论，分别建立模型进行了求解；4.在问题三的求解过程中，本文成功的画出了每个测控点辐射圆锥在地球表面的投影，并且在计算重叠区域的面积时，本文使用了图像分割、定积分求解的方法精确的计算出了重叠区域面积。6.2模型的缺点及改进1.在问题一中我们假设了观测站是均匀分布的，这与实际情况是不一致的。因为观测站的选址要考虑地形、气候等多方面影响，理论上适合建立观测站的位置可能由于地形、气候、人文条件等因素的影响而无法实现；2.问题二中，如果测控站投影直径远小于宽带的宽度，那么问题将会转化成如何用最少的圆去覆盖固定矩形的类似于最小覆盖问题，本文并没有对此做更深入的探讨，这是更进一步改进的一个重点方面。3.在问题三中，虽然本文将发射过程和运行过程的测控进行了分开考虑，但是由于火箭发射过程中行走的轨道是完全不可知的（这是航天专家组长期研究制定的），而且这段时间火箭受到各种作用而不断的改变飞行姿态和角度，因而仿真也是无法实现的，所以极难建立准确模型求解。庆幸的是，我们查阅资料得知了发射过程时测控率是100%的。如果能有升空阶段轨道的微分方程，那么问题便会变得迎刃而解。附 录1. 卫星监控站的理想分布图t=.0:.01:2*pi;x=4*cos(t);y=4*sin(t);x1=10*cos(t);y1=10*sin(t);x2=3*cos(t);y2=18*sin(t)-4;plot(x,y,x1,y1,x2,y2)axis offtitle(卫星监控站的理想分布状态)2. 卫星近圆轨道的动力行为模拟x=2.5; y=0.2;dt=0.001;A=1.5;B=3.6;for i=1:15000 x1=x+(A-(B+1)*x+x2*y)*dt; y1=y+(B*x-x2*y)*dt; plot(x1,y1); hold on; x=x1;y=y1;end3. 卫星的动力系统函数文件function dx=myfun(t,x)dx=zeros(2,1);dx(1)=x(2);dx(2)=6*sin(2*t)/(1+3*sin(t)*sin(x(2)2-0.25*sqrt(1+sin(2*t)*sin(t);4. 求解卫星动力系统的代码clear;x45 y45=ode45(myfun,0:20:500,0.0 0.0);plot(x45,y45,.-)ylabel(飞船离地面高度/千米)xlabel(飞船发射后各个时刻/秒)5. 同一平面内均匀分布的监测站和卫星高度关系H=6375:.1:10000;N=pi./(87*pi/180-asin(6371*sin(93*pi/180)./H);plot(H,N*10(-15),.-)6. 第二问中至少所需的监控站数目R=6300;n1=2:30;H1=R./(cot(93*pi/180)*sin(pi./n1)+cos(pi./n1)-R;n2=3:30;H2=R./(cot(93*pi/180)*sin(pi./n2)+cos(pi./n2)-R;subplot(1,2,1),plot(n1,H1,*-);xlabel(至少所需的观测站数目n)ylabel(卫星离地面的高度H)grid onsubplot(1,2,2),plot(n2,H2,-);xlabel(至少所需的观测站数目n)ylabel(卫星离地面的高度H)grid on7. 第三问中监控站监控范围和卫星轨迹投影区域示意图%由纬度计算纵坐标guancezhan=81.9 81.4667 79.1167 84 78.05 76.2333 66.85 38.7167 8.5667 66.4333 19.0500;s = 6300*2*pi*guancezhan/360%由经度计算横坐标，以0度经线为纵轴jingdu=116+23/60 75+59/60 79 101+10/60 120+22/60 109+30/60 67+2/60 40+5/60 250-27/60 118+4/60 18+29/60;weidu=39+54/60 39+28/60 37+7/60 42 36+3/60 34+14/60 24+51/60 3+17/60 33+26/60 24+26/60 22+57/60;hengzuobiao=2*pi*6300*cos(weidu*pi/180).*jingdu/360plot(hengzuobiao,s,*);hold onfor i=1:11 t=0:.01:2*pi; x11=hengzuobiao(i)+1269.2*cos(t); y11=s(i)+1269.2*sin(t); plot(x11,y11); hold onendhold onm=0:100:25000;n=9236.3*ones(1,size(m);plot(m,n,.-)xlabel(南纬42度)ylabel(0度经线)grid on23