本教程将解释如何初始化Rugged并使用它来对卫星图像进行地理定位。假设传感器是一个具有2000个像素和20°视场的单行成像仪,偏离天顶角10°。GPS和AOCS辅助数据可用,并为我们提供了在获取期间记录的位置、速度和姿态四元数的列表。通过将所有这些信息传递给Rugged,我们将能够精确地定位图像上的每个点在地球上的位置。嗯,并不是完全精确,因为这个第一个教程没有使用数字高程模型,而是将地球视为一个椭球体。DEM将在第二个教程中添加:使用DEM进行直接定位。这里的目标仅限于解释如何初始化Rugged对传感器和获取的所有信息。
让我们从定义成像仪开始。传感器模型由其观测几何(即每个物理像素的视线)和其日期模型描述。
我们需要定义卫星坐标系中每个传感器像素的视线(LOS)向量坐标;X轴与没有指向的情况下的速度向量平行,Z指向地球,Y使得X、Y、Z形成一个右手坐标系。
为此,我们需要以下包
import org.hipparchus.geometry.euclidean.threed.Rotation;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import java.util.ArrayList;
import java.util.List;
import org.orekit.rugged.los.LOSBuilder;
import org.orekit.rugged.los.FixedRotation;
import org.orekit.rugged.los.TimeDependentLOS;
像素i相对于仪器的原始视线方向由向量定义:
List<Vector3D> rawDirs = new ArrayList<Vector3D>();
for (int i = 0; i < 2000; i++) {
// 20°视场,2000个像素
rawDirs.add(new Vector3D(0d, i*FastMath.toRadians(20)/2000d, 1d));
}
仪器相对于轨道的X轴偏离10°,我们需要旋转视线方向以获得卫星坐标系中的视线
LOSBuilder losBuilder = new LOSBuilder(rawDirs);
losBuilder.addTransform(new FixedRotation("10-degrees-rotation", Vector3D.PLUS_I, FastMath.toRadians(10)));
在这里,我们考虑到视线方向在时间上是恒定的,也可以通过使用其他变换来获得时间相关的视线。还可以在原始方向和最终视线之间附加多个变换。
TimeDependentLOS lineOfSight = losBuilder.build();
日期模型给出图像中的行号作为时间的函数,反之亦然。这里我们使用线性日期模型(行采样率是恒定的)
我们使用Orekit处理时间和日期,并使用Rugged定义日期模型:
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScalesFactory;
import org.orekit.rugged.linesensor.LinearLineDatation;
AbsoluteDate absDate = new AbsoluteDate("2009-12-11T10:49:55.899994", TimeScalesFactory.getGPS());
LinearLineDatation lineDatation = new LinearLineDatation(absDate, 1d, 20);
LinearLineDataion类使用3个参数实例化:第一个是参考日期,第二个是参考日期的行号,最后一个参数是行速率(此示例中为每秒20行)。
请注意,Orekit可以通过TimeScalesFactory处理许多不同的时间尺度,因此我们不需要担心转换。
现在,我们已经定义了LOS和数据,可以在Rugged中初始化一个线传感器对象:
import org.orekit.rugged.linesensor.LineSensor;
LineSensor lineSensor = new LineSensor("mySensor", lineDatation, Vector3D.ZERO, lineOfSight);
第一个参数是传感器的昵称。这是必要的,因为我们可以定义多个传感器(对于具有多个光谱通道或探测器阵列的成像仪很有用)。第三个参数是传感器相对于卫星参考系中心的位置。在这里设置为0。
作为Rugged的输入,我们需要传递足够的信息来描述在获取期间平台的位置和方向。在我们的示例中,位置、速度和四元数的列表是硬编码的。在现实生活中,我们将从卫星辅助遥测中提取GPS和AOCS数据。
请注意,为了模拟目的,我们也可以使用Orekit来模拟轨道。它提供了非常方便的功能,用于传播具有偏航补偿的太阳同步轨道(地球观测卫星的典型轨道)。
Rugged期望位置(单位:米)、速度(单位:米/秒)和四元数以惯性坐标系表示。Orekit实现了所有标准的惯性和地球坐标系,因此我们只需要指定我们使用的坐标系并进行必要的转换。
从惯性坐标系到地球旋转坐标系的转换对用户来说是透明的,它基于最新的岁差/章动模型,并应用了IERS发布的修正。IERS公告和其他物理数据提供在orekit数据文件夹中。有几种配置Orekit以使用这些数据的方法。更多信息请参见这里
在我们的应用程序中,我们只需要知道我们正在使用的坐标系的名称。位置和速度以ITRF地球坐标系给出,而四元数以EME2000惯性坐标系给出。
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.utils.IERSConventions;
Frame eme2000 = FramesFactory.getEME2000();
boolean simpleEOP = true; // 我们不想计算毫米级别的微小潮汐效应
Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, simpleEOP);
姿态四元数被分组在一个TimeStampedAngularCoordinates对象的列表中:
import org.orekit.utils.TimeStampedAngularCoordinates;
List<TimeStampedAngularCoordinates> satelliteQList = new ArrayList<TimeStampedAngularCoordinates>();
每个姿态样本(四元数、时间)都被添加到列表中,
AbsoluteDate attitudeDate = new AbsoluteDate(gpsDateAsString, TimeScalesFactory.getGPS());
Rotation rotation = new Rotation(q0, q1, q2, q3, true); // q0是标量项
Vector3D rotationRate = Vector3D.ZERO;
Vector3D rotationAcceleration = Vector3D.ZERO;
TimeStampedAngularCoordinates pair = new TimeStampedAngularCoordinates(attitudeDate, rotation, rotationRate, rotationAcceleration);
satelliteQList.add(pair);
其中,例如gpsDateAsString设置为“2009-12-11T10:49:55.899994”
同样,位置和速度将被设置在一个TimeStampedPVCoordinates
的列表中。在添加到列表之前,它们必须被转换为EME2000坐标系:
import org.orekit.utils.TimeStampedPVCoordinates;
import org.orekit.utils.PVCoordinates;
import org.orekit.frames.Transform;
List<TimeStampedPVCoordinates> satellitePVList = new ArrayList<TimeStampedPVCoordinates>();
AbsoluteDate ephemerisDate = new AbsoluteDate(gpsDateAsString, TimeScalesFactory.getGPS());
Vector3D position = new Vector3D(px, py, pz); // 在ITRF坐标系中,单位:米
Vector3D velocity = new Vector3D(vx, vy, vz); // 在ITRF坐标系中,单位:米/秒
PVCoordinates pvITRF = new PVCoordinates(position, velocity);
Transform transform = itrf.getTransformTo(eme2000, ephemerisDate);
PVCoordinates pvEME2000 = transform.transformPVCoordinates(pvITRF);
satellitePVList.add(new TimeStampedPVCoordinates(ephemerisDate, pvEME2000.getPosition(), pvEME2000.getVelocity(), Vector3D.ZERO));
最后我们可以初始化Rugged。代码如下:
import org.orekit.rugged.api.AlgorithmId;
import org.orekit.rugged.api.BodyRotatingFrameId;
import org.orekit.rugged.api.EllipsoidId;
import org.orekit.rugged.api.InertialFrameId;
import org.orekit.rugged.api.Rugged;
import org.orekit.rugged.api.RuggedBuilder;
import org.orekit.rugged.errors.RuggedException;
import org.orekit.utils.AngularDerivativesFilter;
import org.orekit.utils.CartesianDerivativesFilter;
Rugged rugged = new RuggedBuilder().
setAlgorithm(demAlgoId).
setDigitalElevationModel(demTileUpdater, nbTiles).
setEllipsoid(EllipsoidId.WGS84, BodyRotatingFrameId.ITRF).
setTimeSpan(acquisitionStartDate, acquisitionStopDate, tStep, timeTolerance).
setTrajectory(InertialFrameId.EME2000,
satellitePVList, nbPVPoints, CartesianDerivativesFilter.USE_P,
satelliteQList, nbQPoints, AngularDerivativesFilter.USE_R).
addLineSensor(lineSensor).
build();
哎呀,听起来很复杂。其实并不难,因为我们已经定义了大部分所需的内容。让我们逐行描述Rugged实例的构建过程:
由于Rugged作为顶层对象将用于所有用户交互,它相当复杂,并且可以以几种不同的方式构建,因此采用了构建器模式。构建器本身通过调用专用的setter方法来配置,每个概念(数字高程模型、交点算法、轨迹等)都有对应的setter方法。由于所有这些概念都可以链接在一起,setter方法本身实现了流畅接口模式,这意味着每个setter方法都返回构建器实例,因此可以直接调用另一个setter方法。
setAlgorithm方法指定要使用的交点算法。由于本教程旨在非常简单,我们选择直接使用椭球体而不是真实的数字高程模型,因此可以将AlgorithmId.IGNORE_DEM_USE_ELLIPSOID
作为该方法的参数。
setDigitalElevationModel方法指定数字高程模型。实际上,由于我们决定在本教程中忽略数字高程模型,我们可以省略这个调用,程序仍然可以正常工作。我们选择保留这个调用,以防止用户忘记为真正使用数字高程模型的交点算法设置数字高程模型。由于将忽略该模型,我们可以将该方法的参数设置为null
和0
。当然,如果选择了另一个算法,null参数显然不会起作用,这在另一个教程中有解释:使用数字高程模型进行直接定位。
setEllipsoid方法定义了椭球体的形状和方向。我们使用简单的预定义枚举值:EllipsoidId.WGS84
,InertialFrameId.EME2000
,但也可以使用自定义的椭球体。
setTimeSpan方法用于定义搜索算法(直接和反向定位)的全局时间跨度。有四个参数用于此目的:acquisitionStartDate(获取开始日期),acquisitionStopDate(获取停止日期),tStep(预计算帧变换缓存将填充的步长),timeTolerance(反向定位期间允许外推的边际,以秒为单位)。tStep参数是实现良好性能的关键参数,因为帧变换将在整个时间跨度内使用此时间步长进行预计算。这些计算是昂贵的,因为它们涉及地球岁差/章动模型的评估。因此,变换被预计算和缓存,而不是为每个像素重新计算。然而,如果预计计算的变换预计在缩短的时间跨度内进行直接和反向定位(例如几十秒),那么以毫秒为单位的一半轨道上的变换预计算将是一种计算能力的浪费。因此,典型的值是尽可能限制时间跨度,以正确覆盖预期的直接和反向定位调用,并使用一毫秒到一秒之间的步长,具体取决于所需的精度。要使用的确切值取决于任务。最后的timeTolerance参数只是在最终预计算变换之前和之后使用的边际,以允许轻微外推,如果在搜索期间间隔略微超出。典型的值是允许几个图像行,例如5行的容差将意味着计算容差为:timeTolerance = 5 / lineSensor.getRate(0))。
setTrajectory方法定义了航天器的演化。参数是带有时间戳的位置和速度的列表,以及它们相对于惯性坐标系的定义和插值选项:要使用的点数和导数的滤波器类型。插值多项式对于没有任何导数的nbPVPoints(CartesianDerivativesFilter.USE_P
的情况:仅使用位置,不使用速度)具有度数nbPVPoints - 1。在包括速度的计算中(CartesianDerivativesFilter.USE_PV
的情况),插值多项式具有度数2 * nbPVPoints - 1。如果位置/速度数据质量良好且相隔几秒钟,可以选择仅使用少量点进行插值,但同时使用位置和速度进行插值;在其他情况下,可以选择更多点,但仅使用位置进行插值。对于姿态四元数也有类似的参数。
最后一个使用的setter方法,addLineSensor,注册了一个线传感器。从其前缀(add
而不是set
)可以推断出,可以多次调用该方法以注册多个传感器,所有这些传感器都将在构建的Rugged实例中可用。在这里我们只调用了一次该方法,因此我们将只使用一个传感器。
在最后一个setter被调用之后,我们调用build()
方法来真正构建Rugged实例(而不是RuggedBuilder实例)。
Rugged默认考虑了一些更准确位置的修正:
不补偿延迟或速度组合主要用于针对不补偿的系统进行验证。当像素的视线已经包括光的畸变修正时,显然必须禁用该修正。
如果应该忽略这些修正,则在调用build()
之前必须插入其他一些setter:
setXxxx().
setLightTimeCorrection(false).
setAberrationOfLightCorrection(false).
build();
各种setter可以以任何顺序调用。唯一重要的是,一旦通过调用build()
方法创建了一个Rugged实例,它就与构建器无关,因此后续对setter的调用不会更改构建实例。实际上,可以创建一个构建器,然后调用其build()
方法创建第一个Rugged实例,然后通过再次调用一些setter来修改构建器配置,并从新配置中构建第二个Rugged实例。这允许在同一个程序中对两个不同的配置进行比较,而无需重新创建所有内容。例如,可以按照以下三个步骤进行:
RuggedBuilder ruggedBuilder = new RuggedBuilder().
setAlgorithm(xxxx).
[ ... some setters ...]
setTrajectory(xxxx);
然后添加所需的传感器:
ruggedBuilder.addLineSensor(lineSensor);
并创建Rugged实例:
Rugged rugged = ruggedBuilder.build();
最后一切都准备就绪,可以开始真正的工作了。让我们尝试定位地球上的一个点,即左上角的点(第一行,第一个像素):
import org.orekit.bodies.GeodeticPoint;
Vector3D position = lineSensor.getPosition(); // 这将返回一个零向量,因为我们将传感器相对于卫星的位置设置为0。
AbsoluteDate firstLineDate = lineSensor.getDate(0);
Vector3D los = lineSensor.getLOS(firstLineDate, 0);
GeodeticPoint upLeftPoint = rugged.directLocation(firstLineDate, position, los);
行号可以为负值;相关日期必须属于在构建Rugged时给定的时间范围内,即setTimeSpan(acquisitionStartDate, acquisitionStopDate, ...)。否则将抛出RuggedException异常。
像素号必须在0和(传感器像素大小-1)之间;在我们的例子中为0到1999。否则将抛出ArrayIndexOutOfBoundsException异常。
源代码可在DirectLocation.java中找到(位于src/tutorials下的fr.cs.examples包中)