一种GNSS-R镜面反射点精确估计方法

俞永庆，杜 毅

(1.中国石油化工股份有限公司 胜利油田分公司，山东 东营 257000；2.北京航空航天大学 电子信息工程学院，北京 100191)

0 引言

全球导航卫星系统反射信号(Global Navigation Satellite System Reflectometry,GNSS-R)技术是一种新型遥感技术，以导航卫星信号作为信号源，通过处理目标物反射信号获得反射面物理参数[1]。与传统遥感技术相比，GNSS-R具有高时空分辨率、宽刈幅、低功耗和应用领域广泛等优势[2]，成为近些年的研究热点。

镜面反射点是反射区的中心，是确定反射区位置的参考点[3]，精确估算镜面反射点对反演目标以及码相位偏移[4]的确定具有重要意义。迄今为止，镜面反射点估计方法主要有S.C.Wu算法[5]、Gleason算法[6]和线段二分法[7]。其中线段二分法利用二分查找原理对S.C.Wu算法进行改进，获得较高的收敛速度，可快速估算镜面反射点。但是，线段二分法采用地球椭球模型，使镜面反射点收敛到基准椭球面之上。基准椭球面与实际地表之间存在较大偏差，特别是随着土壤湿度[8]、植被覆盖[9]等陆面GNSS-R应用的发展，反射点处的海拔高度已成为镜面反射点估计时不可忽视的误差源。针对上述问题，在线段二分法的基础上引入Google Elevation API提供的海拔信息作为修正量，使镜面反射点最终收敛到实际地形之上，提高镜面反射点的估计精度。

1 GNSS-R几何关系

GNSS-R几何关系如图1所示，其中O为地心，R为接收机位置，T为GNSS卫星位置，S为镜面反射点位置，M为OS延长线与线段RT的交点，αr为向量SR与向量SM之间的夹角，αt为向量ST与向量SM之间的夹角。由几何光学理论可知，S为镜面反射点的充要条件是αr=αt。

图1 GNSS-R几何关系Fig.1 GNSS-R geometric relationship

2 传统镜面反射点估计方法误差分析

大地水准面、基准椭球面和实际地面三者之间的关系如图2所示。

图2 海拔误差示意Fig.2 Schematic diagram of elevation error

P为实际地面上一点，其大地高度h是该点到基准椭球面的法线距离，其海拔高度H是该点到大地水准面的法线距离，大地高度h与海拔高度H存在如下近似关系[10]：

h≈H+Nh，

(1)

式中，Nh是大地水准面高度，即大地水准面高出基准椭球面的法线距离，由于基准椭球面是一个最接近于大地水准面的椭球面，因此P点的大地高度值可以近似视为海拔高度值，即：

h≈H。

(2)

在本文后续分析中，认为大地高度h与海拔高度H相等。由Kirchhoff几何光学模型[11-13]可知，收敛到实际地表的镜面反射点pe与收敛到基准椭球面上的镜面反射点po水平方向偏差Δd满足如下关系：

(3)

式中，a，b为常数，由卫星、接收机和镜面反射点三者所构成的几何关系决定；h为镜面反射点处的海拔高度。由式(3)可知，传统镜面反射点估计算法的估计误差与反射点处的海拔高度呈正相关。

3 修正海拔误差的镜面反射点估计算法

3.1 Google Elevation API

本文在基于线段二分法的镜面反射点估计算法基础上引入Google Elevation API提供的海拔高度进行误差修正，使镜面反射点最终收敛到实际地表。Elevation API可以提供地表任意一点的海拔数据，该接口的使用方式如图3所示，用户只需要输入目标点的纬度以及经度，便可以得到该点的海拔高度。

图3 Google Elevation API使用方式Fig.3 Application of Google Elevation API

3.2 基于线段二分法的镜面反射点估计算法

基于线段二分法的镜面反射点估计算法利用折半查找原理减少迭代次数，其具体工作流程如图4所示。首先，初始化搜索区间两端点a=R、b=T，计算ab中点M，并通过M点坐标得到其星下点S；然后，分别计算向量SR与向量SM之间的夹角αr以及向量ST与向量SM之间的夹角αt；最后，判断αr与αt的大小关系，若αr=αt则退出迭代并输出S点坐标；若αr<αt则令搜索区间端点b=M并进行下一轮计算；若αr>αt则令搜索区间端点a=M并进行下一轮计算。

图4 线段二分法工作流程Fig.4 Work flow of dichotomy of line segment

3.3 修正海拔误差的镜面反射点二分查找估计算法

访问Google Elevation API存在大约100 ms的时延，在迭代过程中频繁地调用API会严重影响计算速度。为了尽量减少API的调用次数，本文将镜面反射点估计方法分为预估计和误差修正2个阶段。其中预估计阶段与上述基于线段二分法的镜面反射点估计算法完全一致，该阶段输出收敛到基准椭球面上的镜面反射点坐标。误差修正阶段的工作流程如图5所示，首先获取预估计阶段输出镜面反射点S的纬度以及经度，访问API获取S点处的海拔高度h。然后判断h是否为零，若是则直接输出S；否则，重新进行二分查找，在确定S点坐标时将其大地高度设置为h。

图5 误差修正阶段的工作流程Fig.5 Work flow of error correction stage

4 反射事件时空分布统计

为了分析引入海拔误差修正对镜面反射点估计精度的影响，本文对传统线段二分法和引入海拔误差修正后的线段二分法进行对比分析。

4.1 固定探测目标

首先分析反射面海拔固定，改变接收机和GNSS卫星位置时的镜面反射点估算精度。反演目标设置为中国最大的内陆湖泊——青海湖，其湖面海拔约3 194 m。接收机为轨道高度695 km的LEO卫星，GNSS卫星为轨道高度21 528 km的北斗MEO卫星。选取不同的接收机和GNSS卫星位置组合，共得到5组仿真数据，具体坐标设置如表1所示。

表1 接收机及GNSS卫星位置Tab.1 Positions of receiver and GNSS satellite

分别使用未进行海拔误差修正的线段二分法和引入海拔误差修正的线段二分法估算镜面反射点，具体结果如表2所示。

表2 2种方法估计的镜面反射点坐标Tab.2 Coordinates of specular points estimated by two methods

可以看出，二者在水平方向以及竖直方向上均存在一定的偏差，其中竖直方向偏差为镜面反射点处的海拔高度。

为了进一步分析接收机和卫星位置对镜面反射点估计偏差的影响，计算引入海拔误差修正线段二分法得到的镜面反射点处的卫星高度角，并计算相同条件下2种方法输出镜面反射点之间的水平方向偏差和竖直方向偏差，结果如图6所示。由图可知，水平方向偏差与高度角呈一次函数关系，通过最小二乘拟合将数据点拟合为直线，得到水平方向偏差随高度角增加的变化率为-49.02，而竖直方向偏差与镜面反射点处的海拔高度有关，与卫星及接收机位置无关。

图6 误差与卫星高度角关系Fig.6 Relationship between error and elevation of satellite

4.2 固定接收机与GNSS卫星之间的相对位置

为了分析海拔变化对镜面反射点估算精度的影响，假设接收机运行在北纬39°线上空，轨道高度为695 km，GNSS卫星运行在北纬30°线上空，轨道高度为21 528 km，且始终保持二者位于同一经度。令接收机和GNSS卫星绕地球运行一周，以5°为步进估算镜面反射点。计算得到镜面反射点处海拔高度、2种算法输出的镜面反射点之间水平方向偏差以及竖直方向偏差，其对应关系如图7所示。由图可知，镜面反射点估算精度的提升幅度与反射点处海拔高度呈正相关。

图7 反射点处海拔高度和水平误差与经度的关系Fig.7 Relationship between elevation of specular point and horizontal error with longitude

5 结束语

在传统线段二分法的基础上引入Google Elevation API提供的海拔信息对镜面反射点进行修正，提出了一种精确的镜面反射点估计方法。通过对比实验验证了该方法可消除海拔偏差，获得更高精度的镜面反射点。该方法继承了线段二分法迭代次数少、运算时间短的优点，将镜面反射点收敛到真实地表，使得镜面反射点估算精度得到提升。但是，该方法依赖于Google Elevation API所提供数据的准确性和稳定性，并且访问该接口有约100 ms的时延，导致整体处理速度较慢，未来可通过建立离线数据库的方式改善接口调用延迟。