一套工程可直接实现、完全对应哨兵1号 GRD 产品、从像素(x,y)→经纬度(lon,lat)的标准解法 + 计算公式 + 实现步骤，欧空局Sentinel-1官方定的几何模型，你拿到姿态、轨道、产品XML/标注文件就能直接编码计算。

已知：GRD图像像素坐标(X,Y) + 卫星星历/姿态 + 产品辅助数据 → 求(WGS84经纬度/大地坐标)。

一、先明确：哨兵1号 GRDH 是什么几何投影

Sentinel-1 Level-1 GRD（含GRDH）是：

• Ground Range，地距投影

• 零多普勒几何（Zero-Doppler Geometry）

• 已经做过：多视、幅度检测、地距重采样、burst拼接、几何校正到WGS84椭球（无DEM，仅椭球校正）

它不是斜距(Range-Doppler)平面，而是地距-方位(GR-Az)均匀采样的图像，因此定位不能直接用SLC的R-D模型硬套，必须用GRD产品规定的地理定位模型：

基于零多普勒面 + 地距/方位时间映射 + 星历插值 + 椭球求交。

二、核心思路

给定GRD图像像素坐标 (x, y)（x=距离向像素，y=方位向像素），流程是：

1. 像素(x,y) → 映射到方位向时间 ta、地距向采样点 rg

1. 由ta 插值卫星位置 Sta 、速度 Vta

1. 由地距 rg  → 换算斜距 R （GRD是地距，必须转斜距）

1. 在零多普勒面上，联立：

• 斜距方程：|P-Sta|=R

• WGS84椭球方程

1. 联立求解得到地面点 大地坐标 (B,L,H) → lat, lon

P表示大地坐标点位置，R表示斜距矢量。

GRD产品里已经提供了大量关键系数（行时间、地距起始、地距间隔、近距斜距、曲率校正系数等），不用你从头建模。

图 1 SAR图像目标定位分析图

注：1 姿态角会影响速度V的方向，因此需要求R*V，R为卫星轨道坐标系到天线坐标系旋转矩阵。注意中间还有卫星本体坐标系。

三、符号与坐标系约定（与哨兵1号产品文档一致）

• 图像坐标系：

• x：Ground Range 距离向像素，从左到右增大

• y：Azimuth 方位向像素，从上到下增大

• 输出：WGS84椭球，经纬度（lon, lat），单位度

• 卫星状态：位置S=XYZECEF ，速度V=VXVYVZECEF ，均为地心地固ECEF系

• 零多普勒条件：目标P在卫星瞬时速度矢量垂直的平面上

四、完整可代码化计算步骤（最关键）

步骤1：从GRD产品中读取定位关键参数

这些都在 S1A_IW_GRDH_1SDV_xxx.xsd.xml 或产品的 annotation XML 里：

1）方位向（Azimuth）关键参数

• linesPerBurst：每个burst行数

• numBursts：burst总数

• azimuthTimeStart, azimuthTimeEnd：每个burst的起止时间

• azimuthPixelSpacing：方位向像元间距（或直接用行时间差）

• 每个像素y对应的绝对方位时间 ta

2）地距向（Ground Range）关键参数

• firstGroundRangeSample：第一个距离向像素对应的地距

• groundRangePixelSpacing：地距像元间距 Δrg

• slantRangeToFirstGroundRangeSample：第一个地距像素对应的斜距

• groundRangeToSlantRange 多项式/查找表：GRD产品自带的地距↔斜距转换（由地球曲率、入射角决定）

重点：GRD是地距等间隔采样，斜距不是等间隔，必须用产品提供的转换，不能自己简单除cosθ。

步骤2：像素(x,y) → 方位时间 tₐ

对GRD，每一行y对应唯一的方位向时间（零多普勒时间）。

2.1 找到像素y属于哪个burst

在对应burst内，线性插值得到该行的绝对时间：

ta(y) = taz_start+yin_burst ˙ Δtaz

其中：

• tazstart ：该burst第一行的方位时间

• Δtaz ：方位向脉冲重复时间对应的行采样间隔

哨兵1号GRD的方位时间是严格逐行线性的，直接线性插值即可，精度足够。

2.2 由 tₐ 插值卫星位置与速度

使用提供的星历点（orbit state vectors），用拉格朗日插值或三次样条，插值得到：

S(ta)=[XS(ta), YS(ta), ZS(ta)]ECEF

V(ta)=[VX(ta), VY(ta), VZ(ta)]ECEF

（速度必须同时刻插值，不能近似）

步骤3：像素x → 地距 r_g → 斜距 R

3.1 像素x → 地距

rg(x)=rg0+x⋅Δrg​

• rg0 ：x=0对应地距

• Δrg ：groundRangePixelSpacing

3.2 地距 r_g → 斜距 R（核心！GRD专用）

Sentinel-1 GRD产品在XML里提供了地距到斜距的映射，一般两种形式：

1. 分段线性表

1. 多项式：

Rrg=a0+a1rg+a2rg2+a3rg3

严禁自己用 R = r_g / cos(θ) 近似，误差会到米级~十米级，不符合产品几何。

拿到R，就是该像素在零多普勒时刻对应的斜距。

步骤4：联立零多普勒 + 斜距 + 椭球，解地面点P

这是SAR几何定位的标准闭式，GRD同样遵守零多普勒定位模型。

我们要求地面点P（ECEF）满足3个方程：

方程1：斜距方程

∥P−S(ta)∥=R

即：

(PX−XS)2+(PY−YS)2+(PZ−ZS)2=R2                       (1)

方程2：零多普勒方程

GRD所有像素都在零多普勒面上：

(P−S(ta))⋅V(ta)=0                                      (2)

方程3：WGS84椭球方程

(PX2+PY2)/a2+PZ2/b2 = 1                                （3）

其中：

• a = 6378137.0m （长半轴）

• b = 6356752.314245m （短半轴）

步骤5：求解地面点P（ECEF）

(1)(2)(3) 是3个方程3个未知数，闭式可解，不需要迭代（工程常用闭式解）。

简化求解思路（最常用、最稳定）

由方程(2)，P 位于过S、以V

为法向量的平面。在该平面内建立局部坐标系：

u=V(ta)(法向)

w=∥S(ta)∥S(ta)​(地心径向)

v=w×u

得到局部正交基[u,v,w]

。

将P

表示为：

P=S+αv+βw

代入零多普勒方程(2)自动满足，再代入(1)和(3)，得到一个一元二次方程，直接求根，取地表正根即可。

解出P=PXPYPZ

后，ECEF转经纬度，就得到你要的地理坐标。

步骤6：ECEF → 经纬度 (lon, lat)

使用标准WGS84大地坐标系转换：

1. 经度：

λ=arctan2(PY, PX)

1. 纬度（迭代法，3步收敛）：

h=0，N=a，

循环：

五、GRD与SLC定位的关键区别（你必须注意，否则必错）

1. SLC：斜距等间隔，直接用R-D模型

1. GRD：地距等间隔，斜距不等距，必须用产品自带的gr→sr转换

• 这是90%的人算GRD定位不准的根源

1. GRD已经做过多视、拼接、地距重采样，不能再用SLC的burst原始斜距时间。

1. GRD定位不使用卫星姿态（姿态几乎不影响），只需要轨道+星历

• 因为GRD已经重采样到零多普勒地距，姿态误差在产品生产时已补偿

• 你给的姿态角，在GRD定位中几乎用不上，这是和光学图像最大的不同！

重要结论：

哨兵1号GRD像素定位，主要依赖轨道（位置+速度）+地距/斜距映射，姿态角几乎不参与计算。

光学图像必须用姿态，SAR-GRD不用，这是行业常识。

六、总结

哨兵1号GRDH数据为零多普勒地距投影产品，像素(x,y)到地理坐标的解算流程为：

1. 由像素方位向y插值得到零多普勒时间，并插值卫星位置与速度；

1. 由距离向x结合产品参数得到地距，再通过产品内置转换得到斜距R；

1. 联立斜距方程、零多普勒方程、WGS84椭球方程，闭式解算地面点ECEF坐标；

1. 将ECEF坐标转换为WGS84经纬度地理坐标。

该过程以轨道为核心，姿态角影响可忽略，与光学图像基于共线方程与姿态的定位模型完全不同。