2． 信息工程大学 地理空间信息学院，河南 郑州 450052
2．Institute of Geospatial Information，Information Engineering University，Zhengzhou 450052,China
1 引 言
空间地图投影是专门为各种卫星遥感图像建立的动态地图投影，是卫星遥感制图的动态数学基础。侧视雷达图像的精确处理与制图也需要建立相应的空间地图投影。
自从1974年文献[1]提出空间斜圆柱投影(SOM)概念以来，国内外很多学者涉足该领域的研究。1977年，文献[2]和文献[3]各自独立地导出了SOM投影公式。从20世纪80年代开始，文献[4, 11, 12, 13]围绕SOM投影进行系列研究；文献[5, 8]对等角空间投影进行了研究；文献[6, 7, 17]研究了空间投影的系统理论，并探讨了一些应用；文献[16]研究了等角投影之间变换的微分方法。
侧视雷达是一种主动的微波遥感成像系统，通过向地面发射脉冲波获得数据胶片。制图工作的主要问题是如何记录这些数据胶片的数学模型，并进行几何校正，使得成图误差较小。图像中所包含的几何畸变可以表示为图像上各像元的位置坐标与地图坐标系中的目标地物坐标的差异。SAR图像的几何畸变有诸多因素引起，其中包括卫星姿态、地球自转和弯曲、轨道进动以及图像投影面的选取等。要对侧视雷达图像进行高精度的几何纠正，选择合适的地图投影是至关重要的。
现有的适合侧视雷达图像的数学模型是逐点或逐行建立的瞬时构像方程[9]，不能连续地表示所采集的数据。SOM投影和等角空间投影(CSP)等所研究出来的一系列空间地图投影都是专门为沿星下点轨迹周围扫描成像(如SPOT和TM图像等)的卫星图像而设计的空间地图投影，投影精度只能满足星下点轨迹周围经纬度±1°以内的扫描图像的精确制图。而侧视雷达的侧视区域中心线一般偏离星下点轨迹大约270 km以上，即侧视区域在星下点轨迹周围±1°之外，比如海洋卫星侧视雷达成像系统的侧视区域偏离星下点轨迹±2°之外[10]，显然用SOM投影或等角空间投影等用来表示侧视雷达图像是不合适的，投影精度不能满足制图要求，因此，专门为侧视雷达图像设计一种空间地图投影是非常必要的。
本文目的是模拟侧视雷达成像的物理过程和几何关系，建立一种适合星载侧视雷达图像的空间斜圆锥投影。该投影保持侧视区域中心线不变形，直接建立侧视区域内像素点与地面点之间较近似的对应关系。该投影是为精确处理侧视雷达图像而专门建立的一种动态地图投影。
2 空间斜圆锥投影的几何原理考虑卫星运动、地球自转和轨道进动等因素，假设有一个空间斜圆锥切于地球，如图 1所示，切线为侧视区域中央曲线，该圆锥与卫星和地球之间存在着相对运动关系，按此几何模型来模拟侧视雷达成像的物理过程，从而导出一种投影模型，该投影使SAR图像在一个空间斜圆锥上构成，并且保持侧视区域中心线长度不变形和曲率不变。将侧视区域投影在圆锥上，然后再将圆锥面展成平面，即成空间斜圆锥投影。
|
| 图 1 空间斜圆柱投影几何模型 Fig. 1 Geometric model of space oblique conic projection |
为使空间斜圆锥总是与地球保持瞬间相切，可使圆锥表面绕其自身轴摆动。尽管每条轨道所对应的地面景物变动的经度范围接近26°，但每条侧视区域中心线在投影平面上的投影是固定不变的。因此，即使每一条轨道所侧视的是地球表面的不同部分的图像，但侧视区域条带的投影坐标系是始终保持连续不变的。这表示侧视区域内点的投影坐标是地面点经纬度φ、λ以及时间t的函数，且是t的周期函数，即x=f1(φ,λ,t)，y=f2(φ,λ,t)。
3 侧视区域中心线轨迹投影为了建立空间斜圆锥投影公式，首先要建立侧视区域中心线L在投影面上的投影公式，一般是基于两个条件，即保持侧视区域中心线L长度不变形和曲率不变，如图 2所示。
|
| 图 2 星下点轨迹投影 Fig. 2 Projection of ground track line |
设中心线L上点A(φL,λL,t)投影到地图平面上为(xL,yL,t)。
3.1 等长条件设侧视区域中心线L的弧长为
投影面上中心线L的弧长为 式中，vL(t)是t时刻侧视区域中心线L上的扫描点的瞬时速率。为使长度保持不变，即s=s′，则 3.2 曲率不变条件
为了保形，要求中心线L的投影轨迹的曲率半径等于其在椭球切平面内投影轨迹的瞬时曲率半径ρ(t)，即
3.3 解出星下点的投影轨迹方程
联立式(3)、式(4)解出星下点的投影轨迹方程
由此推得侧视区域中心线L的投影为 4 空间斜圆锥投影公式
根据前面所建立的侧视区域中心线L的投影公式，详细推导出空间斜圆锥投影公式。方法是将侧视区域内的点按正形圆锥投影公式进行投影，再考虑侧视点的动态性，最后得到相应的空间斜圆锥投影公式。
因为卫星运行轨道与赤道是斜交的，为方便地运用正形圆锥投影公式，首先需要进行赤道变换。假设星下点轨迹线是新赤道，卫星运行轨道相对赤道的倾角为i，如图 3所示考虑到地球自转，根据文献[11, 18]，若已知旧经纬度(φ,λ)，求新经纬度(φ′,λ′)的公式为
根据式(7)，运用迭代法即可由旧经纬度求出新经纬度。
|
| 图 3 变换赤道 Fig. 3 Transformed equator |
反之，若已知新经纬度(φ′,λ′)，求旧经纬度(φ,λ)的公式为
如图 4所示，考虑侧视区域内的点在投影面上的投影，因为侧视区域是一个不宽的带形，侧视带内的点可以看做是中心线上点的很小位移，所以点到中心线的距离可以近似地看做是直线距离。
|
| 图 4 空间斜圆锥投影图 Fig. 4 Figure of space oblique conic projection |
假设零多普勒条件成立。如图 4所示，卫星运行到t时刻，与其侧视区域中心线相切的正轴等角圆锥投影几何关系见图 5。此时卫星飞行方向与x轴的夹角为
|
| 图 5 正轴等角圆锥投影 Fig. 5 Conformal conic projection |
在图 5中，把瞬时直角坐标系x1Ay1旋转一个ε角，即得
根据文献[20]，以侧视区域中心线上的点为原点，卫星飞行方向为x轴建立坐标系，如图 5所示，设B(φ′,λ′)是侧视区域内的点，则相对于新赤道的等角正轴圆锥投影公式为
式中
(R为地球半径；s0为侧视区域中心线到星下点轨迹的距离)
把直角坐标系x2Ay2与直角坐标系xOy进行叠加，对于侧视区域内的任一点B，它在直角坐标系xOy中的投影坐标为
式中，t为从初始时刻开始沿轨道运行的时间；vL(t)为t时刻侧视区域中心线上的点A相对地球的瞬时速度；ρ(t)为A点处的瞬时曲率半径。 5 侧视点B的经纬度φ(t),λ(t)与时间t的关系
假定卫星运行t时刻侧视到B(φ,λ)点(即波束指向B(φ,λ))，此时相应的侧视区域中心线上点B0的地理坐标为(φ0(t),λ0(t))，确定φ、λ与时间t的关系。图 6中，根据零多普勒条件，B(φ,λ)和B0(φ0,λ0)都应处在由单位法向量n和脉冲向量w所决定的平面内，即满足恒等式
式中,ΔR为从B0(φ0,λ0)到B(φ,λ)的向量
ΔR=R－R0= (Rx－Rx0)i+(Ry－Ry0)j+(Rz－Rz0)k
|
| 图 6 卫星侧视点位移 Fig. 6 Displacement of satellite side glance point |
根据文献[3]，由零多普勒条件知，卫星运行t时刻后，λ′=λ′0=2πt/P2，φ′0=s0/R为已知常数。把φ′0、λ′0和φ′、λ′分别代入式(7)，通过迭代即可求得φ0、λ0和、λ。
6 反解变换根据式(13)和式(14)，对任意时间t，取初始经纬度为φ(0)=φL(t)，λ(0)=λL(t)，代入式(13)和式(14)，就得到相应的SOC坐标x(0)、y(0)。构造迭代公式
式中，φ(k)、λ(k)代入到式(13)和式(14)可得到相应的SOC坐标x(k)、y(k)；偏导数矩阵也可由式(13)和式(14)求得。
根据式(20)，经过有限次的迭代(一般需要迭代4到5步)就可得到相应的反解经纬度φ(t)、λ(t)。
7 投影变形先求式(13)和式(14)的各个偏导数。分别对这两个式子求偏导数，先求对φ的偏导数
同理，关于λ或t的偏导数就是把上两式中的αφ换成αλ或αt既可，今后用φ→λ,t记之。
下面讨论投影变形公式，首先讨论经纬线长度变形，根据文献[6]，经线变形比为
纬线变形比为
式中
其次再讨论角度长度变形，投影后经纬线夹角为
式中
方位角投影后的变形公式
式中
S=E1r2+F1Mrtanα+2F2Mrtanα(dt/dλ)+ F3M2tan2α(dt/dλ)+E3M2tan2α(dt/dλ)2
8 应用举例以海洋卫星侧视雷达成像系统为例，通过模拟数据进行计算，求SOC正反解坐标。参考文献[10]，假设海洋卫星侧视雷达成像系统雷达图像测量的有关参数为：卫星运行周期为P2=100min；卫星轨道为圆形；卫星轨道面倾角i=108°(逆算)；卫星轨道近升距ω=0.0；升交点赤经Ω=0.0；卫星飞行高度为H=800km；侧视条带宽度为D=100km(见图 7)；视线高度角为α1=16.9°与α2=23.1°之间(见图 8)；侧视区域中心线离星下点的距离s0=292km。并假设卫星运行初始时刻为降交点处，取t0=τ=0，θ0=0；另外选取地球半径R=6378160.0m。
|
| 图 7 卫星扫描示意图 Fig. 7 Sketch of satellite scan |
|
| 图 8 视线高度角 Fig. 8 Angle of view |
根据上述数据可求得侧视偏天底角α=s0/H=20.9°，卫星运动的平均角速度为n=2π/P2=3.6°=0.062832rad/min，侧视区域中心线的变换纬度为φ′0=s0/R=2.623°，t时刻卫星运动的偏近点角E=nt，地球自转角θ=ωet。
考虑地球自转，可求得卫星在地心坐标系中的坐标为
经计算可得
把上述结果代入式(12)和(13)得
式中，ε、d、d0、δ的表达式见式(12)和式(15)～(18)。 9 结 论
为满足卫星侧视雷达图像制图的需要，本文对卫星侧视雷达图像的空间地图投影问题进行了研究，得到了适合卫星侧视雷达图像数据表述的空间斜圆锥投影。
本文研究结果的价值主要体现在：
(1) 本文在前人研究空间投影理论的基础上，提出了探求空间地图投影的一种方法。前人建立空间地图投影一般都是考虑以切平面为基准的空间方位投影和以圆柱为基准的空间圆柱投影。本文以空间切圆锥面为基准，建立相应的空间斜圆锥投影。
(2) 导出了空间斜圆锥投影公式，得出了一种空间地图投影模型，发展了空间地图投影理论。与空间斜墨卡托投影相比，在几何模型选择上有改进，而且适宜于表述的遥感影像种类也不同。空间斜墨卡托投影适宜于表述垂直摄影或扫描所获取的遥感影像数据，空间斜圆锥投影适宜于表述侧视雷达所获取的遥感影像数据。
(3) 为卫星侧视雷达图像数据建立了比较适宜的动态数学基础，数据处理误差要比选用空间斜墨卡托投影作为数学基础要小得多，对侧视雷达图像数据的精确处理和连续遥感制图具有一定的使用价值。
| [1] | COLVOCORESSES A P. Space Oblique Mercator[J]. Pho-togrammetric Engineering, 1974, 40(8): 921-926. |
| [2] | JUNKINS J L,TURNER J D. Formulation of a Space Oblique Mercator Map Projection[R]. Charlottesville:University of Virginia, 1977:1-80. |
| [3] | SNYDER J P. Space Oblique Mercator Projection Mathematical Development[J]. Geological Survey Bulletin, 1977,1518:1-90. |
| [4] | YANG Qihe, JOHN P,WALDO R,et al. Map Projection Transformation—Principles and Applications[M]. London:Taylor and Francis, 1999:49-50. |
| [5] | CHENG Yang. The Conformal Space Projection[J]. Cartography and Geographic Information Systems,1996, 23(1):37-50. |
| [6] | REN Liucheng. Space Projection Theory and Its Applications in Remote Sensing Technology[M]. Beijing:Science Press, 2003:7-14.(任留成. 空间投影理论及其在遥感技术中的应用[M]. 北京:科学出版社,2003:7-14.) |
| [7] | REN L C, KEITH C C, ZHOU C H, et al. Geometric Rectification of Satellite Imagery with Minimal Ground Control Using Space Oblique Mercator Projection Theory[J]. Cartography and Geographic Information Science, 2010, 37(4):261-272. |
| [8] | CHENG Yang. On Conformal Projection Maintaining a Desired Curve on the Ellipsoid without Distortion[R]. Washington: ASPRS/ACSM,1992:294-303. |
| [9] | QIAN Zengbo, LIU Jingyu, XIAO Guochao. Space Photogrammetry[M].Beijing:PLA Press,1992:60-107.(钱曾波,刘静宇,肖国超. 航天摄影测量[M].北京:解放军出版社,1992:60-107.) |
| [10] | FRANZ L. Radar Imaging Applied to Surveying and Mapping of Map and Chart[M]. Proceedings of Remote Sensing.Beijing: Surveying and Mapping Press,1981:50-55.(佛朗茨.成像雷达用于测制地图和海图 [M]. 遥感译文集.北京:测绘出版社,1981:50-55.) |
| [11] | YANG Xiaomei,YANG Qihe, YANG Yubin.Three Types of Latitudinal Function and Linear Interpolation Algorithms Used in Surveying and Cartography[J]. Acta Geodetica et Cartographica Sinica, 1997, 26(1):92-95.(杨晓梅,杨启和,杨玉彬. 测量与制图中运用的三种纬度函数和线性插值算法[J].测绘学报,1997, 26(1):92-95.) |
| [12] | YANG Qihe. Perspective Azimuthal Projection of Changing View Points[J].Acta Geodetica et Cartographica Sinica,1987,16(4):105-110.(杨启和.变视点透视方位投影[J].测绘学报,1987, 16(4):105-110.) |
| [13] | YANG Qihe. Search for the Affine Transformation of Equal-area Projection and Its Improved Projections[J].Acta Geodetica et Cartographica Sinica, 1992, 21(1):46-52.(杨启和. 等面积投影的仿射变换及其改良的等面积投影探求[J].测绘学报,1992,21(1):46-52.) |
| [14] | WALDO R T. Polycylindric Map Projections[J].The American Cartographer,1986,13(2): 117-120. |
| [15] | WALDO R T. Bidimensional Regression[J].Geographical Analysis, 1994, 26(2): 186-212. |
| [16] | LI Yongpu.Differential Method of the Transformation between Conformal Projections[J]. Acta Geodetica et Cartographica Sinica, 1994, 23(1):35-42.(李永朴.平面坐标系统转换的微分方程[J].测绘学报,1994,23(1):35-42.) |
| [17] | REN Liucheng,YE Jiankao, ZHAO Qi. Research on the Distortion Theory of Map Projection by Operator Differential[J].Acta Geodetica et Cartographica Sinica, 1999, 28(1):76-80.(任留成,叶建栲,赵琪. 利用算子微分研究投影变形理论[J].测绘学报,1999,28(1):76-80.) |
| [18] | YANG Qihe. Principle and Method of the Transformation of Map Projection[M].Beijing:PLA Press,1989:20-35.(杨启和. 地图投影变换原理与方法[M].北京:解放军出版社,1989:20-35.) |
| [19] | LI Guozao, YANG Qihe, HU Dingquan. Map Projection[M].Beijing:PLA Press,1993:50-70.(李国藻,杨启和,胡定全. 地图投影[M]. 北京:解放军出版社,1993:50-70.) |
| [20] | WU Zhongxing, YANG Qihe. Principle of Mathematic Cartography[M].Beijing:Surveying and Mapping Press,1989:30-60.(吴忠性,杨启和. 数学制图学原理[M]. 北京:测绘出版社,1989:30-60.) |