【Cesium】使用TLE轨道两行数计算轨道信息，并生成CZML格式文件

TLE为轨道两行数，简单的说是用两行数字表示轨道的相关信息，本文即用轨道两行数来计算任一时刻卫星的位置信息和速度信息，并生成CZML文件能够读取的格式

1、satellite.js库简介

简而言之，satellite.js库可以根据TLE轨道两行数，使用SGP4/SDP4方法对任一时刻卫星的位置信息和速度信息进行计算，也可以完成地惯坐标系和地固坐标系下位置信息的转换等其他功能，本篇文章主要使用到satellite.js库的轨道计算功能，感兴趣的读者可以看一下这个库中其他的功能，链接放在下面
satellite.js库介绍
中文翻译
举例：

// 示例TLE轨道两行数var tleLine1 ='1 25544U 98067A   19156.50900463  .00003075  00000-0  59442-4 0  9992',
    tleLine2 ='2 25544  51.6433  59.2583 0008217  16.4489 347.6017 15.51174618173442';// 初始化卫星记录var satrec = satellite.twoline2satrec(tleLine1, tleLine2);//  获取卫星在某一时间点的位置与方向信息 **（注意时间是从TLE记录的初始时刻后的分钟数，可以为负）**var positionAndVelocity = satellite.sgp4(satrec, timeSinceTleEpochMinutes);//  或者使用JavaScript日期var positionAndVelocity = satellite.propagate(satrec,newDate());// 获得地惯坐标系下的当前时间点的位置信息和速度信息var positionEci = positionAndVelocity.position,
    velocityEci = positionAndVelocity.velocity;// 将观察者（固定区域）设为西经122.03度，北纬36.96度（degreesToRadians——角度转弧度）var observerGd ={longitude: satellite.degreesToRadians(-122.0308),latitude: satellite.degreesToRadians(36.9613422),height:0.370};// 将GMT中国时、Julindate时转成格林尼治恒星时var gmst = satellite.gstime(newDate());// 地惯坐标系下的位置信息转换成地固坐标系下的位置信息var positionEcf   = satellite.eciToEcf(positionEci, gmst),// 大地坐标系转地固坐标系
    observerEcf   = satellite.geodeticToEcf(observerGd),//  地惯坐标系转大地坐标系
    positionGd    = satellite.eciToGeodetic(positionEci, gmst),//  获得观测角度
    lookAngles    = satellite.ecfToLookAngles(observerGd, positionEcf),// 获得笛卡尔坐标中的x，y，zvar satelliteX = positionEci.x,
    satelliteY = positionEci.y,
    satelliteZ = positionEci.z;// 角度可以通过' azimuth '， ' elevation '， ' range_sat '属性访问。var azimuth   = lookAngles.azimuth,
    elevation = lookAngles.elevation,
    rangeSat  = lookAngles.rangeSat;// 通过'经度'，'纬度'，'高度'访问大地坐标。var longitude = positionGd.longitude,
    latitude  = positionGd.latitude,
    height    = positionGd.height;//  将弧度转换为度。var longitudeDeg = satellite.degreesLong(longitude),
    latitudeDeg  = satellite.degreesLat(latitude);

2、生成CZML能够识别的格式

使用CZML文件表示卫星，不仅需要知道卫星的位置信息，还需要知道卫星在运行过程中每一圈轨道的运行时间，这样画出来的轨道才能够是一圈一圈独立的

不然会出现轨道交叉的现象

其实根本在于要控制好卫星轨迹提前出现的时间，以及轨迹保留的时间，也即leadIntervalArray和trailIntervalArray这两个参数，具体的代码如下：

const satellite =require('satellite.js')const moment =require('moment')const julian =require('julian')/*
根据卫星显示的起始时间，终止时间，tle轨道两行数得出czml文件，时间为js的Date对象，tles为对象数组，对象格式为
{
name：xx，
tle1：xx，
tle2：xx
}
*/functiontles2czml(startTime, endTime, tles){// 计算起始时间和终止时间相隔的分钟数let minsInDuration =(endTime.getTime()- startTime.getTime())/60000;//mins//设置为开始时间let initialTime =moment(startTime.toISOString()).toISOString();//设置为结束时间
    endTime =moment(endTime.toISOString()).toISOString();// 初始化czml数据，创建场景信息let tempCZML =[];
    tempCZML.push({"id":"document","name":"CZML Point - Time Dynamic","version":"1.0","clock":{"interval":`${initialTime}/${endTime}`,"multiplier":1,"range":"LOOP_STOP","step":"SYSTEM_CLOCK"}},)// 处理每一个satfor(let no =0; no < tles.length; no++){if(!tles[no].name){
            console.log("请输入第"+ no+1+"个卫星的名称");return};if(!tles[no].tle1){
            console.log("请输入第"+ no+1+"个卫星的第一个两行数");return};if(!tles[no].tle2){
            console.log("请输入第"+ no+1+"个卫星的第二个两行数");return};let sat_name = tles[no].name;// 保存位置信息let res =[];let satrec 
        satrec = satellite.twoline2satrec(tles[no].tle1, tles[no].tle2);//satrec.no：以弧度/分钟为单位的平均运动，一天有1440分钟，一弧度是0.159155圈// to go from RAD/DAY -> REV/DAY: rad * 1440 * 0.159155//to go from REV/PER DAY to MINS/REV -> 1440/RevPerDaylet totalIntervalsInDay = satrec.no *1440*0.159155;//1440 = min && 0.159155 = 1turn// 获得运行一圈的分钟数let minsPerInterval =1440/ totalIntervalsInDay;// mins for 1 revolution around earth// intervalTime 取结束时间 格式为2008-09-20T12:25:40.104Zlet intervalTime = endTime

        let leadIntervalArray =[];let trailIntervalArray =[];
        console.log("Setting intervals...");// 注意：这里之所以要倒过来求leadInterval和trailInterval是因为如果正着求，很有可能在终止时刻卫星并没有运行完一圈，导致轨道只显示一半for(let i = minsInDuration; i >=0; i -= minsPerInterval){if(i <= minsPerInterval){// intial interval let currentOrbitalInterval ={"interval":`${initialTime}/${intervalTime}`,"epoch":`${initialTime}`,"number":[0, minsPerInterval *60,
                        minsPerInterval *60,0]}let currTrail ={"interval":`${initialTime}/${intervalTime}`,"epoch":`${initialTime}`,"number":[0,0,
                        minsPerInterval *60, minsPerInterval *60]}
                leadIntervalArray.push(currentOrbitalInterval);
                trailIntervalArray.push(currTrail);}else{//not initial so make intervals let previousIntervalTime =moment(intervalTime).add(-minsPerInterval,'m').toISOString();let currentOrbitalInterval ={"interval":`${previousIntervalTime}/${intervalTime}`,"epoch":`${previousIntervalTime}`,"number":[0, minsPerInterval *60,
                        minsPerInterval *60,0]}let currTrail ={"interval":`${previousIntervalTime}/${intervalTime}`,"epoch":`${previousIntervalTime}`,"number":[0,0,
                        minsPerInterval *60, minsPerInterval *60]}
                intervalTime =moment(intervalTime).add(-minsPerInterval,'m').toISOString();
                leadIntervalArray.push(currentOrbitalInterval);
                trailIntervalArray.push(currTrail);}}// Seconds between current time and epoch timelet sec =(startTime - julian.toDate(satrec.jdsatepoch))/1000;
        console.log(startTime, julian.toDate(satrec.jdsatepoch), sec);for(let i = sec; i <= sec + minsInDuration *60; i++){//每60秒计算一个位置信息，最后采用拉格朗日插值法处理数据// 根据当前时间距tle两行数历元时刻的分钟数，计算当前卫星位置和速度let positionAndVelocity = satellite.sgp4(satrec, i *0.0166667);// 0.0166667min = 1sec// 地惯坐标系let positionEci = positionAndVelocity.position;
            positionEci.x = positionEci.x *1000;
            positionEci.y = positionEci.y *1000;
            positionEci.z = positionEci.z *1000;// let velocityEci = positionAndVelocity.velocity;// velocityEci.x = velocityEci.x * 1000;// velocityEci.y = velocityEci.y * 1000;// velocityEci.z = velocityEci.z * 1000;

            res.push(i - sec, positionEci.x, positionEci.y, positionEci.z);}let initialCZMLProps ={"id":`${sat_name}`,"name":`${sat_name}`,"availability":`${initialTime}/${endTime}`,"label":{"fillColor":{"rgba":[255,0,255,255]},"font":"11pt Lucida Console","horizontalOrigin":"LEFT","outlineColor":{"rgba":[0,0,0,255]},"outlineWidth":2,"pixelOffset":{"cartesian2":[12,0]},"show":true,"style":"FILL_AND_OUTLINE","text":`${sat_name}`,"verticalOrigin":"CENTER"},"path":{"show":[{"interval":`${initialTime}/${endTime}`,"boolean":true}],"width":3,"material":{"solidColor":{"color":{"rgba":[// 随机生成轨道颜色
                                Math.floor(255* Math.random(0,1)), Math.floor(255* Math.random(0,1)), Math.floor(255* Math.random(0,1)),255]}}},"resolution":120,// The time ahead of the animation time, in seconds, to show the path."leadTime": leadIntervalArray,// The time behind the animation time, in seconds, to show the"trailTime": trailIntervalArray
            },"model":{"show":true,"gltf":"./111.gltf","minimumPixelSize":50,},"position":{// 采用拉格朗日插值法"interpolationAlgorithm":"LAGRANGE",// 1为线性插值，2为平方插值"interpolationDegree":2,// 参考坐标系，地惯坐标系"referenceFrame":"INERTIAL","epoch":`${initialTime}`,"cartesian": res
            }}
        tempCZML.push(initialCZMLProps);}return tempCZML;}exportdefault tles2czml

上面的方法参考了tle2czml库的写法，但是有一点需要注意，如果直接调用tle2czml库的话，最后一圈的卫星轨道往往是不完整的，本文所写的方法是根据时间倒序求leadInterval和trailInterval，这样即便到了终止时间，卫星的轨道显示也是完整的一圈。
如果对CZML文件格式有疑问的可以看 这里
欢迎交流~