汪俊涛 朱勇超,2,3,4 高 飞 李江洋 陶庭叶

1 合肥工业大学土木与水利工程学院，合肥市屯溪路193号，230009 2 东华理工大学江西省数字国土重点实验室，南昌市广兰大道418号，330013 3 武汉大学地球空间环境与大地测量教育部重点实验室，武汉市珞喻路129号，430079 4 武汉大学测绘遥感信息工程国家重点实验室，武汉市珞喻路129号，430079

全球导航卫星系统反射技术(global navigation satellite system reflectometry, GNSS-R)是一种新型的遥感探测技术[1]，该技术以GNSS卫星的L波段为信号源，可用于研究海洋测高[2]、海面风速[3]、海冰[4]、土壤湿度[5]、植被[6]及水体[7]等地表参数。

Rajabi等[8]提出基于GNSS-R的水体分布探测方法，但由于散点图不能直观地表现水体区域的变化，且CYGNSS卫星在内陆实际测量的反射信号信噪比变化很大，难以在全球范围内捕捉到水体，尤其是小型水体。因此，本文在前人研究基础上提出一种新的基于星载GNSS-R的水体探测方法，利用CYGNSS卫星的信噪比数据获取地表反射率，并以刚果盆地为例进行研究。

1 CYGNSS与研究区域

1.1 CYGNSS

飓风全球导航卫星系统(CYGNSS)是美国航空航天局于2016-12发射的由8个低地球轨道航天器组成的卫星星座，每颗卫星上都有一个GNSS-R接收机，可以同时跟踪和处理4个GPS信号，并利用从地表反射后获取的GPS L1 C/A信号生成延迟多普勒图(DDM)。相较于单基站结构，CYGNSS卫星对地表粗糙度的敏感性较低[8]，其低成本、低功耗的无源传感器可提供更短的重访时间[9]。单个CYGNSS卫星的信号覆盖范围取决于入射角及信号漫反射和镜面散射之间的相对关系，在陆地上取决于不同的入射角，信号覆盖范围横向分辨率大约为0.5～1 km，沿轨道纵向分辨率大约为5 km；对于表面非常粗糙的海洋，其空间分辨率约为25 km×25 km[10]。GPS信号对陆地上水体的相干反射导致CYGNSS卫星的高地表反射率信号，相干散射体与非相干散射体之间反射信号强度的巨大差异增强了卫星对小范围平坦水域的分辨能力，因此可以分辨出那些较小的水体，如面积为0.1 km×0.1 km的小区域水体的散射信号依然比水体周围粗糙表面的非相干散射信号强约16 dB[11]。

1.2 研究区域

刚果盆地位于非洲西部，是非洲最大的盆地[12]，由内陆湖形成，包含刚果河流域的大部分区域，拥有世界第二大热带雨林，物产资源丰富。盆地位于10°S～10°N范围内，地处赤道低气压带，属于热带雨林气候，通常全年高温多雨，年降水量可达2 000 mm以上。本文研究区域位于刚果盆地西南部，范围为17°～21°E、1°S～3°N，总面积约为199 700 km2。

2 基于CYGNSS卫星数据的星载GNSS-R水体探测算法

2.1 CYGNSS卫星数据选取

选择2018-12～2019-12刚果盆地区域的CYGNSS卫星数据产品2.1版本L1数据集，数据格式为NetCDF。CYGNSS卫星每日观测数据集都存储在NetCDF文件中，本文使用由CYGNSS卫星获取的延迟多普勒图得到信噪比数据，根据发射机功率、GPS接收机和发射机在镜面反射点的天线增益及接收机与发射机至镜面反射点的距离等主要参数对信噪比进行校正[13]，消除数据对不同参数的依赖性，得到地表反射率(surface reflectance，SR)的表达式为：

(1)

2.2 星载GNSS-R水体探测算法

为探索基于CYGNSS卫星数据的星载GNSS-R技术在水体分布探测领域的能力，提出一种利用地表反射率生成0.01°×0.01°空间分辨率格网的水体掩膜图像处理算法。该算法分为地表反射率数据网格化(A)、去除异常值孤立像素簇(B)、剔除无效值(C)、地表反射率图转换为标准偏差图(D)、地表反射率图转换为方差图(E)、随机游走图像分割(F)等6个步骤，最终获得CYGNSS水体掩膜，具体流程如图1所示。

图1 算法流程

2.2.1 数据格网化

为了能从地表反射率格网图中更直观地观察到水体的分布情况，首先对原始地表反射率值进行处理，减去最低5%的地表反射率值的平均值；再合并2018-12～2019-12刚果盆地内CYGNSS卫星的地表反射率数据，并将研究区划分格网。对于单颗卫星，地表反射率数据会落入不同的格网单元，CYGNSS卫星星座由8颗微型卫星组成，因此不同卫星的多个数据会落入同一格网单元，再利用区域均值算法重新计算各个格网单元的反射率，并将结果作为研究对象。本文研究范围为17°～21°E、1°S～3°N，按0.01°×0.01°的空间分辨率划分为160 000个网格，有约64%的格网单元内有数据，没有数据的格网单元有57 783个，将无数据格网单元设为无效值(NaN)。

2.2.2 去除孤立像素簇

由于CYGNSS卫星的空间采样特性，采集到的某些地表反射率值非常高，这与地球表面的特性无关，可能是GPS功率变化导致的[14]。为了去除这些异常值，需识别出由异常值组成的孤立像素簇，并将其移除。首先选择一个阈值threshold用于区分水体和非水体，记为Th，并根据阈值划分将反射率格网图转换为二值化图像，本文利用SciPy库中的measurements.label函数识别异常值的像素簇。调用函数的功能是实现连通区域标记，对于二值化图像来说，每个像素点的值一般为0/1，如果2个像素点位置相邻且取值相同，那么这2个像素点处于同一相互连通的区域内。一个连通区域可能包含多个像素点，如果将同一连通区域的像素点用同一数值进行标记，则称为连通区域标记。Scikit-image是Python中的图像处理工具，使用measure模块下的label函数实现连通区域标记。

处理过程中，参数input表示需要进行处理的二值化图像，输出labels为一个从0开始标记的数组。本文将地表反射率值高于Th的像素设为1，低于Th的像素设为0，函数能识别地表反射率值格网图中所有单独的像素簇，并使用不同的整数标记每个像素簇，简单计算出每个像素簇中元素的数量。然后选择一个像素簇阈值pixel cluster-size，记为Ps，将元素数量低于Ps的像素簇设为无效值。表1为算法的参数及取值范围。

表1 算法参数及取值范围

2.2.3 剔除无效值

由于部分格网单元内没有数据或在移除异常地表反射率值像素簇时遗留下大量无法参与数值计算的无效值，为便于进行地表反射率值格网图的转换，使用SciPy库中最近邻插值法为每个包含无效值的网格单元赋值，对地表反射率值格网图中的空白部分进行填充，以剔除无效值。最近邻插值法是简单的灰度值插值法，是将目标图像像素的灰度值等于其最近邻输入像素的灰度值，但插值后图像效果并不好，放大后的图像画质劣化明显，为能更直观地观察到水体，使水体区域区分度更高，将地表反射率格网图转换为标准偏差图(STD)和方差图(VAR)。

2.2.4 标准偏差图

标准偏差是衡量数据个体间变化大小的指标，反映整个数据样本相对于平均值的离散程度，具体计算公式为：

(2)

式中，μ为总体X的均值。为完成标准偏差图的转换，先定义一个阈值STD box-size，记为Ss，以每个格网单元为中心建立一个大小为Ss的矩形框，确定矩形框内的平均地表反射率值。为获得更好的结果，取2倍标准差作为极限误差，对于每个格网单元分配值的确定，需比较矩形框的中心反射率值与格网内平均地表反射率值及标准差的大小关系，因此每个网格单元会被分配为-2、-1、0、1或2。转换过程会在图像中产生大量新的孤立像素簇，即异常的标准偏差值像素簇，因此需要再次去除异常的标准偏差值像素簇，并填充去除异常值像素簇遗留下的无效值，完成标准偏差图的转换。当去除异常值像素簇时，设置阈值为0，移除高于平均值的任何一个小的像素簇。

2.2.5 方差图

方差可以用来度量数据个体变量与均值之间的偏离程度，具体计算公式为：

(3)

方差图的转换与标准偏差图的方法类似，需先定义一个阈值VAR box-size，记为Vs，即以每个格网单元为中心建立大小为Vs的矩形框，确定矩形框内的平均值，将每个网格单元分配为0、1、2、3或4。转换过程中也需要去除新产生的孤立像素簇，去除异常方差值像素簇，并对无效值进行填充，完成地表反射率格网图向方差图的转换。

2.2.6 随机游走图像分割

因本文的标准偏差图显示效果比方差图好，故选择标准偏差图作为衡量指标，最后使用Python的scikit-image库中random-walker图像分割函数，通过识别图像中的多个单元将标准偏差图分割为水体和旱地。本文选择随机游走(random-walker)图像分割算法，因为在接收CYGNSS卫星信号的过程中伴有噪声误差，而随机游走算法特别擅长于分割有噪声的图像。

随机游走算法是一种基于图像的分割算法，属于一种交互式的图像分割，其分割思想是以图像的像素为图的顶点，相邻像素之间的四邻域或八邻域关系为图的边界，并根据像素属性及相邻像素之间特征的相似性定义图中各边的权值，以构建网络图；然后通过指定阈值标记图像前景和背景，即前景和背景物体的种子像素(本文中的前景和背景物体分别为水体和旱地)；再以边上的权重为转移概率，未标记的像素节点为初始点，计算每个未标记节点首次到达各种子像素的概率，根据概率的大小划分未标记节点，得到最终的图像分割结果。

本文选择高阈值Th和低阈值Tl分别标记水体和旱地，每个大于等于Th的像素都会被标记为Th标签，每个小于等于Tl的像素都会被标记为Tl标签，在此设置Th为1，Tl为0。然后允许标记像素各向异性扩散，每个未标记的像素都被分配到最先到达它的标记的标签，得到标准偏差图的二值化图像，最终得到CYGNSS水体掩膜。

3 结果分析

为获得表1中各参数的最佳取值，选取2019年刚果盆地区域MODIS水体掩膜产品，在具体实验中调整Th、Ps、Ss和Se的取值并进行多组实验，以获得最符合地表反射率值图的水体掩膜。将生成的CYGNSS水体掩膜与MODIS水体掩膜进行对比，发现最佳参数集为Th=10、Ps=8、Ss=150、Se=140。图2为刚果盆地中2条小河流的算法演变，可以看出，2条小河流从左往右流动。

图2 刚果盆地一段河流演变

从图3看出，在刚果盆地的大片湖泊流域进行测试时，图像处理算法效果也非常好(图3)。

对比图3中的MODIS水体掩膜与CYGNSS水体掩膜发现，MODIS水体掩膜捕捉到的河流有更为精细的细节，其分辨率达到250 m，而CYGNSS水体掩膜则相对比较粗糙，分辨率不高，但能够明确识别出更多的支流。从图3(c)可以看出，有些支流在MODIS水体掩膜中出现缺失，如18°～19°E、0°～1°N区域缺失了一条河流，主要原因是MODIS影像产品是基于光学遥感，其观测会被云层或植被阻挡，而CYGNSS卫星以L波段为信号源，能够穿透云层和稠密的植被，识别被覆盖范围的水体信息。

图3 刚果盆地CYGNSS水体掩膜与MODIS水体掩膜

2019-10～12刚果盆地进入雨季，基本每天都会下雨，气候潮湿闷热，12月大量降水形成洪水，整个流域水量明显增加。图4(a)为2019-07～09洪水发生前通巴湖的CYGNSS水体掩膜，图4(b)为2019-10～12通巴湖的CYGNSS水体掩膜。可以看出，CYGNSS卫星水体探测算法也能监测洪水前后水淹区域的变化。

图4 通巴湖雨季前后对比

4 结 语

为探索基于CYGNSS数据的星载GNSS-R技术在水体探测方面的能力，本文研究了2019年刚果盆地的CYGNSS卫星数据，提出基于相干散射校正后的地表反射率作为主要研究参数的算法，将研究区域划分为0.01°×0.01°空间分辨率格网，并将数据格网化，经过去除孤立异常地表反射率值像素簇、填充空白网络、转换为标准偏差图和方差图及图像分割等处理后，提取研究区域的CYGNSS水体掩膜。通过与MODIS水体掩膜进行对比发现，本文算法探测水体的效果非常显著，能够识别出MODIS水体掩膜中因植被或云层遮挡而缺失的支流，并能识别内陆水体的位置，绘制洪水发生前后水体区域的变化情况，监测季节性洪水等短时间尺度的水文现象。

尽管验证了星载GNSS-R技术在探测水体领域的能力，但仍存在一些问题。如本文提出的算法在刚果盆地的测试效果非常显著，但在一些干旱地区，本文算法很难区分植被较少的灌溉农田与湖泊。为解决这一问题，需要根据多期CYGNSS卫星数据绘制长期水体地图，根据多个季节的数据集对比分析湖泊在农田中的位置。对于内陆水体的监测，在地表粗糙度较高的地区进行水体探测时，CYGNSS卫星数据的空间分辨率仍需提高，以达到制图产品规定的标准。本文算法可应用于年、季及月等短时间尺度的水体探测，但尚未在更小的时间分辨率上(如日变化)进行测试，后续将进行更多的工作来完善该算法。

致谢:感谢NASA提供的CYGNSS卫星数据和MODIS水体掩膜产品。