基于遙感的北部灣茅尾海岸灘時空變化研究

楊鈺文 ，盧 遠 ，2*，黃 萍

（1. 南寧師范大學 地理科學與規(guī)劃學院， 廣西 南寧 530001；2. 北部灣環(huán)境演變與資源利用教育部重點實驗室，廣西 南寧 530001）

水資源是連接自然環(huán)境與人類社會的關鍵紐帶[1]。受全球氣候變化與人類活動等因素的影響，世界上許多河流的入海水沙通量急劇減少[2-5]。這不僅影響近岸海域的鹽度、懸沙濃度和營養(yǎng)物質濃度，而且還會引起河口淺灘淤漲速率改變、河口水下三角洲發(fā)育、海岸線變遷，最終造成海灣生態(tài)環(huán)境的變化[6-9]。這些情況成為了當前地球科學，尤其是國際地圈——生物圈計劃之“海岸帶陸海相互作用”（IGBP-LOICZ）研究的熱點。海岸帶是人類社會經濟高度集中分布的地區(qū)，海岸線的變化對港口的安全、海岸帶生態(tài)環(huán)境的改變、濱海土地的侵蝕都有著極其重要的影響。潮灘是河流來沙堆積的主要場所，是潛在的土地資源，對入海水沙變化的反應十分敏感[10-11]。因此，在變化環(huán)境下探討岸線以及潮灘的演變規(guī)律尤為重要。

目前，常用的海岸線變化特征分析方法有疊加分析法、多重緩沖區(qū)法、動態(tài)分割法、基于點的計算和基于剖面的計算等[12]。其中，基于剖面的計算方法可同時在多層空間尺度上進行，對岸線變化特征的刻畫更為深刻與全面，其具體的速率計算方法又分為端點速率、平均速率、線性回歸速率與加權線性回歸速率4 種。端點法的計算量最小，僅需兩個時期的歷史海岸線就可計算岸線的長期變化特征。線性回歸法根據統(tǒng)計原則求得對應所有數據的最佳趨勢線，可較好地反映岸線變化的總體趨勢。因此，將端點法和線性回歸法相結合，能夠很好地計算出岸線的時空變化特征[13]。因潮灘地理位置特殊，受潮位影響較大，大部分為測繪的“盲區(qū)”。基于遙感的潮灘范圍獲取大多是對瞬時水邊線的分析、校正來獲取的[14-16]。基于遙感瞬時水邊線獲取潮灘范圍的方法操作較復雜、效率較低。為此，本文依據灘涂與水體的光譜差異，建立了新遙感指數—灘涂指數（Tidal Flat Index， TFI），用以自動獲取出露水面潮灘范圍。該方法無須判定水邊線，操作簡便，為潮灘的長時間、大面積變化監(jiān)測提供了新思路。

在自然與人類活動雙重影響的背景下，廣西海岸環(huán)境呈現(xiàn)脆弱性[17]，岸線沖淤交互存在、平直化趨勢明顯[18]，灘涂面積遞減[19-20]，近海水質惡化[21]等等。這一系列變化直接影響著岸線、潮灘沖淤的演變規(guī)律，脅迫海岸帶地區(qū)動植物生長，最終將導致海岸帶生態(tài)環(huán)境破壞[17]。北部灣茅尾海地處廣西海岸帶北部，深入內地，受人類活動影響強烈。近年來，茅尾海沿海流域的河流水文一致性發(fā)生顯著改變，入海水沙通量呈減少趨勢[22-23]。為探討在變化的環(huán)境下茅尾海海灣對河流來水來沙減少的響應，以茅尾海岸線和河口潮灘作為研究對象，利用DSAS 軟件中的端點法和線性回歸法分析岸線的時空分布，采用新遙感指數反演多年河口潮灘面積，并輔以潮位數據揭示了潮灘的變化趨勢，最后分析總結了在變化環(huán)境下茅尾海岸線與河口潮灘的演變規(guī)律，由此為茅尾海未來發(fā)展規(guī)劃提供科學依據。

1 研究區(qū)和數據來源

1.1 研究區(qū)概況

茅尾海地處廣西壯族自治區(qū)南部，為欽州灣內灣（108°28′～108°37′E，21°46′～21°54′N），其東面、西面、北面均與陸地相鄰，如圖1 所示。據《中國海灣志》記載[24]，茅尾海海岸線總長約115.15 km，沿岸以淤泥質海岸為主，兼有基巖和沙質、沙泥質岸的共同特征，其中東側海岸（欽江區(qū)域）受侵蝕，而西側海岸（茅嶺江區(qū)域）發(fā)生沉積。茅尾海入灣徑流包括欽江和茅嶺江，兩江攜帶來的泥沙在河口區(qū)附近沉積并且不斷向海推進發(fā)展，形成大片沙質和泥質淺灘，孕育了豐富的海岸灘涂資源。茅尾海海域潮流屬不規(guī)則全日潮，因地形較平坦，潮灘在高潮時部分被淹沒，低潮時出露。由于潮灘分布廣，灣內生長有大片的紅樹（Rhizophora apiculata Bl），林間夾長有 茳 芏 （Cyperus malaccensis）， 白 鷺（Egretta garzetta），蓑衣鶴（Anthropoides virgo），以及大量種類繁多的魚蝦貝類等生物生活在其中，盛產大蠔（Ostrea gigas thunberg），青蟹（Scylla serrata），對蝦（Penaeus orientalis）和石斑魚（Epinephelus spp）。同時，江海兩水交融，導致灣內鹽度較低，有利于漁業(yè)生態(tài)旅游的發(fā)展，例如灣內擁有自治區(qū)級重要的大蠔養(yǎng)殖基地。2010 年以來，欽州市出臺了《茅尾海綜合整治規(guī)劃》，通過綜合整治進一步保護茅尾海的海洋資源和生態(tài)環(huán)境，建設宜商宜居海灣城市，進一步推進欽州港口、新城、產業(yè)的快速發(fā)展。

圖1 研究區(qū)Landsat 影像圖

1.2 數據來源

本文采用的遙感數據共計149 期，其中Landsat系列有110 期，環(huán)境星有39 期。Landsat 系列衛(wèi)星包括 TM、ETM+、OLI，時間分辨率為 16 d，空間分辨率為30 m，時間序列為1990—2018 年，數據來源于美國地質調查局（USGS）和地理空間數據云（http：//www.gscloud.cn/）；環(huán)境星作為影像數據插補，時間分辨率為4 d，空間分辨率為30 m，時間序列為2009—2018 年，數據來源于中國資源衛(wèi)星應用中心（http：//www.cresda.com/CN/）。此外，還用到茅尾海灣內龍門港歷史實測潮位數據，時間跨度為1990—2018 年，參見海洋出版社整編的《潮汐表》[25]；以及覆蓋茅尾海的歷史海圖資料（來源于中國人民解放軍海軍海道測量局），該海圖為1965—2016 年測量數據整理得出，包含灣內水深、地形等信息。

2 研究方法

2.1 海岸線提取與分析方法

選取23 期質量較高、時間跨度為1990—2018年的Landsat 影像，結合實地調查資料、潮位信息，采用人機交互的方法對茅尾海海岸線進行解譯。最后，利用數字化海岸線分析系統(tǒng)（Digital Shoreline Analysis System，DSAS）對茅尾海歷史海岸線變化趨勢進行分析。

DSAS 是ESRI ARCGIS 的一個附加組件，它可以利用多個歷史海岸線數據來計算出海岸線變化速率，其基本原理是通過歷史海岸線到基線的距離不同來描述海岸線的變化情況。DSAS 包含13 個計算海岸線變化速率的指標，本文選取線性回歸變化速率和端點變化速率來描敘茅尾海岸線的時空變化。該DSAS 模型中參數設置如下：因岸線較曲折，共設置6 條基線，均位于海上（Offshore）；兩條剖面線之間的平均間隔為100 m，剖面線長度為1 800 m，其他設置默認。依據海岸線的剖面線與研究區(qū)地理事物空間位置關系，將研究區(qū)域海岸線劃分為3 部分：西灣茅嶺江區(qū)域（剖面線序號：1 號—152 號），灣頂北部區(qū)域（剖面線序號：153 號—466 號）以及東灣欽江區(qū)域（剖面線序號：467 號—638 號），如圖2 所示。

圖2 茅尾海岸線分區(qū)示意圖

2.1.1 線性回歸變化速率（Linear Regression，LRR）LRR 是以時間為自變量，岸線到基線的距離為因變量，利用最小二乘法來擬合出線性方程，最終計算出岸線的變化速率。LRR 能夠擬合出所有年份海岸線變化情況。其方程可以表示為[26]：

式中：y 是一個因變量，表示海岸線到基線的距離；x 為自變量，表示年份；a 為擬合的常數截距；b是回歸斜率，表示每個單位x 變化所對應的y 變化，即LRR。

2.1.2 端點變化速率（End Point Rate，EPR） EPR是指兩個時期岸線之間距離的變化速率，它可以對兩條海岸線直接進行計算分析。公式如下：

式中：i 表示最遠年份海岸線與剖面線相交的點，j 表示最近年份海岸線與剖面線相交的點；EPR(i,j)表示自最遠年份海岸線到最近年份海岸線的終點變化率，即EPR；di和dj分別表示最遠年份和最近年份的海岸線到基線的距離；ΔY(j,i)表示最近年份與最遠年份之間的時間間隔。

2.2 灘涂指數（Tidal Flat Index，TFI）

制作覆蓋研究區(qū)域內邊長為1 000 m ×1 000 m的漁網，以每個漁網的中心作為取樣點，結合影像數據、對應潮位信息、歷史海圖資料對取樣點進行解譯并判定為潮灘或水體，取樣點分布如圖3 所示。將在室內判讀的取樣點進行實地考察，最終判定水體為87 個樣點，潮灘為61 個樣點，提取取樣點所對應的藍（Blue）、綠（Green）、紅（Red）、近紅外波段（NIR）的反射率，結果如圖4～圖5。潮灘取樣點中，NIR 的值總體上均大于 Red、Green、Blue 的值；在水體取樣點中，NIR 的值總體上均小于Red，Green，Blue 的值。基于此建立遙感指數 TFI（TFI1，TFI2，TFI3），公式如下：

圖3 取樣點分布圖

圖4 水體取樣點各波段反射率對比

圖5 潮灘取樣點各波段反射率對比

式中：ANIR為近紅外波段反射率值；BBlue為藍光波段反射率值；BGreen為綠光波段反射率值；BRed為紅光波段反射率值。

圖6 取樣點TFI 值分布

計算影像的TFI 值，將取樣點TFI 值導出，如圖6所示，潮灘取樣點的TFI 值整體上均大于1，水體取樣點的TFI 值整體上小于1，因此，設置TFI 閾值為1 以區(qū)分潮灘和水體。

對比 TFI1≥1，TFI2≥1，TFI3≥1 這 3 個表達式提取潮灘分布的效果。利用潮灘誤差與水體誤差作為誤差評判標準，潮灘誤差即將潮灘錯分為水體的點占總量的百分比，水體誤差即將水體錯分為潮灘的點占總量的百分比。分別計算3 個表達式出現(xiàn)誤差的點，誤差信息如表1 所示，可知表達式TFI1≥1 提取潮灘的精度最佳。

表1 潮灘提取表達式精度評價表

3 結果分析

3.1 茅尾海岸線時空變化分析

圖7 茅尾海岸線變化率

圖8 茅尾海岸線變化強度分布圖

利用EPR 和LRR 指標對1990—2018 年的茅尾海岸線進行分析研究。結果表明：茅尾海岸線總體上呈現(xiàn)出向海推進的趨勢，尤其東西兩側的欽江區(qū)域和茅嶺江區(qū)域；北部區(qū)域則保持相對穩(wěn)定的狀態(tài)，詳見圖7。以LRR 為指標，茅尾海岸線變化強度分為：高變化（LRR≤-20 m/a）、中變化（-19 m/a≤LRR≤-10 m/a）、低變化（-9 m/a≤LRR≤4 m/a）。高變化主要集中在欽江區(qū)域中部，中變化大部分分布在茅嶺江區(qū)域中下部以及欽江區(qū)域中部，詳見圖8。

3.1.1 欽江區(qū)域岸線 對欽江區(qū)域岸線的高、中變化區(qū)域進行分析，計算每條岸線到基線的距離再求平均值，即找出岸線變化率開始突變的年份，結果如組圖9 所示，得出欽江高變化區(qū)域岸線在2013年出現(xiàn)大幅度向海推進的趨勢；中變化區(qū)域岸線在2001 年和2015 年出現(xiàn)幅度向海推進。結合岸線突變年份和已有海岸線數據將欽江區(qū)域岸線變化時期分為 1992—2011 年、2012—2013 年、2014—2018年；經計算，這3 個時期EPR 的年平均變化率分別為：-0.15 m/a，-0.78 m/a，-0.98 m/a。因此，欽江區(qū)域岸線集中在2013 年起大幅度向海推進，岸線變化率顯著增大。

圖9 1992-2018 年欽江區(qū)域岸線到基線距離變化情況

3.1.2 茅嶺江區(qū)域岸線 對茅嶺江區(qū)域岸線的高、中變化區(qū)域進行分析，計算每條岸線到基線的距離再求平均值，結果如圖10 所示，得出茅嶺江高變化區(qū)域岸線在1998 年出現(xiàn)大幅度向海推進的趨勢中變化區(qū)域岸線在2001 年出現(xiàn)幅度向海推進。結合岸線突變年份和已有海岸線數據，將茅嶺江區(qū)域岸線變化時期分為 1992—1998 年，1999—2001 年 ，2002—2018年，經計算，3 個時期EPR 的年平均變化率分別為：-0.043 m/a，-0.76 m/a，-0.18 m/a。因此，茅嶺江區(qū)域岸線集中在2001 年前后向海推進速率增大，岸線變化率發(fā)生突變。

圖10 1992—2018 年茅嶺江區(qū)域岸線到基線距離變化情況

3.2 茅尾海潮灘演變趨勢分析

3.2.1 茅尾海潮灘時空變化 因茅尾海地形坡度較平緩，潮灘受潮位影響較大，本文所提取的潮灘面積僅為出露水面部分，即對研究區(qū)域進行海陸分離而獲取潮灘面積，并輔以廣西龍門港驗潮站數據進行分析。為剔除岸線沖淤影響，潮灘研究部分采用統(tǒng)一海岸線作為河口潮灘范圍。通過計算廣西龍門港驗潮站1990—2018 年的潮位資料，茅尾海平均大潮高潮位為510.1 cm，茅尾海全部有效影像中最高潮位影像為 2005 年 10 月 11 日 2：59：05（UTC）的TM 影像，對應潮位501.0 cm，由此繪制茅尾海河口潮灘研究范圍，如圖1 所示。

利用ENVI 批處理工具按TFI1公式對149 期影像進行潮灘面積提取，結果顯示如圖11 所示。潮灘面積和潮位呈明顯的反關系，高潮位對應低潮灘面積，1991 年以來茅尾海的潮灘面積一直表現(xiàn)為劇烈波動，2013—2018 年的數據組合波動形態(tài)發(fā)生明顯改變。

圖11 茅尾海歷史潮灘面積與潮位變化趨勢圖

目前針對這種時間間隔不統(tǒng)一的數據波動形態(tài)尚未有有效的定量評價方法，故選取潮位相近的兩組影像對進行定量分析，結果如表2、圖12 所示。1991—2016 年茅尾海潮灘面積明顯下降，在潮位均為197 cm 水平下，1991 和2016 年茅尾海潮灘面積分別為60.92 km2和36.28 km2，25 年間茅尾海潮灘面積減少了40.45%。結合潮位193.0 cm 和191.6 cm水平的一對影像可以看出，茅尾海潮灘面積變化集中于 2010—2015 年，5 年內潮灘面積減少了41.04%，同時也印證了圖11 中數據波動形態(tài)在2013 年后出現(xiàn)明顯變化的論點。

3.2.2 茅尾海潮灘變化驅動因素 從圖12 可知，欽江區(qū)域潮灘形態(tài)人類活動痕跡較明顯，而茅嶺江河口三角洲處潮灘沖刷形狀較自然。2010 年欽州市出臺了《茅尾海綜合整治規(guī)劃》，對欽江河口淤泥灘進行了三期工程的清淤整治，欽江河口淤泥灘的沖刷形態(tài)主要來源于人為清淤工程[27]。

表2 茅尾海相近潮位影像對信息表

圖12 茅尾海相近潮位歷史影像對比圖

從時間和空間上剔除清淤工程的影響，選取2013 年之前的兩對潮位相近的影像進行對比分析，詳見圖13。茅尾海潮灘面積在緩慢減少，尤其在茅嶺江處的河口三角洲，受到明顯沖刷。如表3 所示，自1995—2009 年，14 年間潮灘面積減少 2.27%；2002—2007 年，5 年內潮灘面積減少3.74%。由此可知，茅嶺江河口三角洲處受河流來沙減少的影響，三角洲形態(tài)發(fā)生變化，面積在緩慢減少。

表3 茅尾海相近潮位影像對信息表

圖13 茅尾海相近潮位歷史影像對比圖

4 討論

4.1 欽江區(qū)域岸灘變化因素分析

20 世紀90 年代以來，欽江區(qū)域岸線，自2013年起大幅度向海推進；欽江區(qū)域潮灘也是在2013年起面積顯著減少，欽江區(qū)域岸線與潮灘變化具有一致性。從岸灘的突變年份來看，其變化主要是受人類活動的影響。根據2010 年欽州市出臺的《茅尾海綜合整治規(guī)劃》，茅尾海東部區(qū)域設置了濱海新城啟動區(qū)以及人工沙灘建設區(qū)，同時為改善水動力條件，加大了航道水深以及船舶的通行能力，茅尾海進行了航道及淺灘清淤工程，分三期實施[27]。因此，受填海造陸以及清淤工程影響，欽江區(qū)域岸線在2013 年起大幅度向海推進；茅尾海潮灘面積在2013 后急劇減少。

4.2 茅嶺江區(qū)域岸灘變化因素分析

西部茅嶺江河口處為紅樹林保護區(qū)，生態(tài)良好、人為活動較弱。從水動力的角度來分析，其岸灘變化主要受河向來沙減少的影響。茅尾海為半封閉式的內海海灣，北部有欽江與茅嶺江匯入，南部出海口狹窄，灣內泥沙主要為河向來沙，受海外作用影響較小。灣內動力以徑流和潮流為主，其中落潮流大于漲潮流，夏季大于冬季，灣內懸沙輸送方向與余流方向一致，無論冬夏季均有向茅嶺江區(qū)域運輸[28]，因此，茅嶺江區(qū)域岸線處一直處于淤積狀態(tài)。據資料記載，欽江和茅嶺江兩河1955—1982 年間的年平均徑流量和年平均輸沙量分別為48.60 億m3和101.80 萬t（珠江流域水文資料，1955—1982）；中國海灣志中記載欽江和茅嶺江兩河年平均徑流量和年平均輸沙量分別為27.3 億m3和58.85 萬t（中國海灣志編委會，1993）。所以，欽江和茅嶺江入海水沙通量呈減少趨勢。當河流來沙減少時，來沙量小于海洋動力攜沙量，三角洲遭受沖刷，沖刷后的淤泥因灣內徑流、潮流作用等影響堆積至茅嶺江岸線處，導致茅嶺江區(qū)域岸線變化率發(fā)生變化。

綜上所述，在變化的環(huán)境背景下人類活動對河口岸灘的影響遠遠大于其他自然因素的影響。河口潮灘對河流水沙響應的敏感性大于岸線，只有當時間和空間尺度足夠大時，河流水沙通量減少對河口岸線、潮灘的影響才能得到充分體現(xiàn)。因此，在一個入海河流系統(tǒng)中，當河流水沙通量發(fā)生變化時，河口岸線和潮灘的演變規(guī)律就會發(fā)生改變（這一改變存在一定的滯后性），海灣生態(tài)系統(tǒng)隨之遭受影響。

5 結論

本文基于Landsat 以及國產環(huán)境星，通過數字化海岸線分析系統(tǒng)和建立的新遙感指數：灘涂指數，結合歷史潮位數據的分析和討論，可得以下結論：（1）1990—2018 年，北部灣茅尾海岸線總體呈現(xiàn)出向海推進的趨勢。受填海造陸影響，欽江區(qū)域岸線在2013 年向海推進速率增大；因潮流和波浪等作用，茅嶺江區(qū)域岸線一直處于淤積狀態(tài)，2001 年前后受河流來沙減少的影響，岸線向海推進速率發(fā)生突變；（2）受清淤工程影響，茅尾海潮灘受到明顯沖刷，2013 年后面積顯著減少。從時間和空間上剔除清淤工程的影響來分析，茅尾海潮灘受河流來沙減少的影響相對較弱，受河流來沙減少影響較明顯區(qū)域在茅嶺江河口三角洲處；（3）依據潮灘和水體光譜的差異性建立新遙感指數—灘涂指數（包含TFI1，TFI2，TFI3），實現(xiàn)自動化提取出露水面潮灘范圍。該方法與現(xiàn)有方法相比具有工作量小、效率高、數據易獲取、提取精度較高等優(yōu)點，為潮灘長時間序列、大面積范圍的研究提供新的方法。該方法的不足就是只能提取出露水面部分的潮灘，水下部分潮灘的提取有待深入研究。