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

汪俊濤 朱勇超,2,3,4 高 飛 李江洋 陶庭葉

1 合肥工業(yè)大學(xué)土木與水利工程學(xué)院，合肥市屯溪路193號(hào)，230009 2 東華理工大學(xué)江西省數(shù)字國土重點(diǎn)實(shí)驗(yàn)室，南昌市廣蘭大道418號(hào)，330013 3 武漢大學(xué)地球空間環(huán)境與大地測量教育部重點(diǎn)實(shí)驗(yàn)室，武漢市珞喻路129號(hào)，430079 4 武漢大學(xué)測繪遙感信息工程國家重點(diǎn)實(shí)驗(yàn)室，武漢市珞喻路129號(hào)，430079

全球?qū)Ш叫l(wèi)星系統(tǒng)反射技術(shù)(global navigation satellite system reflectometry, GNSS-R)是一種新型的遙感探測技術(shù)[1]，該技術(shù)以GNSS衛(wèi)星的L波段為信號(hào)源，可用于研究海洋測高[2]、海面風(fēng)速[3]、海冰[4]、土壤濕度[5]、植被[6]及水體[7]等地表參數(shù)。

Rajabi等[8]提出基于GNSS-R的水體分布探測方法，但由于散點(diǎn)圖不能直觀地表現(xiàn)水體區(qū)域的變化，且CYGNSS衛(wèi)星在內(nèi)陸實(shí)際測量的反射信號(hào)信噪比變化很大，難以在全球范圍內(nèi)捕捉到水體，尤其是小型水體。因此，本文在前人研究基礎(chǔ)上提出一種新的基于星載GNSS-R的水體探測方法，利用CYGNSS衛(wèi)星的信噪比數(shù)據(jù)獲取地表反射率，并以剛果盆地為例進(jìn)行研究。

1 CYGNSS與研究區(qū)域

1.1 CYGNSS

颶風(fēng)全球?qū)Ш叫l(wèi)星系統(tǒng)(CYGNSS)是美國航空航天局于2016-12發(fā)射的由8個(gè)低地球軌道航天器組成的衛(wèi)星星座，每顆衛(wèi)星上都有一個(gè)GNSS-R接收機(jī)，可以同時(shí)跟蹤和處理4個(gè)GPS信號(hào)，并利用從地表反射后獲取的GPS L1 C/A信號(hào)生成延遲多普勒?qǐng)D(DDM)。相較于單基站結(jié)構(gòu)，CYGNSS衛(wèi)星對(duì)地表粗糙度的敏感性較低[8]，其低成本、低功耗的無源傳感器可提供更短的重訪時(shí)間[9]。單個(gè)CYGNSS衛(wèi)星的信號(hào)覆蓋范圍取決于入射角及信號(hào)漫反射和鏡面散射之間的相對(duì)關(guān)系，在陸地上取決于不同的入射角，信號(hào)覆蓋范圍橫向分辨率大約為0.5～1 km，沿軌道縱向分辨率大約為5 km；對(duì)于表面非常粗糙的海洋，其空間分辨率約為25 km×25 km[10]。GPS信號(hào)對(duì)陸地上水體的相干反射導(dǎo)致CYGNSS衛(wèi)星的高地表反射率信號(hào)，相干散射體與非相干散射體之間反射信號(hào)強(qiáng)度的巨大差異增強(qiáng)了衛(wèi)星對(duì)小范圍平坦水域的分辨能力，因此可以分辨出那些較小的水體，如面積為0.1 km×0.1 km的小區(qū)域水體的散射信號(hào)依然比水體周圍粗糙表面的非相干散射信號(hào)強(qiáng)約16 dB[11]。

1.2 研究區(qū)域

剛果盆地位于非洲西部，是非洲最大的盆地[12]，由內(nèi)陸湖形成，包含剛果河流域的大部分區(qū)域，擁有世界第二大熱帶雨林，物產(chǎn)資源豐富。盆地位于10°S～10°N范圍內(nèi)，地處赤道低氣壓帶，屬于熱帶雨林氣候，通常全年高溫多雨，年降水量可達(dá)2 000 mm以上。本文研究區(qū)域位于剛果盆地西南部，范圍為17°～21°E、1°S～3°N，總面積約為199 700 km2。

2 基于CYGNSS衛(wèi)星數(shù)據(jù)的星載GNSS-R水體探測算法

2.1 CYGNSS衛(wèi)星數(shù)據(jù)選取

選擇2018-12～2019-12剛果盆地區(qū)域的CYGNSS衛(wèi)星數(shù)據(jù)產(chǎn)品2.1版本L1數(shù)據(jù)集，數(shù)據(jù)格式為NetCDF。CYGNSS衛(wèi)星每日觀測數(shù)據(jù)集都存儲(chǔ)在NetCDF文件中，本文使用由CYGNSS衛(wèi)星獲取的延遲多普勒?qǐng)D得到信噪比數(shù)據(jù)，根據(jù)發(fā)射機(jī)功率、GPS接收機(jī)和發(fā)射機(jī)在鏡面反射點(diǎn)的天線增益及接收機(jī)與發(fā)射機(jī)至鏡面反射點(diǎn)的距離等主要參數(shù)對(duì)信噪比進(jìn)行校正[13]，消除數(shù)據(jù)對(duì)不同參數(shù)的依賴性，得到地表反射率(surface reflectance，SR)的表達(dá)式為：

(1)

2.2 星載GNSS-R水體探測算法

為探索基于CYGNSS衛(wèi)星數(shù)據(jù)的星載GNSS-R技術(shù)在水體分布探測領(lǐng)域的能力，提出一種利用地表反射率生成0.01°×0.01°空間分辨率格網(wǎng)的水體掩膜圖像處理算法。該算法分為地表反射率數(shù)據(jù)網(wǎng)格化(A)、去除異常值孤立像素簇(B)、剔除無效值(C)、地表反射率圖轉(zhuǎn)換為標(biāo)準(zhǔn)偏差圖(D)、地表反射率圖轉(zhuǎn)換為方差圖(E)、隨機(jī)游走圖像分割(F)等6個(gè)步驟，最終獲得CYGNSS水體掩膜，具體流程如圖1所示。

圖1 算法流程

2.2.1 數(shù)據(jù)格網(wǎng)化

為了能從地表反射率格網(wǎng)圖中更直觀地觀察到水體的分布情況，首先對(duì)原始地表反射率值進(jìn)行處理，減去最低5%的地表反射率值的平均值；再合并2018-12～2019-12剛果盆地內(nèi)CYGNSS衛(wèi)星的地表反射率數(shù)據(jù)，并將研究區(qū)劃分格網(wǎng)。對(duì)于單顆衛(wèi)星，地表反射率數(shù)據(jù)會(huì)落入不同的格網(wǎng)單元，CYGNSS衛(wèi)星星座由8顆微型衛(wèi)星組成，因此不同衛(wèi)星的多個(gè)數(shù)據(jù)會(huì)落入同一格網(wǎng)單元，再利用區(qū)域均值算法重新計(jì)算各個(gè)格網(wǎng)單元的反射率，并將結(jié)果作為研究對(duì)象。本文研究范圍為17°～21°E、1°S～3°N，按0.01°×0.01°的空間分辨率劃分為160 000個(gè)網(wǎng)格，有約64%的格網(wǎng)單元內(nèi)有數(shù)據(jù)，沒有數(shù)據(jù)的格網(wǎng)單元有57 783個(gè)，將無數(shù)據(jù)格網(wǎng)單元設(shè)為無效值(NaN)。

2.2.2 去除孤立像素簇

由于CYGNSS衛(wèi)星的空間采樣特性，采集到的某些地表反射率值非常高，這與地球表面的特性無關(guān)，可能是GPS功率變化導(dǎo)致的[14]。為了去除這些異常值，需識(shí)別出由異常值組成的孤立像素簇，并將其移除。首先選擇一個(gè)閾值threshold用于區(qū)分水體和非水體，記為Th，并根據(jù)閾值劃分將反射率格網(wǎng)圖轉(zhuǎn)換為二值化圖像，本文利用SciPy庫中的measurements.label函數(shù)識(shí)別異常值的像素簇。調(diào)用函數(shù)的功能是實(shí)現(xiàn)連通區(qū)域標(biāo)記，對(duì)于二值化圖像來說，每個(gè)像素點(diǎn)的值一般為0/1，如果2個(gè)像素點(diǎn)位置相鄰且取值相同，那么這2個(gè)像素點(diǎn)處于同一相互連通的區(qū)域內(nèi)。一個(gè)連通區(qū)域可能包含多個(gè)像素點(diǎn)，如果將同一連通區(qū)域的像素點(diǎn)用同一數(shù)值進(jìn)行標(biāo)記，則稱為連通區(qū)域標(biāo)記。Scikit-image是Python中的圖像處理工具，使用measure模塊下的label函數(shù)實(shí)現(xiàn)連通區(qū)域標(biāo)記。

處理過程中，參數(shù)input表示需要進(jìn)行處理的二值化圖像，輸出labels為一個(gè)從0開始標(biāo)記的數(shù)組。本文將地表反射率值高于Th的像素設(shè)為1，低于Th的像素設(shè)為0，函數(shù)能識(shí)別地表反射率值格網(wǎng)圖中所有單獨(dú)的像素簇，并使用不同的整數(shù)標(biāo)記每個(gè)像素簇，簡單計(jì)算出每個(gè)像素簇中元素的數(shù)量。然后選擇一個(gè)像素簇閾值pixel cluster-size，記為Ps，將元素?cái)?shù)量低于Ps的像素簇設(shè)為無效值。表1為算法的參數(shù)及取值范圍。

表1 算法參數(shù)及取值范圍

2.2.3 剔除無效值

由于部分格網(wǎng)單元內(nèi)沒有數(shù)據(jù)或在移除異常地表反射率值像素簇時(shí)遺留下大量無法參與數(shù)值計(jì)算的無效值，為便于進(jìn)行地表反射率值格網(wǎng)圖的轉(zhuǎn)換，使用SciPy庫中最近鄰插值法為每個(gè)包含無效值的網(wǎng)格單元賦值，對(duì)地表反射率值格網(wǎng)圖中的空白部分進(jìn)行填充，以剔除無效值。最近鄰插值法是簡單的灰度值插值法，是將目標(biāo)圖像像素的灰度值等于其最近鄰輸入像素的灰度值，但插值后圖像效果并不好，放大后的圖像畫質(zhì)劣化明顯，為能更直觀地觀察到水體，使水體區(qū)域區(qū)分度更高，將地表反射率格網(wǎng)圖轉(zhuǎn)換為標(biāo)準(zhǔn)偏差圖(STD)和方差圖(VAR)。

2.2.4 標(biāo)準(zhǔn)偏差圖

標(biāo)準(zhǔn)偏差是衡量數(shù)據(jù)個(gè)體間變化大小的指標(biāo)，反映整個(gè)數(shù)據(jù)樣本相對(duì)于平均值的離散程度，具體計(jì)算公式為：

(2)

式中，μ為總體X的均值。為完成標(biāo)準(zhǔn)偏差圖的轉(zhuǎn)換，先定義一個(gè)閾值STD box-size，記為Ss，以每個(gè)格網(wǎng)單元為中心建立一個(gè)大小為Ss的矩形框，確定矩形框內(nèi)的平均地表反射率值。為獲得更好的結(jié)果，取2倍標(biāo)準(zhǔn)差作為極限誤差，對(duì)于每個(gè)格網(wǎng)單元分配值的確定，需比較矩形框的中心反射率值與格網(wǎng)內(nèi)平均地表反射率值及標(biāo)準(zhǔn)差的大小關(guān)系，因此每個(gè)網(wǎng)格單元會(huì)被分配為-2、-1、0、1或2。轉(zhuǎn)換過程會(huì)在圖像中產(chǎn)生大量新的孤立像素簇，即異常的標(biāo)準(zhǔn)偏差值像素簇，因此需要再次去除異常的標(biāo)準(zhǔn)偏差值像素簇，并填充去除異常值像素簇遺留下的無效值，完成標(biāo)準(zhǔn)偏差圖的轉(zhuǎn)換。當(dāng)去除異常值像素簇時(shí)，設(shè)置閾值為0，移除高于平均值的任何一個(gè)小的像素簇。

2.2.5 方差圖

方差可以用來度量數(shù)據(jù)個(gè)體變量與均值之間的偏離程度，具體計(jì)算公式為：

(3)

方差圖的轉(zhuǎn)換與標(biāo)準(zhǔn)偏差圖的方法類似，需先定義一個(gè)閾值VAR box-size，記為Vs，即以每個(gè)格網(wǎng)單元為中心建立大小為Vs的矩形框，確定矩形框內(nèi)的平均值，將每個(gè)網(wǎng)格單元分配為0、1、2、3或4。轉(zhuǎn)換過程中也需要去除新產(chǎn)生的孤立像素簇，去除異常方差值像素簇，并對(duì)無效值進(jìn)行填充，完成地表反射率格網(wǎng)圖向方差圖的轉(zhuǎn)換。

2.2.6 隨機(jī)游走圖像分割

因本文的標(biāo)準(zhǔn)偏差圖顯示效果比方差圖好，故選擇標(biāo)準(zhǔn)偏差圖作為衡量指標(biāo)，最后使用Python的scikit-image庫中random-walker圖像分割函數(shù)，通過識(shí)別圖像中的多個(gè)單元將標(biāo)準(zhǔn)偏差圖分割為水體和旱地。本文選擇隨機(jī)游走(random-walker)圖像分割算法，因?yàn)樵诮邮誄YGNSS衛(wèi)星信號(hào)的過程中伴有噪聲誤差，而隨機(jī)游走算法特別擅長于分割有噪聲的圖像。

隨機(jī)游走算法是一種基于圖像的分割算法，屬于一種交互式的圖像分割，其分割思想是以圖像的像素為圖的頂點(diǎn)，相鄰像素之間的四鄰域或八鄰域關(guān)系為圖的邊界，并根據(jù)像素屬性及相鄰像素之間特征的相似性定義圖中各邊的權(quán)值，以構(gòu)建網(wǎng)絡(luò)圖；然后通過指定閾值標(biāo)記圖像前景和背景，即前景和背景物體的種子像素(本文中的前景和背景物體分別為水體和旱地)；再以邊上的權(quán)重為轉(zhuǎn)移概率，未標(biāo)記的像素節(jié)點(diǎn)為初始點(diǎn)，計(jì)算每個(gè)未標(biāo)記節(jié)點(diǎn)首次到達(dá)各種子像素的概率，根據(jù)概率的大小劃分未標(biāo)記節(jié)點(diǎn)，得到最終的圖像分割結(jié)果。

本文選擇高閾值Th和低閾值Tl分別標(biāo)記水體和旱地，每個(gè)大于等于Th的像素都會(huì)被標(biāo)記為Th標(biāo)簽，每個(gè)小于等于Tl的像素都會(huì)被標(biāo)記為Tl標(biāo)簽，在此設(shè)置Th為1，Tl為0。然后允許標(biāo)記像素各向異性擴(kuò)散，每個(gè)未標(biāo)記的像素都被分配到最先到達(dá)它的標(biāo)記的標(biāo)簽，得到標(biāo)準(zhǔn)偏差圖的二值化圖像，最終得到CYGNSS水體掩膜。

3 結(jié)果分析

為獲得表1中各參數(shù)的最佳取值，選取2019年剛果盆地區(qū)域MODIS水體掩膜產(chǎn)品，在具體實(shí)驗(yàn)中調(diào)整Th、Ps、Ss和Se的取值并進(jìn)行多組實(shí)驗(yàn)，以獲得最符合地表反射率值圖的水體掩膜。將生成的CYGNSS水體掩膜與MODIS水體掩膜進(jìn)行對(duì)比，發(fā)現(xiàn)最佳參數(shù)集為Th=10、Ps=8、Ss=150、Se=140。圖2為剛果盆地中2條小河流的算法演變，可以看出，2條小河流從左往右流動(dòng)。

圖2 剛果盆地一段河流演變

從圖3看出，在剛果盆地的大片湖泊流域進(jìn)行測試時(shí)，圖像處理算法效果也非常好(圖3)。

對(duì)比圖3中的MODIS水體掩膜與CYGNSS水體掩膜發(fā)現(xiàn)，MODIS水體掩膜捕捉到的河流有更為精細(xì)的細(xì)節(jié)，其分辨率達(dá)到250 m，而CYGNSS水體掩膜則相對(duì)比較粗糙，分辨率不高，但能夠明確識(shí)別出更多的支流。從圖3(c)可以看出，有些支流在MODIS水體掩膜中出現(xiàn)缺失，如18°～19°E、0°～1°N區(qū)域缺失了一條河流，主要原因是MODIS影像產(chǎn)品是基于光學(xué)遙感，其觀測會(huì)被云層或植被阻擋，而CYGNSS衛(wèi)星以L波段為信號(hào)源，能夠穿透云層和稠密的植被，識(shí)別被覆蓋范圍的水體信息。

圖3 剛果盆地CYGNSS水體掩膜與MODIS水體掩膜

2019-10～12剛果盆地進(jìn)入雨季，基本每天都會(huì)下雨，氣候潮濕悶熱，12月大量降水形成洪水，整個(gè)流域水量明顯增加。圖4(a)為2019-07～09洪水發(fā)生前通巴湖的CYGNSS水體掩膜，圖4(b)為2019-10～12通巴湖的CYGNSS水體掩膜?？梢钥闯觯珻YGNSS衛(wèi)星水體探測算法也能監(jiān)測洪水前后水淹區(qū)域的變化。

圖4 通巴湖雨季前后對(duì)比

4 結(jié) 語

為探索基于CYGNSS數(shù)據(jù)的星載GNSS-R技術(shù)在水體探測方面的能力，本文研究了2019年剛果盆地的CYGNSS衛(wèi)星數(shù)據(jù)，提出基于相干散射校正后的地表反射率作為主要研究參數(shù)的算法，將研究區(qū)域劃分為0.01°×0.01°空間分辨率格網(wǎng)，并將數(shù)據(jù)格網(wǎng)化，經(jīng)過去除孤立異常地表反射率值像素簇、填充空白網(wǎng)絡(luò)、轉(zhuǎn)換為標(biāo)準(zhǔn)偏差圖和方差圖及圖像分割等處理后，提取研究區(qū)域的CYGNSS水體掩膜。通過與MODIS水體掩膜進(jìn)行對(duì)比發(fā)現(xiàn)，本文算法探測水體的效果非常顯著，能夠識(shí)別出MODIS水體掩膜中因植被或云層遮擋而缺失的支流，并能識(shí)別內(nèi)陸水體的位置，繪制洪水發(fā)生前后水體區(qū)域的變化情況，監(jiān)測季節(jié)性洪水等短時(shí)間尺度的水文現(xiàn)象。

盡管驗(yàn)證了星載GNSS-R技術(shù)在探測水體領(lǐng)域的能力，但仍存在一些問題。如本文提出的算法在剛果盆地的測試效果非常顯著，但在一些干旱地區(qū)，本文算法很難區(qū)分植被較少的灌溉農(nóng)田與湖泊。為解決這一問題，需要根據(jù)多期CYGNSS衛(wèi)星數(shù)據(jù)繪制長期水體地圖，根據(jù)多個(gè)季節(jié)的數(shù)據(jù)集對(duì)比分析湖泊在農(nóng)田中的位置。對(duì)于內(nèi)陸水體的監(jiān)測，在地表粗糙度較高的地區(qū)進(jìn)行水體探測時(shí)，CYGNSS衛(wèi)星數(shù)據(jù)的空間分辨率仍需提高，以達(dá)到制圖產(chǎn)品規(guī)定的標(biāo)準(zhǔn)。本文算法可應(yīng)用于年、季及月等短時(shí)間尺度的水體探測，但尚未在更小的時(shí)間分辨率上(如日變化)進(jìn)行測試，后續(xù)將進(jìn)行更多的工作來完善該算法。

致謝:感謝NASA提供的CYGNSS衛(wèi)星數(shù)據(jù)和MODIS水體掩膜產(chǎn)品。