利用Sentinel-1A數據的納木錯湖面月際變化監測

常亞茹，張繼賢，韓文立，丁夏萌，張越

(1.山東科技大學 測繪與空間信息學院，山東 青島 266590；2.國家測繪產品質量檢驗測試中心，北京 100036；3.蘭州交通大學 測繪與地理信息學院，蘭州 730070)

0 引言

氣候變化已成為目前全球面臨的一大關鍵性問題，而湖泊作為參與水文循環和氣候變化過程的重要組成部分，是全球環境變化的敏感地帶和典型區域，它的形成與消失、擴張與收縮均受季節氣候模式和長期環境變化影響[1]。青藏高原被譽為“全球氣候變化的驅動機與放大器”“亞洲水塔”，分布著地球上海拔最高、數量最多、面積最大的高原湖泊群[2]，因此受到學者們的高度關注。2019年7月在亞洲水塔國際研討會上，中國科學家發起以青藏高原冰川和湖泊為核心的第三極水塔計劃，第三極變化情況已成為國際水資源與生態環境研究熱點之一。青藏高原地區海拔較高、人口相對稀少，湖泊的自然條件保持較好，湖面的面積變化受人類活動影響較小[3]，能夠較為真實地反映區域氣候變化情況，對監測冰川融化[4]影響和反映氣候模式變化[5]有著重要作用。納木錯湖是青藏高原第二大湖泊，是一個封閉的湖盆區，空氣稀薄，溫差較大，生態環境基本為天然狀態，人為擾動因素較小，其湖泊的變化主要受自然環境和氣候變化的影響[6]。因此，分析和研究納木錯湖面積變化趨勢與氣候、地形等因素的關系，對反饋高海拔區域生態環境和氣候變化有著深遠意義[7]。

近年來，許多學者逐漸開展對青藏高原湖泊變化的研究。相較于傳統人工測量手段成本高、數據缺乏等問題，利用遙感技術能更全面、更快速高效地獲取湖泊動態變化情況[8]，尤其在自然條件惡劣、缺乏地面觀測數據和水文觀測不連續的青藏高原高海拔地區。然而光學遙感觀測受天氣和時間的影響較大[9]，因此目前研究主要是利用有限的幾景光學遙感影像對高原湖泊提取[10]，缺乏時間上的連續性，在一定程度上難以反映高原湖泊逐年、逐季度甚至是逐月的動態變化情況。合成孔徑雷達(synthetic aperture radar，SAR)的優勢在于不受云、雨、霧的影響，具有全天候和全天時的探測能力，且能獲取不同于光學傳感器的遙感地物信息，可以實現對湖泊連續性變化監測研究[11]。目前，基于SAR影像提取水體的方法已有廣泛應用，Gong等[12]利用灰度和紋理特征的變化對洪水檢測進行研究，取得了較高的精度；龐科臣等[13]利用改進的Otsu法分割DEM，去除山體陰影對水體提取，能有效降低水體提取的虛警率；鄧瀅等[14]基于結合紋理與極化分解的面向對象方法對極化SAR圖像識別水體區域，能夠提取完整水陸邊界和提高小型水體檢測率；王慶等[15]引入紋理特征和不同極化通道之間的極化差、極化比等參數與第一主成分閾值相結合，研究鄱陽湖水體的水域面積動態變化情況；湛南渝等[16]提出將簡化的SLIC超像素分割算法運用到SAR影像水體提取，實現基于SAR數據2019年株洲洪水監測。由于高原地區受云霧和海拔高的影響較大，現有研究多為長時間序列的高原湖泊面積變化，缺少針對一年間連續12月湖泊面積變化的研究，而星載SAR能夠獲取高原多云地區的短周期數據。基于此，本文利用2018年12期Sentinel-1A雷達影像數據，以青藏高原納木錯湖為研究區域，采用面向對象分割方法，通過最大類間方差法分割閾值，并結合灰度共生矩陣提取的紋理特征，提取湖泊邊界信息，深入研究2018年納木錯湖泊逐月面積變化情況及驅動因素，為相似氣候條件的高海拔湖泊群水體提取和變化監測研究提供依據。

1 研究區域及實驗數據

1.1 研究區域

納木錯湖位于西藏自治區中部(30°30′N～30°55′N，90°16′E～91°03′E)，面積約1 920 km2，湖面平均海拔達4 718 m，最大水深超過120 m，是世界上海拔最高的大型湖泊。納木錯東部和南部是終年積雪的念青唐古拉山脈和岡底斯山脈，北部和西部為地形起伏緩和的藏北高原丘陵和湖群，湖區地勢略微低于周邊區域(圖1)。納木錯湖處于水流的聚集中心處，沿湖水系發達，有羅薩、查哈蘇太河、打爾古藏布等河流。該地區屬于亞寒帶季風半干旱氣候區，環境寒冷干燥多風，光、熱、水資源充足，年日照時數達3 000小時左右，日照強烈，雨、旱季節分明，每年6月至9月為雨季，多年平均降水量為410 mm，11月至翌年5月為旱季，對氣候波動較為敏感[17]。

圖1 納木錯湖高程分布

1.2 數據及其預處理

1)研究數據源。本文主要利用了Sentinel-1A雷達數據、DEM數據和氣象資料數據。Sentinel-1A數據來源于美國地質調查局網站(http：//glovis.usgs.gov/)，選取2018年1月至12月的12期Sentinel-1干涉寬幅模式(IW)的地距(GRD)影像為研究數據源，極化方式為VV/VH，其中GRD數據為經過多視處理、采用WGS-84橢球投影到地距的影像。圖2為該影像示例。Sentinel-1A數據及參數見表1。DEM 數據使用SRTMDEM 30 m的高程數據(數據來源于地理空間數據云)。氣象數據包括納木錯湖附近的西藏當雄、班戈、那曲和申扎四個氣象站2018年的日平均氣溫、日降水量、日平均風速、日日照時數、日平均蒸發量氣象數據，用四個氣象站點的氣象數據取平均數值表示納木錯湖地區的氣象值。

圖2 納木錯湖泊2018年Sentinel-1A影像(以1、4、7、11月為例)

表1 Sentinel-1A數據及參數表

2)數據預處理。本文利用SARscape工具對12景Sentinel-1A影像進行預處理。原始的SAR影像數據為幅度數據，需對其進行輻射定標、地形校正、裁剪、濾波和地理編碼等預處理，將實驗區域的雷達強度圖像轉換為后向散射系數圖像。

2 研究方法

2.1 最大類間方差法

表面近乎光滑的水體在雷達影像中后向散射值小，表現為暗區，灰度閾值分割法依據此特征可以解算得到圖像灰度直方圖的極值，從而確定水體的分割閾值，將圖像分為背景和目標兩部分，小于閾值部分的圖像歸為水體，大于閾值部分的歸為背景，形成二值圖。最大類間方差法是由日本學者 Otsu提出來的一種全局最優閾值確定方法，也稱為Otsu法[18]，是應用較廣泛的灰度閾值分割法之一。Otsu法的基本思想是利用圖像的灰度值，將原圖像分為背景和目標兩部分，以兩者之間的最大方差確定圖像的分割閾值。由于方差是圖像像元灰度分布均勻性的一種度量，因此，方差值越大說明圖像中目標和背景兩部分區域的差別也越大，所對應的分割效果就越好[19]。

2.2 紋理特征提取

在已有的紋理特征提取方法中，灰度共生矩陣(gray-level co-occurrence matrix，GLCM)應用最為廣泛，且其提取精度高，因此采用灰度共生矩陣提取水體紋理信息并結合水體幾何信息特征，在一定程度上能夠減少其他地物對提取水體信息的干擾。灰度共生矩陣是利用窗口內特征值作為該窗口中心像素特征值的統計方法，其通過研究灰度的空間相關特性來描述紋理特征[20]，反映了圖像中任意兩點灰度的相關性和紋理特征的統計性質。在SAR影像上，灰度共生矩陣反映圖像空間結構特征，表現區域內一致性和區域間相對性[21]，其意義是，細紋理灰度空間變化很快，粗紋理隨距離增大僅有細微變化。GLCM紋理特征值較多，其中，均值反映圖像灰度均勻性和紋理規則程度，水域、道路、建筑物區域的灰度較為均勻，在均值圖像上表現為較高亮度，適用于水體提取；均值越小表示相關地物的紋理雜亂無章、難以描述；反之值越大代表紋理規律性強、易于描述。綜上，本研究選取均值特征描述紋理信息。

2.3 結合紋理特征與閾值分割法提取湖泊

目前基于SAR數據提取水體的方法主要有灰度閾值分割法、機器學習法、濾波法和結合輔助信息的提取方法等[22]，本文結合灰度共生矩陣提取紋理特征和最大類間方差法閾值分割的面向對象分割方法提取水體。首先，對Sentinel-1A影像數據進行預處理(圖3)；其次，通過最大類間方差法確定每個影像的最佳閾值，得到水體和其他背景物二值圖；最后，將初水體提取二值圖和由灰度共生矩陣提取的紋理特征結合，補充添加的紋理信息，彌補了只采用最大類間方差法產生的相干斑噪聲和圖像各處灰度信息不均勻的缺陷。另外，青藏高原地區由于地形的原因會在SAR圖像上產生地形陰影干擾信息，結合DEM高程輔助數據，以消除由于山體陰影造成的錯分混淆現象。

圖3 研究技術路線

3 結果與分析

3.1 納木錯湖面積提取結果與評價

1)湖泊面積提取結果。本實驗中對預處理后的12景SAR影像通過最大類間方差法確定的最佳閾值分別為：-14.5、-15、-14.7、-16.5、-18.5、-18、-15.8、-16、-16.5、-15.8、-14.8、-14，根據這些閾值分別得到湖泊提取的二值圖。由圖4(a)中Otsu法提取湖泊結果發現，水體和陸地未能區分，湖泊的輪廓不清晰，存在虛警率較大的問題。值得注意的是，閾值分割方法簡單、運算速度較快，能進行水體的快速提取，但此方法缺乏空間特征考量，在區分與水體散射特征相似的地物時較容易發生誤分的情況，精度和魯棒性較差，同時由于受相干斑噪聲和圖像各處灰度信息不均勻現象的影響，難以獲取準確的閾值，尤其是高海拔地區湖泊表面有湖冰存在，只通過簡單的閾值法無法完全正確地提取湖泊的整體邊界，提取的水體邊緣線不夠光滑連續，水體和河岸邊緣線也不能明確分離，存在邊緣誤差特征。

隨著SAR影像的空間分辨率不斷提高，圖像亮度范圍隨之增大，紋理結構信息也更加豐富。本文在提取水體時引入紋理特征[23]，增強圖像的空間信息，抑制SAR圖像的斑點噪聲。將紋理灰度共生矩陣與圖像中的灰度信息相結合，使提取結果的“椒鹽現象”明顯減少。從圖4(b)湖泊提取結果發現，基于GLCM-均值方法水體提取結果邊界清晰但仍存在誤提現象，尤其在影像的右上側部分明顯，未能提取完整湖體邊界。通過最大類間方差法確定水體分割初閾值，并結合灰度共生矩陣提取紋理信息，彌補只通過閾值法在處理大面積影像時獲取的結果主觀性強、可追蹤性差的缺點，湖泊提取結果如圖4(c)所示，該方法能更好地區分出水域和其他地形，提取邊界清晰明了。

圖4 不同方法提取結果

2)精度分析與評價。

(1)不同水體提取方法精度分析。為了驗證方法的有效性和精度，采用地理國情數據人工標注的解譯水體范圍作為真實水體值，將本文方法、Otsu閾值法和基于紋理特征提取的水體范圍與真實水體提取值對比分析。引用誤差統計方法進行精度評估，準確率F較大，說明提取方法效果較好；準確率F較小代表提取方法效果較差。

如表2所示，以納木錯湖為研究區域，Otsu方法提取湖泊面積準確率最低，僅有54.73%，這是由于Otsu法未能將水體和陸地區分開，存在大量誤提現象；基于紋理特征的提取精度88.26%，基本能提取湖泊的完整邊界，但噪聲影響依舊存在；結合紋理特征與Otsu方法提取精度明顯提高，達到95.12%，能較好地獲取湖泊提取結果。

表2 不同水體提取方法精度對比

(2)不同季節的水體提取精度分析。根據SAR影像的輻射特性，水體的后向散射系數由波長、斜視角、復介電常數、粗糙度等雷達參數決定。在水體提取分析中，Sentinel-1A SAR影像主要考慮地表粗糙度對后向散射系數的影響，表面越粗糙，后射回波信號越強，后向散射系數越大；平靜的水面屬于平滑面，主要表現為鏡面反射，其后向散射強度低。圖5為納木錯不同月份的水體后向散射系數變化圖，各階段后向散射強度層次明顯，夏季風速較小且湖面處于完全消融狀態，近似為鏡面反射，其后向散射強度值明顯低于其他季節；冬季后向散射系數值最大，由于溫度降低湖面開始結冰，且冬季風速較大，導致納木錯冬季時湖面呈粗糙形態，后向散射回波信號增強。基于此，以1、4、7、10月為例，分析不同季節、不同湖面形態的水體提取精度，準確率分別為：80.15%、83.62%、89.37%、85.26%。夏季7月份時水體提取精度最高，湖面狀態在一年中最為平靜，提取邊界清晰；秋季10月份精度次之，溫度逐漸降低至零下，湖面開始結冰，表面形態發生變化，影響水體提取精度；春季4月份時精度由于受風速的影響不是很高；冬季1月份水體提取精度最低，湖面既是結冰狀態又受大風影響，導致水體邊界有所誤提。

圖5 納木錯水體后向散射系數時變特征

3.2 納木錯湖湖泊面積變化監測分析

基于本文提取水體方法得到2018年12個月份納木錯湖的空間分布范圍，分析了納木錯湖在不同月份湖面面積的變化情況，具體如下。

由圖6發現，2018年12個月內納木錯湖的面積變化呈現起伏狀態，主要經歷三個階段：從3—6月迅速增加，到6—10月平穩變化，最終10月至次年2月收縮減小。從數值上來看，9月份納木錯湖面積達到峰值，為2 020.639 km2，3月份時湖泊面積最小，只有1 980.179 km2，二者相差40.46 km2。由于納木錯地處藏北高原，每年冰封期長達5個月(完全封凍時間近3個月)，結合影像數據可直觀發現1—3月份湖泊大部分處于冰凍狀態，直到4月份開始消融，此時湖水面積開始明顯增加；夏季由于溫度上升導致冰川融水和湖冰自身消融，加之降水量增加，使得9—10月份湖泊的面積達到峰值；從11月進入冬季，溫度降低降水減少，湖泊面積發生收縮變化。

圖6 2018年納木錯湖面積統計圖

通過圖7納木錯湖12個月份提取的湖面矢量疊加圖可以觀察得到，納木錯湖泊周圍都存在邊界擴張的現象，其中有四處明顯變化區域，分別在湖區的東北部、西北部和西南部，主要是由南部的岡底斯山脈和東部的念青唐古拉山脈冰雪消融補給，以及周圍的羅薩、查哈蘇太河、打爾古藏布等河流注入引起變化。

圖7 2018年納木錯湖空間分布及面積變化

3.3 湖面面積變化驅動因素分析

1)氣象資料分析。結合2018年納木錯湖周圍的當雄、班戈、那曲和申扎四個氣象站點的日均氣溫、日均降水量、日均風速、日均蒸發量、日均日照時數等氣象數據(圖8)，分析了納木錯湖區2018年的氣候變化。從圖8氣象數據和其變化趨勢線能夠看出，2018年期間氣溫的變化起伏最為明顯，1—4月氣溫均處于0 ℃以下，5—10月溫度均處于0 ℃以上，6—9月日均溫度都可達到10 ℃，11月起氣溫再次降到0 ℃以下，且溫度與湖泊的面積的變化走向大致相似。2018年月降水的變化趨勢也較為明顯，1—3月幾乎無降水，5—9月的降水明顯增加進入雨季，其中7—8月降水最為豐富達到日均5 mm，10—12月再次進入枯水期幾乎無降水。納木錯地區風資源豐富，2018年期間大風日主要集中在12月至翌年5月，大風日占全年的一半時間以上，最大日均風速可達到7 m/s。夏季蒸發量最大，冬季蒸發量最小，蒸發量整體波動不大。湖區及附近區域日照強烈，2018年日照時數接近3 000小時，年均日照率大于65%。

圖8 2018年納木錯區域氣象數據

2)湖面面積變化與氣象資料相關性分析。通過對兩個或多個具有相關性的變量進行分析，能夠得到兩個變量之間的線性關系程度，即Pearson相關性分析。納木錯湖泊水位的變化直接反映為湖水表面積的變化，而影響納木錯湖水位變化的因素有各支流、雪山、氣溫和降雨等原因。通過將統計得到的納木錯湖湖面面積統計值與各氣象數據之間進行相關性分析，由面積-氣象數據相關性分析表3的結果分析得出：納木錯湖面積變化與氣溫(相關系數r=0.680，p<0.05)和降水(相關系數r=0.638，p<0.05)存在較強的正相關性，與平均風速存在強的負相關(相關系數r=-0.796，p<0.01)；湖面面積變化與蒸發量的相關系數為r=-0.300，p=0.344(p>0.05)，與日照時數的相關系數為r=-0.171，p=0.596(p>0.05)，因此湖面面積變化與蒸發量和日照時數變化的相關性不顯著。

表3 面積-氣象數據相關性分析表

由氣象數據和相關性分析結果可得到：納木錯湖區雨、旱季節分明，熱季和雨季均為每年的6—9月，湖水水位在9—10月達到年內穩定最大值；多風是納木錯流域及其附近地區氣候的顯著特點，風速增大加快湖泊水面蒸散發。隨著氣溫的不斷升高，可能會使湖泊補給流域的上游冰川加速消融和湖冰自身融化，溫度升高是納木錯湖泊面積增加的主要原因，其中，納木錯湖東部和南部的念青唐古拉山脈、岡底斯山脈雪冰融化為湖泊增加了水域補給；雨季降水量增加也是影響湖面積變化的重要因素，風速對湖泊面積變化起負相關作用，蒸發量和日照時數兩個變量對納木錯湖面積變化影響較小。

3)地形和相對水深分析。對納木錯湖及周圍區域地形分布進行分析發現：湖區西北側是人類居住活動區域，地勢相對湖泊較高；東南側是念青唐古拉山脈，四周多雪山，湖泊位于地勢最為低洼的區域，是水流的匯集中心，且不易流散。疊加DEM高程圖中可以觀察到，納木錯湖地處區域地勢相對較低，因此能夠推斷湖水水位增高部分會掩蓋地勢相對較低地區，水位繼續上漲，納木錯湖仍會向地勢相對較低的區域擴張。

利用2018年3月Landsat 8 OLI遙感影像數據反演得到整個納木錯湖泊的相對水深，從相對水深反演結果(圖9)發現：相對水深最小的位置是陸地，值為0為圖中白色區域，值越大表示相對水深越大，即在圖中呈現的顏色越深；對納木錯湖水深反演結果進行不同的顏色分級顯示，可以清晰明了地看出整個納木錯湖的面積變化過程，即從深紅色變到橘色再到淡黃色最后到綠色，邊界較為清晰；納木錯湖泊向水深相對較淺的區域擴張消融，主要的擴張邊界在西側和東北側方向；相對水深反演結果與2018年納木錯湖面積變化過程(圖7)基本一致。

圖9 納木錯湖相對水深反演

4)湖泊景觀形狀指數變化分析。景觀形狀指數用來描述自然景觀受外界條件包含人為因素干擾的難易程度，也能反映景觀的復雜程度和易損程度，此值越小代表湖泊景觀越容易收到外界干擾，反之亦然。其中，湖泊景觀形狀指數(lake landscape shape index，LLSI)的計算如式(1)所示。

(1)

式中：LLSI為湖泊景觀形狀指數；L為湖泊的周長；A為湖泊的面積。

由圖10得到，1、2、3、4、12月LLSI值較小，由于期間溫度降低湖面結冰和受大風影響，湖泊易受到外界干擾，LLSI值均低于2.6；5—11月LLSI值較大，是由于期間湖泊的氣候環境和湖面狀態相對穩定，受外界干擾影響降低。

圖10 2018年納木錯區域LLSI變化

4 結束語

本文基于Sentinel-1A合成孔徑雷達數據，以納木錯湖為研究區域，在Otsu閾值法對影像粗分割提取的基礎上，利用灰度共生矩陣分析影像紋理特征，并結合DEM輔助數據消除山體陰影造成的錯分混淆現象，對粗提取結果進一步優化整合，實現納木錯湖泊的準確提取，最終獲得湖泊2018年逐月面積變化情況。結合周圍四個氣象站點的氣象數據，分析湖區月際變化趨勢及原因，結論如下。

1)通過實驗對比分析得到，Otsu法提取精度為54.73%，處理速度快但精度較低；基于紋理特征方法提取精度為88.26%，可以明顯獲取湖泊邊界，但存在邊界不清晰且誤提的現象；結合紋理特征與閾值分割的方法提取精度達到95.12%，效果得到明顯提升。實驗結果證明，SAR影像結合本文方法可用于高海拔地區湖泊水體信息的提取，并能實現逐年、逐月的常態化動態變化監測。

2)2018年內，納木錯湖面積呈現起伏變化，3至6月面積增加，6至10月平穩變化，10月至次年2月收縮減小；在9至10月份時面積達到了峰值2 020.639 km2，3月份時湖泊面積最小，只有1 980.179 km2。通過湖泊面積變化與氣象數據的相關性分析，氣溫升高和降水增加等因素是納木錯湖泊面積變化的主要正相關因素，風速是影響湖面面積變化的主要負相關因素，蒸發量和日照時數變量對湖泊逐月面積變化結果影響較小。另外，對整個湖區進行相對水深反演和對湖區周圍地形分析，反演結果與面積變化過程很好地對應，在西北、西南和東北方向有明顯擴張。

綜上所述，利用合成孔徑雷達SAR數據對納木錯湖觀測，可以有效地獲取其水體短期內的連續變化信息，掌握其變化規律，分析影響湖面面積變化的因素。本文仍存在不足之處，目前只針對納木錯湖2018年逐月數據進行了水體提取，今后需擴大時間跨度，研究年度、季度甚至月度的湖泊面積變化情況，并進一步加強對相似氣候條件下湖泊的研究及對比分析。