SAR入射角与方位角全解：元数据、轨道解算与逐像元栅格

做地形归一化、做后向散射角度改正、算局部入射角做土壤水分反演，这几件事里第一个卡住人的地方，往往不是算法，而是"这两个角度到底从哪来"。我第一次处理 Sentinel-1A 的 IW 数据时，想当然地以为入射角会像 Landsat 的太阳角度那样，明明白白写在元数据的某个字段里，结果翻遍了 manifest、翻遍了 annotation，只找到几十个点组成的网格和一个 incidenceAngleMidSwath ，前者稀得没法直接用，后者只是整幅影像的中点粗值。后来才慢慢摸清楚：Sentinel-1 的角度信息有三条获取路径，精度、成本、适用场景完全不同，选错路径会让后面的工作全白做。这篇就把入射角和方位角这两件事从头拆开，包括它们各自的定义陷阱、从注释文件里能直接抠出什么、怎么用轨道状态矢量自己算，以及逐像元出栅格的几条工程路线。适合已经能打开 S1 影像、但卡在角度这一环的读者，也适合被"方位角"这个词搞晕过一次的人。

1. 三个名字相近的角度，定义差了十万八千里

1.1 入射角：椭球法线和真实地表法线，是两个世界

先把最容易混的一个概念摆正。SAR 影像里说的"入射角"（incidence angle），标准定义是 雷达视线方向与当地地表法线的夹角 。注意这里的关键词是"当地地表法线"。如果把地表当成一个光滑的地球椭球，那法线就是椭球法线，算出来的角度叫 椭球入射角 （ellipsoid incidence angle）；如果考虑真实地形，法线换成 DEM 起伏面的法线，算出来的叫 本地入射角 （local incidence angle）。

这两个值的差别在平原上可能只有零点几度，在山地能拉开到 ±15° 甚至更多。原因很直观：斜坡朝向卫星时，视线与坡面法线的夹角变小，本地入射角小于椭球入射角；斜坡背向卫星时反之。做过辐射地形校正的人都知道，山区影像一明一暗的条纹，本质上就是本地入射角在剧烈变化。

注意：很多人拿椭球入射角去做地形归一化，做完发现山区还是花的，问题不在算法，在用错了角度输入。

还有一个更细的分支叫 投影本地入射角 （projected local incidence angle），它是在地图投影坐标系下测量视线与地表法线的夹角，会额外引入投影变形带来的偏差。做辐射改正时一般用非投影的本地入射角更干净，做统计分析和跟其他投影数据比对时，投影版本更好对齐。

至于 仰角 （elevation angle），它和入射角是互补关系：仰角 = 90° − 入射角，指的是视线与当地水平面的夹角。有些文献里两个词混着用，看到数值先判断一下量级——入射角在 Sentinel-1 的 IW 模式下大致落在 29°~46°，仰角则在 44°~61°，靠取值范围就能反推对方用的是哪套定义。

1.2 方位角：一个词，两种完全不同的东西

方位角这个词是真正的重灾区。在 Sentinel-1 的语境里，它至少有两个含义：

第一是 平台航向角 （heading angle），描述卫星速度矢量在水平面上的投影与正北方向的夹角，顺时针为正。这个值跟成像几何无关，只跟轨道有关，在整景影像里基本是常数（同一轨道弧段内变化很小）。

第二是 视线方位角 （line-of-sight azimuth），描述从地面目标指向卫星的方向，投影到水平面后与正北的夹角，同样顺时针为正。这个值才是做几何分析、形变分解、坡向匹配时真正需要的。

两者在零多普勒、正侧视的理想情况下，差值正好是 90°。Sentinel-1 是右视成像，天线指向卫星前进方向的右侧，所以视线方位角 ≈ 航向角 + 90°。反过来，如果你在某些工具里看到的角度跟航向角基本一致，那它给的一定是 heading，不是 LOS azimuth，别拿去算坡向。

这个差别在实际工作里非常要命。做 InSAR 形变分解成东西向和南北向分量时，用的是视线方位角；做叠掩和阴影分析时，需要的是视线矢量在水平面上的投影方向。拿错一个，结果会整体偏 90°，而且由于对称性，图上看起来还挺"正常"，很难一眼发现。

1.3 一张表把定义钉死

名称	几何含义	Sentinel-1 IW 典型取值	主要用途
椭球入射角	视线与椭球法线夹角	29°~46°（随距离向变化）	辐射定标、后向散射角改正的基准
本地入射角	视线与 DEM 坡面法线夹角	椭球值 ±15°，山地波动大	地形辐射校正、分类特征
投影本地入射角	投影坐标系下的本地入射角	与本地值接近，含投影变形	与其他投影数据联合分析
仰角	视线与水平面夹角	44°~61°	与入射角互补，注意别搞反
平台航向角	速度矢量投影与正北夹角	升轨约 60°~100°，降轨约 240°~280°	轨道分析、几何粗校验
视线方位角	目标指向卫星的水平方向	航向角 + 90°（右视）	形变分解、坡向匹配、阴影分析

表里航向角那一行的范围之所以写得宽，是因为它随纬度变化：轨道倾角约 98°，在赤道附近航向接近东西向，越往高纬，速度方向的南北分量越大，航向角会明显偏离东西向。所以只能用这个范围做粗筛，不能当精确值用。

2. 注释文件里到底埋了什么，够不够用

2.1 geolocationGrid：几十个点的稀疏真相

解开一个 S1 产品的压缩包，annotation 目录下的 XML 里藏着一块 geolocationGrid 。用任意文本编辑器搜这个关键词，你会看到一串 geolocationGridPoint ，每个点带着这些字段：

• azimuthTime ：该点的方位向时间
• slantRangeTime ：该点的斜距时间
• latitude / longitude / height ：该点对应的地面位置
• incidenceAngle ：入射角
• elevationAngle ：仰角
• azimuthAngle ：方位角（这里给的是航向角口径）

关键问题是密度。对于 IW 模式的 SLC 产品，这一组网格点在距离向大约每 100 个像元一个、方位向大约每 10 个 burst 边界一个，整景算下来也就一两百个点。用这么稀疏的采样去描述一个随距离向单调变化、随方位向缓慢变化的角度场，理论上插值是可行的，因为入射角在方位向的变化确实很平缓（同一轨道弧段内几十公里才变零点几度）。

但也有个致命细节： 这些点的 azimuthAngle 是航向角，不是视线方位角 。官方文档对这个字段的描述在不同版本的产品说明里措辞有过调整，不同第三方库的解读也不完全一致。我的做法是：拿到这个字段后，先跟同一景的状态矢量算出来的航向角做比对，如果两者吻合到 0.1° 以内，就确认它是 heading；如果不吻合，再考虑是不是做了某种转换。这一步校验花不了五分钟，但能避免后面全盘错位。

2.2 一个字段解决不了整幅影像

imageAnnotation 下面还有一个 incidenceAngleMidSwath ，顾名思义是测绘带中点的入射角。这个值拿来写论文的方法章节、做数据描述没问题，但它只是一个常数，反映不了距离向从近端到远端十几度的跨度。有人图省事直接拿它当整景入射角去做辐射改正，结果近距端和远距端的后向散射偏差能到 1 dB 以上，做时间序列分析时会出现虚假的趋势。

正确的用法是把它当成一个 sanity check：你自己算出来的入射角场，中间那一行的均值应该跟它对得上，对不上说明哪里出了问题。

2.3 用插值把这些点用起来

如果你的需求只是"每个像元有个大致合理的入射角"，插值 geolocationGrid 是最省事的路线。思路是把网格点的 (slantRangeTime, azimuthTime) 作为自变量、 incidenceAngle 作为因变量，构建一个二维的规则网格，然后对任意像元坐标做双线性插值。

import xml.etree.ElementTree as ET
import numpy as np

def read_geolocation_grid(xml_path):
    """从 S1 annotation XML 中读取地理定位网格点"""
    root = ET.parse(xml_path).getroot()
    fields = ['azimuthTime', 'slantRangeTime', 'latitude', 'longitude',
              'height', 'incidenceAngle', 'elevationAngle', 'azimuthAngle']
    data = {f: [] for f in fields}
    for pt in root.findall('.//geolocationGridPoint'):
        for f in fields:
            txt = pt.findtext(f)
            data[f].append(txt)
    # 时间和角度转成数值
    out = {}
    for f in fields:
        if f == 'azimuthTime':
            out[f] = np.array([np.datetime64(s.replace('Z', '')) for s in data[f]])
        else:
            out[f] = np.array([float(s) for s in data[f]])
    return out

拿到这些点之后，还要注意两个工程细节。一是网格点的 azimuthTime 是 UTC 字符串，格式类似 2023-05-01T10:23:45.123456 ，转成 datetime64 后要做成单调递增的数组，有些产品在 burst 边界处会出现重复时间戳，插值前必须去重。二是 slantRangeTime 在 GRD 产品里没有意义（GRD 已经投影到地距），这时候要用 latitude/longitude 作为网格坐标来做插值，或者干脆跳过插值路线，直接走下一节的几何计算。

插值路线的精度上限就在那儿：方位向几百个像元一个采样点，遇到地形剧烈变化时，插出来的本地入射角基本没法用。所以这条路线只适合椭球入射角，本地入射角必须老老实实逐像元算。

3. 拿轨道状态矢量自己算，可控性最高

3.1 为什么元数据的角度不够用

除了稀疏，元数据角度还有一个隐性短板：它绑定在产品的固有网格上。如果你的工作流是把 S1 影像重采样到了别的投影、别的分辨率，或者只裁了一个子区、做了一次配准，元数据里的角度网格就跟你的影像对不上了，必须跟着做同样的几何变换，误差会一层层累积。

自己算的好处是所有中间量都在手里：卫星位置、速度、地面点坐标、视线矢量，想输出什么就输出什么，想算哪个像元就算哪个像元，还能顺带算出斜距、局部入射角、甚至视线矢量的三个分量，做形变分解时直接拿来用。

3.2 轨道插值：Lagrange 比线性靠谱得多

S1 产品注释文件里的 orbitList 通常每 10 秒左右给一组 ECEF 下的位置和速度矢量。要算任意时刻的卫星位置，必须插值。线性插值在这里是不够的——轨道是连续光滑的曲线，线性插值会在采样点之间引入位置误差，量级能到米级甚至更大，对入射角的影响虽然不大，但对精确定位、对逐像元的斜距计算是实打实的偏差。

常用做法是 8 阶或 9 阶 Lagrange 多项式插值，取目标时刻前后各 4~5 个点拟合：

def lagrange_interp(t, ts, ys):
    """对轨道状态矢量做 Lagrange 多项式插值
    t: 目标时刻（秒，相对某个基准）
    ts: 已知时刻数组（秒）
    ys: 已知状态矢量，shape = (n, 3)
    """
    n = len(ts)
    result = np.zeros(ys.shape[1], dtype=float)
    for i in range(n):
        L = 1.0
        for j in range(n):
            if i != j:
                L *= (t - ts[j]) / (ts[i] - ts[j])
        result += L * ys[i]
    return result

实际调用时不要拿全部轨道点去拟合，取离目标时刻最近的 9 个点就够，高阶多项式用全量数据会出龙格现象。另外要注意，插值前必须把时间统一成相对某个参考时刻的秒数（用 float 而不是 datetime ），否则会掉进浮点精度陷阱。

3.3 从 ECEF 转到本地 ENU，几何就豁然开朗了

有了卫星的 ECEF 位置，再有了地面点的经纬度和高程，剩下的就是坐标变换。把卫星相对地面点的矢量从 ECEF 转到以地面点为原点的东-北-天（ENU）坐标系：

def ecef_to_enu_matrix(lon, lat):
    """构造 ECEF -> ENU 的旋转矩阵，lon/lat 为弧度"""
    sl, cl = np.sin(lon), np.cos(lon)
    sp, cp = np.sin(lat), np.cos(lat)
    return np.array([
        [-sl,      cl,      0.0],
        [-sp * cl, -sp * sl, cp],
        [ cp * cl,  cp * sl, sp]
    ])

def los_geometry(sat_ecef, gnd_ecef, lon, lat):
    """输入卫星和地面点的 ECEF 坐标，输出入射角、视线方位角、斜距"""
    r = sat_ecef - gnd_ecef            # 地面指向卫星
    M = ecef_to_enu_matrix(lon, lat)
    e, n, u = M @ r                    # 东、北、天三个分量
    rng = np.linalg.norm(r)
    inc = np.degrees(np.arccos(np.clip(u / rng, -1.0, 1.0)))
    az = np.degrees(np.arctan2(e, n)) % 360.0
    return inc, az, rng

这段代码里有两个点值得多说一句。第一， u / rng 就是视线与天顶方向的夹角余弦，而这个夹角正好等于入射角（在把地表当椭球的前提下）。第二，方位角用 arctan2(e, n) 而不是 arctan2(n, e) ，顺序反了会得到一个跟真实值差 90° 的结果，而且因为角度还在 0~360 范围内，肉眼很难发现异常。我见过有人在项目里用了半年才发现方位角全错，就是因为这一行写反了。

3.4 视线方位角和航向角，靠多普勒约束闭环

如果你只想要视线方位角，还有个更简单的近似：前面说过右视情况下视线方位角 ≈ 航向角 + 90°。航向角可以用速度矢量转到 ENU 后算：

def heading_from_velocity(vel_ecef, lon, lat):
    M = ecef_to_enu_matrix(lon, lat)
    e, n, _ = M @ vel_ecef
    return np.degrees(np.arctan2(e, n)) % 360.0

这个近似在零多普勒中心频率为零时是精确的。实际上 S1 的成像多普勒中心并不严格为零，注释文件里有 dopplerCentroid 的多项式系数，可以据此算出残余偏角。这个偏差通常在千分之几度量级，做一般性的几何分析可以忽略，但如果你的应用对视线方向特别敏感（比如高精度的形变分解），还是把多普勒项补上更稳。

还有一个常被忽略的点：速度矢量必须是 ECEF 系下的，不能直接用惯性系速度。ECEF 随地球自转，惯性系速度转过来要减去自转带来的分量，这一步漏了会让航向角偏出好几度，而且偏的方向跟纬度有关，很难察觉。

4. 逐像元出栅格：三条路线怎么选

4.1 三条路线的成本与精度对照

到了要出一张跟影像完全对齐的角度栅格这一步，工程上有三条主流路线，各有各的适用边界：

路线	精度	上手成本	适合场景
注释文件插值	中低，方位向稀疏	极低，几十行代码	快速预览、椭球入射角粗估
现成工具导出	高，内置地形校正	中，需要熟悉参数	常规业务、批量生产
自建几何解算	高，完全可控	高，需要 range-Doppler 知识	特殊投影、子区、科研定制

我的建议是：第一次做这件事用工具导出，跑通之后拿结果去校验自己写的插值代码，两者对得上再考虑自建流程。上来就手搓 range-Doppler 解算，容易卡在几何收敛上出不来。

4.2 用 SNAP 导出，最省事但要防版本陷阱

SNAP 的地形校正算子（Terrain-Correction）能在做几何校正的同时输出本地入射角、投影本地入射角和椭球入射角三个波段。图形界面里对应几个勾选项，命令行用 gpt 也可以批处理：

gpt Terrain-Correction \
  -Ssource=./S1A_IW_GRDH_1SDV_xxx.zip \
  -PdemName="SRTM 1Sec HGT" \
  -PmapProjection=EPSG:4326 \
  -PpixelSpacingInMeter=10.0 \
  -PsaveLocalIncidenceAngle=true \
  -PsaveProjectedLocalIncidenceAngle=true \
  -t ./out_terrain_corrected.tif

注意： saveLocalIncidenceAngle 这类参数名在不同 SNAP 版本之间改过，有的版本叫法不同，有的版本要求先设置 -PoutputComplex=false 。如果命令报参数不识别，别去猜，直接在图形界面里勾选好对应项，然后用"File → Export → Graph XML"把整条处理链导出成 XML，再拿这个 XML 喂给 gpt。这个套路我用了很多次，比查文档快得多，而且导出的 XML 里参数全是当前版本认可的写法。

另外提醒一点：地形校正用的是 DEM，DEM 的分辨率和垂直基准直接影响本地入射角的精度。用 30 米 DEM 和用 90 米 DEM，在山区算出来的本地入射角差异可以有 2°~3°，做定量分析时要交代清楚用的是哪套。垂直基准也要留意，DEM 给的一般是正高，而几何解算用的是椭球高，两者差一个大地水准面差距，量级在几十米。这个差异对入射角的影响很小（米级的水平位移对应的角度变化微乎其微），但对绝对定位有影响，做高精度配准时需要补上。

4.3 自建 range-Doppler 流程的核心思路

如果你要处理的是非常规投影、或者只关心某个小区域想做到极致，就得自己解 range-Doppler 方程。核心是三个约束条件联立：

1. 斜距方程 ：卫星到目标的距离等于该像元的斜距。
2. 多普勒方程 ：相对速度在视线方向的投影等于该像元的多普勒中心频率。
3. 椭球方程 ：目标点落在地球椭球面上（或者 DEM 给定的高程面上）。

实际操作时不用真的解三元方程组，用一个迭代就够：先假设目标点在一个初始高程上，用前两个方程解出经纬度，再用 DEM 插值出新的高程，代回去重算，迭代三四次就能收敛到厘米级。收敛之后，卫星位置减去目标位置得到的矢量转 ENU，就是精确的视线方向和入射角。

这条路的代码量大概在三百行左右，关键难点不在公式，在数据准备：轨道要插值、DEM 要重采样到合适的分辨率、时间基准要统一。我自己踩过的一个坑是 DEM 瓦片的拼接，如果目标区跨了两块 DEM，接缝处的高程跳变会让本地入射角出现一条明显的线，看着像影像的问题，其实是 DEM 的问题。

4.4 地形校正前后，角度差异到底有多大

很多人关心一个问题：GRD 产品已经做过地距投影，几何已经"平"了，还有没有必要再算本地入射角？

答案是看用途。如果你只是做简单的辐射定标、看大范围的后向散射分布，椭球入射角就够了，GRD 自带的几何信息完全能支撑。但如果要做地形辐射校正、做山区的定量反演，本地入射角是绕不过去的，而 GRD 产品本身不提供这个信息，必须用 DEM 重新算。

差异的量级可以参考：在坡度 20° 左右的丘陵地带，本地入射角相对椭球入射角能偏 5°~8°；换算成后向散射，按经验散射模型估算，大约对应 0.5~1.5 dB 的偏差。这个量级对于区分地物类型可能不影响，但对于做时序变化检测，足以淹没真实的信号。

5. 实测中最容易翻车的几处

5.1 单位、顺序、约定，三个低级但高频的错误

角度相关的 bug 有个共同特点：结果看起来"合理"，所以不会被发现。我整理了几个踩过的：

• 弧度与度混用 。 np.arctan2 返回弧度，忘了转度就直接往栅格里写，出来的值域在 -3.14 到 3.14 之间，一看就不对。但如果只是拿去做了某个三角函数运算，就会静默地错下去。
• arctan2 的参数顺序 。前面提到过， arctan2(e, n) 和 arctan2(n, e) 差 90°。判断方法很简单：S1 是右视，视线方位角应该在航向角顺时针 90° 附近，拿一景数据算一遍看看是否符合这个关系。
• 正负号和方向约定 。方位角是从正北顺时针还是从正东逆时针，是从目标指向卫星还是从卫星指向目标。这四个组合里只有一个是对的。写代码时把约定写在注释里，几个月后回头看能救命。

5.2 IW 子条带和 burst 重叠区

IW 模式由三个子条带拼成，相邻子条带之间有重叠，同一景里不同 burst 之间也有约 1~2 秒的时间重叠。这意味着 同一个地面位置可能对应两组不同的几何参数 ，因为卫星在不同时刻、不同视角下扫过它。

处理方式取决于你的目的。如果是生成整景的角度栅格，通常取第一个有效值或者取平均，但要记录清楚用的是哪个，因为重叠区的入射角差异可能有零点几度。如果是做 burst 级的干涉处理，那每个 burst 必须用自己那套参数，绝不能跨 burst 复用。

还有一个隐蔽的坑：不同子条带的入射角范围是错开的，拼成一整景之后，如果直接看统计直方图，会看到三个峰。这不是数据错误，是正常的。看到这种图上做分析时，记得按子条带分开统计。

5.3 轨道精度与时间基准

产品自带的轨道（restituted orbit）精度大约在几米到十几米，精密轨道（precise orbit）能到厘米级。对于入射角计算，几米的轨道误差带来的角度偏差在千分之几度量级，基本可以忽略；但如果你要做的是绝对定位配准、或者需要极高精度的斜距，就一定要用精密轨道。

用外部精密轨道有个前置动作：确认时间基准一致。S1 的时间是 UTC，精密轨道文件的时间标签也是 UTC，理论上不会有问题。但在做长时序处理时，如果跨了闰秒调整的年份，个别情况下会出现 1 秒级的时间偏移，表现为整景结果的系统性偏差。我遇到过两次，排查了半天才定位到时间对齐上。所以养成一个习惯：拿到外部轨道后，先算一个已知点的位置，跟产品注释文件里的轨道点比对一下，差得太离谱就知道时间对错了。

5.4 nodata 与边界

角度栅格的边界处理也值得说一句。地形校正后影像四周会有 nodata 区，如果角度栅格没做同样的掩膜，后续做统计时 nodata 会被当成 0 度参与计算，把均值拉低一大截。更隐蔽的是，有些工具输出的 nodata 值不是常见的 -9999 而是 0，而 0 度恰好不是一个合法的入射角值（入射角不会到 0），所以从值域上能筛掉，但如果你的处理流程没做这个检查，就会带着一堆 0 继续往下跑。

我的习惯是在生成角度栅格后立刻做三个检查：值域范围、nodata 占比、均值与 incidenceAngleMidSwath 的偏差。三项都过了再进入下一步。

6. 怎么确认你算出来的角度是对的

角度这个量最大的问题是"没有参照物"。影像对不对，肉眼能看出来；分类对不对，有样本能验证；但角度算错了 5°，图上什么异常都没有。所以必须自己造几道校验关卡。

第一道：值域与分布 。S1 的 IW 模式入射角在 29°~46° 之间，如果你的结果超出了这个范围，直接回头查代码。直方图应该是单峰、偏态、随距离向单调变化的形状，如果出现双峰或者明显的台阶，多半是子条带拼接或者插值边界出了问题。

第二道：物理关系校验 。视线方位角应该等于航向角 ± 90°（右视为加），这个关系在整景内应该稳定成立，偏差在 1° 以内。如果偏差随位置系统性变化，说明坐标变换里某一步的参考系搞错了。

第三道：与其他工具交叉验证 。拿同一景数据，用现成工具导出本地入射角，跟自己的结果做逐像元差值。平原区的差值应该在 0.1° 以内，山区因为 DEM 不同可能有 1°~2° 的差异，这个量级是正常的。如果整景差出 10° 以上，先检查是不是一个用了椭球入射角、一个用了本地入射角。

第四道：跟实际地物对照 。找一块平坦的开阔区域（大面积水体最理想，水面接近水平），那里的本地入射角应该跟椭球入射角几乎一致。如果水面上出现明显的角度起伏，说明 DEM 在那块区域有噪声（水体在 DEM 里经常被处理成平地或者异常值），这时候得到的本地入射角是不可信的。

最后分享一个小技巧：如果懒得分四次检查，可以把这些校验写成一个函数，输入角度栅格和注释文件，输出一份简单的报告，包括值域、均值、方差、跟 incidenceAngleMidSwath 的偏差、以及四个象限的统计。跑一次不到十秒，但能挡住绝大多数低级错误。我在几个项目里都保留了这个习惯，尤其是在换数据源或者换版本之后，先跑一遍报告再动手，比出了问题回头找要省事得多。