一套工程可直接实现、完全对应哨兵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方位向像素流程是像素(x,y) → 映射到方位向时间ta、地距向采样点rg由ta插值卫星位置Sta、速度Vta由地距rg→ 换算斜距RGRD是地距必须转斜距在零多普勒面上联立斜距方程|P-Sta|RWGS84椭球方程联立求解得到地面点大地坐标(B,L,H) → lat, lonP表示大地坐标点位置R表示斜距矢量。GRD产品里已经提供了大量关键系数行时间、地距起始、地距间隔、近距斜距、曲率校正系数等不用你从头建模。图 1 SAR图像目标定位分析图注1姿态角会影响速度V的方向因此需要求R*VR为卫星轨道坐标系到天线坐标系旋转矩阵。注意中间还有卫星本体坐标系。三、符号与坐标系约定与哨兵1号产品文档一致图像坐标系xGround Range距离向像素从左到右增大yAzimuth方位向像素从上到下增大输出WGS84椭球经纬度lon, lat单位度卫星状态位置SXYZECEF速度VVXVYVZECEF均为地心地固ECEF系零多普勒条件目标P在卫星瞬时速度矢量垂直的平面上四、完整可代码化计算步骤最关键步骤1从GRD产品中读取定位关键参数这些都在S1A_IW_GRDH_1SDV_xxx.xsd.xml或产品的annotationXML 里1方位向Azimuth关键参数linesPerBurst每个burst行数numBurstsburst总数azimuthTimeStart,azimuthTimeEnd每个burst的起止时间azimuthPixelSpacing方位向像元间距或直接用行时间差每个像素y对应的绝对方位时间ta2地距向Ground Range关键参数firstGroundRangeSample第一个距离向像素对应的地距groundRangePixelSpacing地距像元间距ΔrgslantRangeToFirstGroundRangeSample第一个地距像素对应的斜距groundRangeToSlantRange多项式/查找表GRD产品自带的地距↔斜距转换由地球曲率、入射角决定重点GRD是地距等间隔采样斜距不是等间隔必须用产品提供的转换不能自己简单除cosθ。步骤2像素(x,y) →方位时间tₐ对GRD每一行y对应唯一的方位向时间零多普勒时间。2.1找到像素y属于哪个burst在对应burst内线性插值得到该行的绝对时间ta(y) taz_startyin_burst˙Δtaz其中tazstart该burst第一行的方位时间Δtaz方位向脉冲重复时间对应的行采样间隔哨兵1号GRD的方位时间是严格逐行线性的直接线性插值即可精度足够。2.2由tₐ插值卫星位置与速度使用提供的星历点orbit state vectors用拉格朗日插值或三次样条插值得到S(ta)[XS(ta), YS(ta), ZS(ta)]ECEFV(ta)[VX(ta), VY(ta), VZ(ta)]ECEF速度必须同时刻插值不能近似步骤3像素x →地距r_g →斜距R3.1像素x →地距rg(x)rg0x⋅Δrgrg0x0对应地距ΔrggroundRangePixelSpacing3.2地距r_g →斜距R核心GRD专用Sentinel-1 GRD产品在XML里提供了地距到斜距的映射一般两种形式分段线性表多项式Rrga0a1rga2rg2a3rg3严禁自己用R r_g / cos(θ)近似误差会到米级~十米级不符合产品几何。拿到R就是该像素在零多普勒时刻对应的斜距。步骤4联立零多普勒斜距椭球解地面点P这是SAR几何定位的标准闭式GRD同样遵守零多普勒定位模型。我们要求地面点PECEF满足3个方程方程1斜距方程∥P−S(ta)∥R即(PX−XS)2(PY−YS)2(PZ−ZS)2R2 (1)方程2零多普勒方程GRD所有像素都在零多普勒面上(P−S(ta))⋅V(ta)0 (2)方程3WGS84椭球方程(PX2PY2)/a2PZ2/b2 13其中a 6378137.0m 长半轴b 6356752.314245m 短半轴步骤5求解地面点PECEF(1)(2)(3) 是3个方程3个未知数闭式可解不需要迭代工程常用闭式解。简化求解思路最常用、最稳定由方程(2)P位于过S、以V为法向量的平面。在该平面内建立局部坐标系uV(ta)(法向)w∥S(ta)∥S(ta)(地心径向)vw×u得到局部正交基[u,v,w]。将P表示为PSαvβw代入零多普勒方程(2)自动满足再代入(1)和(3)得到一个一元二次方程直接求根取地表正根即可。解出PPXPYPZ后ECEF转经纬度就得到你要的地理坐标。步骤6ECEF →经纬度(lon, lat)使用标准WGS84大地坐标系转换经度λarctan2(PY, PX)纬度迭代法3步收敛h0Na循环五、GRD与SLC定位的关键区别你必须注意否则必错SLC斜距等间隔直接用R-D模型GRD地距等间隔斜距不等距必须用产品自带的gr→sr转换这是90%的人算GRD定位不准的根源GRD已经做过多视、拼接、地距重采样不能再用SLC的burst原始斜距时间。GRD定位不使用卫星姿态姿态几乎不影响只需要轨道星历因为GRD已经重采样到零多普勒地距姿态误差在产品生产时已补偿你给的姿态角在GRD定位中几乎用不上这是和光学图像最大的不同重要结论哨兵1号GRD像素定位主要依赖轨道位置速度地距/斜距映射姿态角几乎不参与计算。光学图像必须用姿态SAR-GRD不用这是行业常识。六、总结哨兵1号GRDH数据为零多普勒地距投影产品像素(x,y)到地理坐标的解算流程为由像素方位向y插值得到零多普勒时间并插值卫星位置与速度由距离向x结合产品参数得到地距再通过产品内置转换得到斜距R联立斜距方程、零多普勒方程、WGS84椭球方程闭式解算地面点ECEF坐标将ECEF坐标转换为WGS84经纬度地理坐标。该过程以轨道为核心姿态角影响可忽略与光学图像基于共线方程与姿态的定位模型完全不同。
Sentinel-1 GRD 几何定位计算详解
一套工程可直接实现、完全对应哨兵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方位向像素流程是像素(x,y) → 映射到方位向时间ta、地距向采样点rg由ta插值卫星位置Sta、速度Vta由地距rg→ 换算斜距RGRD是地距必须转斜距在零多普勒面上联立斜距方程|P-Sta|RWGS84椭球方程联立求解得到地面点大地坐标(B,L,H) → lat, lonP表示大地坐标点位置R表示斜距矢量。GRD产品里已经提供了大量关键系数行时间、地距起始、地距间隔、近距斜距、曲率校正系数等不用你从头建模。图 1 SAR图像目标定位分析图注1姿态角会影响速度V的方向因此需要求R*VR为卫星轨道坐标系到天线坐标系旋转矩阵。注意中间还有卫星本体坐标系。三、符号与坐标系约定与哨兵1号产品文档一致图像坐标系xGround Range距离向像素从左到右增大yAzimuth方位向像素从上到下增大输出WGS84椭球经纬度lon, lat单位度卫星状态位置SXYZECEF速度VVXVYVZECEF均为地心地固ECEF系零多普勒条件目标P在卫星瞬时速度矢量垂直的平面上四、完整可代码化计算步骤最关键步骤1从GRD产品中读取定位关键参数这些都在S1A_IW_GRDH_1SDV_xxx.xsd.xml或产品的annotationXML 里1方位向Azimuth关键参数linesPerBurst每个burst行数numBurstsburst总数azimuthTimeStart,azimuthTimeEnd每个burst的起止时间azimuthPixelSpacing方位向像元间距或直接用行时间差每个像素y对应的绝对方位时间ta2地距向Ground Range关键参数firstGroundRangeSample第一个距离向像素对应的地距groundRangePixelSpacing地距像元间距ΔrgslantRangeToFirstGroundRangeSample第一个地距像素对应的斜距groundRangeToSlantRange多项式/查找表GRD产品自带的地距↔斜距转换由地球曲率、入射角决定重点GRD是地距等间隔采样斜距不是等间隔必须用产品提供的转换不能自己简单除cosθ。步骤2像素(x,y) →方位时间tₐ对GRD每一行y对应唯一的方位向时间零多普勒时间。2.1找到像素y属于哪个burst在对应burst内线性插值得到该行的绝对时间ta(y) taz_startyin_burst˙Δtaz其中tazstart该burst第一行的方位时间Δtaz方位向脉冲重复时间对应的行采样间隔哨兵1号GRD的方位时间是严格逐行线性的直接线性插值即可精度足够。2.2由tₐ插值卫星位置与速度使用提供的星历点orbit state vectors用拉格朗日插值或三次样条插值得到S(ta)[XS(ta), YS(ta), ZS(ta)]ECEFV(ta)[VX(ta), VY(ta), VZ(ta)]ECEF速度必须同时刻插值不能近似步骤3像素x →地距r_g →斜距R3.1像素x →地距rg(x)rg0x⋅Δrgrg0x0对应地距ΔrggroundRangePixelSpacing3.2地距r_g →斜距R核心GRD专用Sentinel-1 GRD产品在XML里提供了地距到斜距的映射一般两种形式分段线性表多项式Rrga0a1rga2rg2a3rg3严禁自己用R r_g / cos(θ)近似误差会到米级~十米级不符合产品几何。拿到R就是该像素在零多普勒时刻对应的斜距。步骤4联立零多普勒斜距椭球解地面点P这是SAR几何定位的标准闭式GRD同样遵守零多普勒定位模型。我们要求地面点PECEF满足3个方程方程1斜距方程∥P−S(ta)∥R即(PX−XS)2(PY−YS)2(PZ−ZS)2R2 (1)方程2零多普勒方程GRD所有像素都在零多普勒面上(P−S(ta))⋅V(ta)0 (2)方程3WGS84椭球方程(PX2PY2)/a2PZ2/b2 13其中a 6378137.0m 长半轴b 6356752.314245m 短半轴步骤5求解地面点PECEF(1)(2)(3) 是3个方程3个未知数闭式可解不需要迭代工程常用闭式解。简化求解思路最常用、最稳定由方程(2)P位于过S、以V为法向量的平面。在该平面内建立局部坐标系uV(ta)(法向)w∥S(ta)∥S(ta)(地心径向)vw×u得到局部正交基[u,v,w]。将P表示为PSαvβw代入零多普勒方程(2)自动满足再代入(1)和(3)得到一个一元二次方程直接求根取地表正根即可。解出PPXPYPZ后ECEF转经纬度就得到你要的地理坐标。步骤6ECEF →经纬度(lon, lat)使用标准WGS84大地坐标系转换经度λarctan2(PY, PX)纬度迭代法3步收敛h0Na循环五、GRD与SLC定位的关键区别你必须注意否则必错SLC斜距等间隔直接用R-D模型GRD地距等间隔斜距不等距必须用产品自带的gr→sr转换这是90%的人算GRD定位不准的根源GRD已经做过多视、拼接、地距重采样不能再用SLC的burst原始斜距时间。GRD定位不使用卫星姿态姿态几乎不影响只需要轨道星历因为GRD已经重采样到零多普勒地距姿态误差在产品生产时已补偿你给的姿态角在GRD定位中几乎用不上这是和光学图像最大的不同重要结论哨兵1号GRD像素定位主要依赖轨道位置速度地距/斜距映射姿态角几乎不参与计算。光学图像必须用姿态SAR-GRD不用这是行业常识。六、总结哨兵1号GRDH数据为零多普勒地距投影产品像素(x,y)到地理坐标的解算流程为由像素方位向y插值得到零多普勒时间并插值卫星位置与速度由距离向x结合产品参数得到地距再通过产品内置转换得到斜距R联立斜距方程、零多普勒方程、WGS84椭球方程闭式解算地面点ECEF坐标将ECEF坐标转换为WGS84经纬度地理坐标。该过程以轨道为核心姿态角影响可忽略与光学图像基于共线方程与姿态的定位模型完全不同。
