基于CYGNSS數據的星載GNSS-R水體分布探測方法研究

汪俊濤 朱勇超,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水體掩膜產品。