李霄宇,程亮,沙紅良,肖立敏,孫林云
(1.河北工程大學(xué) 水利水電學(xué)院,河北 邯鄲 056021;2.南京水利科學(xué)研究院 水文水資源與水利工程科學(xué)國家重點實驗室,江蘇 南京 210029;3.長江水利委員會水文局長江下游水文水資源勘測局,江蘇 南京 210011)
沿海灘涂原為我國沿海漁民對淤泥質(zhì)潮間帶的俗稱,廣義上還包括潮上帶和潮下帶可供開發(fā)利用的部分(彭建等,2000)。灘涂不僅是重要的后備土地資源,也是重要的濕地資源(陸永軍等,2011)。21 世紀初期以來,在工業(yè)化快速推進和土地緊縮的形勢下,我國掀起了以城市化、工業(yè)化及港口建設(shè)為主要特征的新一輪大規(guī)模灘涂開發(fā)熱潮。灘涂資源的大規(guī)模開發(fā)利用緩解了沿海地區(qū)土地資源壓力,但也造成了沿海生態(tài)、資源、環(huán)境和災(zāi)害等諸多問題(陸永軍 等,2011;喬吉果等,2018)。分析掌握灘涂開發(fā)利用特征,對于海域空間資源綜合管理、灘涂資源開發(fā)利用及其對生態(tài)環(huán)境和水資源利用影響評估具有重要意義(李京梅等,2012;徐彩瑤等,2018;李飛等,2019)。
灘涂開發(fā)利用增加了土地資源,改變了岸線。通過土地利用和岸線的變化,可以分析揭示灘涂的開發(fā)利用進程與特征。遙感具有空間宏觀性、時間周期性與空間連續(xù)性等特點,被廣泛用于沿海城市土地利用變化分析和海岸線識別。楊小艷等(2019)基于Landsat 8 影像構(gòu)建了面向?qū)ο蟮耐恋乩梅诸惙椒ú?yīng)用在連云港市,結(jié)果表明該方法精度較高。王豐等(2014)基于Landsant TM/ETM+影像識別了天津市濱海新區(qū)土地利用,進而分析了經(jīng)濟、人口等驅(qū)動因素對土地利用變化的影響。王琎等(2016)利用Landsat TM/ETM+與SPOT 影像,監(jiān)測珠江口灣區(qū)海岸線和海岸帶土地利用的時空變化,分析兩者的關(guān)系并探究了區(qū)域海岸線變遷的主要原因。康波等(2017) 利用Landsat 遙感影像,分析了近30 年長島南五島海岸線的時空變化特征。利用遙感影像提取岸線的方法總體上可分為人工目視解譯和計算機自動解譯兩類(梁立等,2018)。目視解譯主要依據(jù)不同類型海岸的遙感解譯標志來進行(馬小峰等,2007;汪翡翠等,2019),具有精度高、海岸線連續(xù)等優(yōu)點,但耗時耗力(馬小峰等,2007);自動解譯因其效率高、可復(fù)用性強等優(yōu)點,已成為研究熱點。常見的自動解譯方法有閾值分割法、邊緣檢測法、面向?qū)ο筇崛》ā⒒诨顒虞喞P头皡^(qū)域生長法等(梁立等,2018)。其中閾值分割法利用水體指數(shù)增強遙感影像的水體信息,通過設(shè)定合適閾值區(qū)分海水與陸地,提取海水邊界得到海岸線。該方法操作較為簡便,但水體指數(shù)構(gòu)建是難點。McFeeters(1996) 利用近紅外波段和綠波段對水體信息具有較強吸收和反射性的差異,提出了歸一化差異水體指數(shù)NDWI。徐涵秋(2005)利用中紅外波段替換近紅外波段,提出了改進水體指數(shù)MNDWI,應(yīng)用結(jié)果表明,與NDWI相比MNDWI 能明顯增強水體與建筑物反差,降低土壤與建筑物噪音,更適于提取城鎮(zhèn)范圍內(nèi)的水體。Wang 等(2019)基于Landsat 遙感影像數(shù)據(jù),采用MNDWI 提取了1975—2015 年共5 個時期中國海岸線長度和海岸帶面積的時空變化信息。李宗梅等(2017)基于Landsat TM/ETM+影像數(shù)據(jù),采用目視解譯與自動解譯方式提取了福建省2002 年、2006 年和2010 年的海岸線,疊加分析海岸線變化特征和程度。孫孟昊等(2019)基于Landsat 數(shù)據(jù)利用MNDWI 提取水邊線,分析了潮汐規(guī)律對海岸水邊線信息提取的影響,結(jié)合海岸高程數(shù)據(jù)提出了潮位校正方法。魏帆等(2019)運用3S 技術(shù),采用MNDWI 提取瞬時水邊線,利用分形維數(shù)法、土地利用矩陣和回歸分析等方法,研究了近40 年圍填海活動下環(huán)渤海海岸線的動態(tài)演變特征。總結(jié)來看,現(xiàn)有研究中同時從岸線與土地利用變化兩個方面對灘涂開發(fā)利用特征進行遙感分析的較少,而且這兩個方面均存在亟待完善的問題。首先,直接由遙感影像得出的岸線,實際上是衛(wèi)星過頂時刻拍攝的海陸分界線,不是我國海域使用管理中常用的由“平均大潮高潮時水陸分界的痕跡線”所定義的海岸線(梁立等,2018;馬小峰 等,2007)。受潮位、海水懸浮泥沙、季節(jié)等影響(馬小峰等,2007),這一分界線無法科學(xué)準確地反映灘涂開發(fā)利用特征與進程。其次,灘涂圍墾、填海造陸等灘涂開發(fā)利用活動,增加了土地資源,形成了灘涂開發(fā)利用區(qū)。但現(xiàn)在通過這個區(qū)域的土地利用變化來反映灘涂開發(fā)利用特征的研究較少。
海岸線由人工岸線與自然岸線組成。其中人工岸線由擋潮防浪的海堤、碼頭港口建設(shè)、灘涂圍墾等組成(汪翡翠等,2019)。灘涂開發(fā)利用或改變已有的人工岸線,或使部分自然岸線變成了人工岸線。在遙感影像中,人工岸線形狀規(guī)則、邊線平直、多呈亮白色或灰白色,與自然岸線區(qū)別明顯,易于識別(馬小峰等,2007;汪翡翠等,2019);而且人工岸線的高程一般高于潮位,識別結(jié)果不受潮位影響,自然岸線則受潮位影響。因此,對人工岸線進行識別分析,可以消除潮位影響,能科學(xué)準確地揭示灘涂開發(fā)利用進程與特征。而且通過人工岸線還能計算圍填海面積,定位灘涂開發(fā)利用區(qū)域,為通過該區(qū)域土地利用變化進一步分析灘涂開發(fā)利用特征奠定基礎(chǔ)。為此,本文從人工岸線和灘涂開發(fā)利用區(qū)土地利用的遙感識別與變化分析兩個方面,提出灘涂開發(fā)利用進程與特征的遙感分析方法,并以天津市濱海新區(qū)為研究對象開展實例研究,對該方法進行展示和檢驗,揭示天津市濱海新區(qū)灘涂的開發(fā)利用特征。
天津市濱海新區(qū)位于渤海之濱,是北方首個自由貿(mào)易試驗區(qū)和全國綜合配套改革試驗區(qū)。行政區(qū)劃面積2 270 km2,2019 年常住人口299 萬。從全國地理信息資源目錄服務(wù)系統(tǒng)的1 頤100 萬全國基礎(chǔ)地理數(shù)據(jù)庫(https://www.webmap.cn/commres.do?method=result100W)中下載了濱海新區(qū)矢量地圖數(shù)據(jù)。該數(shù)據(jù)有多個版本,本研究使用的是最新的2019 版,整體現(xiàn)勢性為2017 年。研究表明濱海新區(qū)灘涂開發(fā)利用從2002 年左右起開始加速(喬吉果等,2018),因此,本文以2002 年作為起始年份。從地理空間數(shù)據(jù)云網(wǎng)站(http://www.gscloud.cn/) 獲取了覆蓋濱海新區(qū)的2002 年7 月10 日、2006 年5 月2 日、2010 年7 月8 日、2014年8 月20 日和2018 年5 月3 日5 個年份的Landsat TM、ETM+、OLI 遙感影像。利用ENVI 軟件對影像進行輻射定標、大氣校正等預(yù)處理,并對2010 年和2014 年遙感影像進行去條帶處理。
1.2.1 人工岸線識別
為了保證精度,并得出適宜的人工岸線識別方法,分別采用目視解譯和自動提取來識別人工岸線。目視解譯的主要依據(jù)是前述提到的人工岸線在遙感影像上的特征。自動提取通過閾值分割法,采用改進的歸一化差值水體指數(shù)(MNDWI) (徐涵秋,2005):

式中:IMNDWI為MNDWI 的數(shù)值,籽green和 籽mir分別為綠光波段與中紅外波段反射率,在Landsat TM 和ETM+影像中分別對應(yīng)于第2 和第5 波段,在OLI影像中分別對應(yīng)于第3 和第6 波段(徐涵秋,2005)。區(qū)分海水與陸地的閾值可通過反映IMNDWI的直方圖來確定。為了驗證這兩種方法的差異與精度,本文補充收集了2017 年遙感影像,并將它們識別出的岸線結(jié)果與2017 年行政區(qū)劃中岸線進行對比分析。
目視解譯時,對于未開發(fā)利用或部分開發(fā)利用的灘涂區(qū)域,將距離海水最近的海堤、圍海養(yǎng)殖岸線、道路等作為人工岸線。防波堤以防御波浪入侵為目的,不屬于灘涂開發(fā)利用,本文未將其作為人工岸線。在識別出不同年份人工岸線位置后,將前后兩年的人工岸線圍成的區(qū)域稱為圍填區(qū)域,把其面積作為圍填海面積。
1.2.2 土地利用分類識別
隨著人工岸線向海域推進,由濱海新區(qū)行政區(qū)劃陸地邊界和人工岸線圍成的空間范圍不斷擴大。為了便于分析,本文以2002 年的人工岸線作為邊界,將這一范圍劃分為陸地區(qū)和灘涂開發(fā)利用區(qū)兩個部分,分別進行土地利用識別。已有陸地區(qū)是由行政區(qū)劃陸地邊界和2002 年人工岸線圍成的區(qū)域,2002 年以后范圍保持不變;而灘涂開發(fā)利用區(qū)則是2002 年與2006—2018 年的人工岸線圍成的區(qū)域,隨灘涂開發(fā)利用而變。
土地利用分類識別前,需要確定土地利用類型分類體系。參考《土地利用現(xiàn)狀分類》 (GB/T 21010—2017) 和相關(guān)文獻(楊小艷 等,2019),并結(jié)合實際情況,將已有陸地區(qū)的土地利用類型分為耕地、建設(shè)用地、林地與草地、水域和未利用地共5 類,灘涂開發(fā)利用區(qū)分為建設(shè)用地、草地、水域和未利用地共4 類。利用ENVI 軟件支持向量機方法進行自動分類,以相應(yīng)時期的Google Earth 影像目視對照作為參考,采用總體分類精度和Kappa系數(shù),對分類結(jié)果進行精度評價。將分類結(jié)果作為原始柵格數(shù)據(jù)導(dǎo)入GIS 軟件,去除遙感影像錯誤分類的背景,繪制不同年份的土地利用圖,并統(tǒng)計不同土地利用類型的面積。
對2017 年人工岸線進行目視解譯和自動提取,并以行政區(qū)劃的岸線為標準,分析比較兩種方法的差異和精度。閾值分割法的步驟一般為圖像預(yù)處理、波段運算水體指數(shù)、去除圖像異常值、選擇合適閾值分離水體與陸地,其中確定水陸分離的閾值是關(guān)鍵。理想情況下,水體指數(shù)大于0 時為水體,但實際上需要根據(jù)不同遙感影像的具體場景對閾值進行調(diào)整,圖1 給出了2017 年遙感影像計算出的IMNDWI分布情況。通過直方圖可初步判斷在斷點-0.06 和0.4 附近像元個數(shù)發(fā)生了較大變化,經(jīng)過多次嘗試并人工剔除分離效果較差的結(jié)果,確定當IMNDWI小于0 時像元代表陸地,大于0.5 時像元代表海水,在0~0.5 之間時,則可能為灘涂。因此,分別選取了0 和0.5 作為閾值。目視解譯與兩種閾值分割出的人工岸線及行政區(qū)劃岸線對比結(jié)果見圖2。

圖1 2017 年MNDWI 頻率直方圖
由圖2 可知,在灘涂圍填、以人工岸線為主的區(qū)域,目視解譯與0.5 和0 閾值分割出的岸線基本重合。在有灘涂、以自然岸線為主的區(qū)域,閾值為0.5 的岸線更靠近海域,形狀蜿蜒曲折,是灘涂與海域的分界線,并非人工岸線,且與另外兩條岸線的位置和形狀差別很大。而閾值取0 時,由于將灘涂視為海域,識別出的岸線實際上是灘涂與陸地分界線,更靠近陸地,因而與目視解譯的人工岸線更貼近,岸線的位置與形狀也相差較小。從岸線長度看,閾值為0.5 和0 分割出的岸線及目視解譯岸線長度分別為318.40 km、332.86 km 和314.25 km。三條岸線雖然位置和形狀相差較大,但由于這些差異長度相互抵消,導(dǎo)致總長度相差不大,說明不能僅通過岸線長度來判斷方法優(yōu)劣。與行政區(qū)劃靠海側(cè)、長度為292.30 km 的岸線相比,閾值為0 與目視解譯岸線的相對誤差分別為8.9%和7.5%,誤差不大。造成誤差的主要原因是行政區(qū)劃將中心漁港、天津港東疆港區(qū)和南港工業(yè)區(qū)尚未完全閉合的堤岸當成了人工岸線,從影像上看這部分區(qū)域目前尚未圍填,而且閾值分割與目視解譯出的岸線曲折度較大,因而長度更長。但這一結(jié)果與現(xiàn)有研究比較接近,汪翡翠等(2019)由2019 年Landsat OLI影像目視解譯出的人工岸線長度為366.93 km,Wang 等(2019)由MNDWI 識別出的2015 年岸線長度為307.80 km。

圖2 2017 年岸線識別結(jié)果對比
濱海新區(qū)2017 年岸線識別結(jié)果對比表明,當灘涂未被全部圍填時,閾值分割法識別出的岸線由自然岸線與人工岸線組成,不完全是人工岸線;當采用較低閾值,把灘涂視為海域,識別出岸線是灘涂與陸地分界線,更貼近人工岸線。濱海新區(qū)由于灘涂開發(fā)利用強度大,以人工岸線為主,人工岸線和陸地與灘涂分界線基本重合,因而由較低閾值分割出的岸線與目視解譯相差不大。但當人工岸線真實位置離陸地與灘涂分界線較遠,由較低閾值分割出的岸線與人工岸線必然相差很大。此外,現(xiàn)有研究表明,閾值分割法結(jié)果易受潮位影響,特別是在濱海新區(qū)這種淤泥質(zhì)海岸。為此,對2018 年岸線也進行了識別,閾值同樣通過MNDWI 樣本點值頻率直方圖來確定,分別取-0.05 和0.55,岸線識別結(jié)果如圖3 所示。對比圖2 和圖3 發(fā)現(xiàn),2018 年影像中灘涂裸露范圍更大,0.55 閾值分割出的自然岸線與行政區(qū)劃距離更遠,明顯向海域移動,說明2018 年影像的潮位更低;在灘涂區(qū)域,-0.05 閾值岸線與目視解譯人工岸線位置、形狀的差別明顯大于2017 年,這印證了現(xiàn)有研究結(jié)論,閾值分割法識別結(jié)果受潮位影響。

圖3 2018 年岸線識別結(jié)果對比
綜合2017 年和2018 年人工岸線識別結(jié)果來看,目視解譯法由于能綜合利用遙感影像、海圖、實地調(diào)研中與岸線相關(guān)信息,結(jié)果精度更高,但也較耗時耗力;而閾值分割法得出人工岸線是初步結(jié)果,還需要利用岸線相關(guān)信息,通過目視解譯,對不合理之處進行調(diào)整,特別是灘涂尚未圍填的區(qū)域。因此,推薦采用閾值分割與目視解譯相結(jié)合方式,以保證精度,并減少人力投入。本次解譯出的2002—2018 年人工岸線如圖4 所示,人工岸線長度和圍填海面積計算結(jié)果見表1。
由圖4 和表1 可知,隨著灘涂的圍墾,人工岸線、圍填海區(qū)域不斷向渤海擴張,岸線長度和圍填面積明顯增加。人工岸線長度由2002 年的146.22 km增加到2018 年的315.69 km,共增加了169.47 km;累積圍填海面積達350.76 km2,已實現(xiàn)了《天津濱海新區(qū)城市總體規(guī)劃(2009—2020 年)》提出的340 km2填海造陸目標。陳天等(2015)按照這一規(guī)劃,繪制了濱海新區(qū)圍海造地規(guī)模圖。經(jīng)過對比,本文識別出的岸線和圍填區(qū)域與該圖基本一致,這說明本文識別結(jié)果是合理的。

圖4 2002—2018 年人工岸線對比圖

表1 2002—2018 年人工岸線長度和圍填海面積
從灘涂開發(fā)區(qū)域的空間變化(圖4) 來看,2002—2006 年圍填海集中在天津港附近,2006—2010 年雖仍然以天津港為主,但中新天津生態(tài)城臨海新城、南港工業(yè)區(qū)、中心漁港等區(qū)域也被圍填,2010—2014 年則以南港工業(yè)區(qū)為主,2014—2018 轉(zhuǎn)移到了臨港中區(qū)。總體上看,圍填區(qū)域主要集中在天津港和南港工業(yè)區(qū)。
從不同時段人工岸線增加長度和圍填海面積角度對比,濱海新區(qū)灘涂開發(fā)從2002 年開始起步,2006—2014 年圍填速度最快,2014—2018 年圍填速度明顯回落。2002—2006 年與2014—2018 年人工岸線增加長度和圍填海面積均遠小于2006—2014 年。2006—2010 年和2010—2014 年,人工岸線長度分別增加了81.02 km 和62.38 km,分別占2002—2018年的43.7%和33.6%;圍填海面積分別為167.42 km2和100.76 km2,分別占2002—2018 年的51.8%和31.1%,灘涂開發(fā)主要集中于2006—2014 年。
2002—2018 年土地利用分類識別精度評價結(jié)果如表2,平均總體精度為95.51%、Kappa 系數(shù)為0.90,滿足了精度要求。土地利用分類結(jié)果如圖5所示,陸地和灘涂開發(fā)利用區(qū)不同土地類型面積如表3 所示。

表2 遙感分類精度評價

圖5 濱海新區(qū)2002—2018 年土地利用圖

表3 濱海新區(qū)不同土地類型的面積單位:km2
由圖5 和表3 可知,隨著灘涂圍填和社會經(jīng)濟發(fā)展,陸地與灘涂開發(fā)利用區(qū)土地利用也發(fā)生了較大變化。陸地區(qū)域建設(shè)用地面積逐漸增大,耕地、水域和未利用地的面積明顯減少,土地利用類型以建設(shè)用地、水域和耕地為主。2002—2018 年,陸地區(qū)域建設(shè)用地面積共增加了294.07 km2,增幅為57.3%,耕地、水域和未利用土地的面積則分別減少了162.85 km2、105.53 km2和27.3 km2,分別降低了22.4%、14.8%和88.4%。2018 年建設(shè)用地、水域和耕地面積比例分別為40.73%、30.60%和28.41%。
灘涂開發(fā)利用區(qū)的建設(shè)用地和未利用土地面積均在持續(xù)增加,形成了以建設(shè)用地和未利用土地為主的利用結(jié)構(gòu)。建設(shè)用地面積在2006—2010 年期間由0.4 km2迅速增加至43.09 km2,隨后緩慢增長,2018 年達到了63.43 km2,占灘涂區(qū)總面積的18.1%;從空間分布上看,建設(shè)用地主要位于天津港的東疆港、南疆港和臨港區(qū)附近,說明天津港及其帶動產(chǎn)業(yè)的發(fā)展是推動建設(shè)用地面積增加的主要原因。由于圍填海面積大幅增加,且增速超過了建設(shè)用地,未利用土地面積大幅增長,其中2006—2010 年和2010—2014 年增速最快,分別增加了87.8 km2和94.5 km2,2014 年以后趨于穩(wěn)定,2018年達到228.05 km2,占灘涂區(qū)總面積的65.0%,大部分土地尚未被開發(fā)利用。部分海域雖被筑堤圍隔,但未完全回填,因此形成了水域。2006 年圍填海面積比較小,水域面積也比較小,僅為7.52 km2;此后,隨著圍填海面積增加,水域面積也在增加,2010 年達到70.5 km2,2014—2018 年維持在59 km2左右。灘涂區(qū)草地包括道路綠化帶和未利用土地上生長的荒草,面積在1.0 km2左右,由于面積小,容易受分類識別誤差影響,識別結(jié)果可能不準確。
灘涂開發(fā)利用為濱海新區(qū)提供了重要土地資源。2002—2018 年,累計圍填海350.76 km2,占總行政區(qū)劃面積(2 270 km2)的15.5%,即為濱海新區(qū)增加了15.5%的土地資源;灘涂開發(fā)利用區(qū)新增的建設(shè)用地面積為63.43 km2,占濱海新區(qū)同期新增建筑面積的20%左右;陸地區(qū)域未利用土地面積僅有3.58 km2,而灘涂開發(fā)利用區(qū)尚有228.05 km2未利用土地,若加上水域面積,則可利用的土地面積達到了286.75 km2,為濱海新區(qū)未來發(fā)展提供了寶貴的土地空間。
從人工岸線和灘涂開發(fā)利用區(qū)土地利用的遙感識別與變化分析兩個方面,提出了灘涂開發(fā)利用進程與特征的遙感分析方法。利用該方法,基于天津市濱海新區(qū)2002—2018 年Landsat 遙感影像,分析了天津市濱海新區(qū)灘涂開發(fā)利用特征,得出如下結(jié)論。
(1)該方法能夠科學(xué)全面地對灘涂開發(fā)利用進程與特征進行遙感分析。首先基于Landsat 遙感影像,采用閾值分割法初步確定人工岸線,綜合利用遙感影像、海圖、實地調(diào)研等岸線相關(guān)信息,通過目視解譯,對不合理之處進行調(diào)整,以得出符合實際的人工岸線;再由識別出的人工岸線,確定灘涂開發(fā)利用區(qū),計算圍填海面積;然后分別對陸地區(qū)和灘涂開發(fā)利用區(qū)開展土地利用分類識別,最終從人工岸線、圍填海、土地利用的時空變化,分析揭示灘涂開發(fā)利用特征。
(2)截至2018 年,濱海新區(qū)人工岸線由2002 年146.22 km 增加至315.69 km,增加了115.9%;累計圍填海350.76 km2,為濱海新區(qū)增加了15.5%的土地資源;灘涂開發(fā)利用區(qū)新增的建設(shè)用地面積為63.43 km2,超過濱海新區(qū)新增建筑面積的20%;陸地區(qū)域未利用土地面積僅3.58 km2,灘涂開發(fā)而尚未利用的土地有228.05 km2,為未來發(fā)展提供了重要后備土地資源。時間上,灘涂開發(fā)利用在2002 年起步,2006—2014 年圍填速度最快,此后明顯回落,主要集中于2006—2014 年;空間上以天津港與南港工業(yè)區(qū)附近為主。