基于最大值合成和最大類間方差法萊州灣濱海灘涂變化研究

李毅， 程麗娜， 魯瑩瑩， 張博淳， 于森， 賈明明

(1.吉林建筑大學測繪與勘查工程學院，長春 130118； 2.中國科學院東北地理與農業生態研究所濕地生態與環境重點實驗室,長春 130102； 3.吉林大學地球科學學院,長春 130061； 4.長春新區北湖英才學校，長春 130000)

0 引言

全面掌握濱海灘涂變化情況的關鍵在于灘涂斑塊信息的高精度完整提取。其困難來源于衛星過境時潮位的不確定性，一定時間段內的遙感影像并不能提取出準確完整的灘涂信息。針對這一問題，國內外學者進行了一系列的相關研究。Murray等[8]利用區域潮汐模型獲取海水高潮與低潮時的遙感影像，以確定瞬時水邊線，并進行疊加，完成了東亞濱海的潮間帶濕地地圖的繪制； Wang等[9]基于GEE(Google Earth Engine)使用Landsat時序影像，通過計算像元的年度水體頻率值，確定最高潮位線和低水線的水體頻率閾值范圍來識別和劃分潮間帶，然后根據AOIs(areas of interest)的植物頻率分布閾值進行灘涂和植被的分類； Jia等[10]利用最大光譜指數合成算法(maximum spectral index composite, MSIC)和最大類間方差算法(Otsu algorithm, OTSU)，基于GEE使用Sentinel-2數據繪制了中國10 m空間分辨率灘涂分布圖； 智超等[11]借助GEE構建時序遙感指數，基于樣本分析遙感指數特征分類水體與植被，根據潮間帶不同植物物候及動力特征進行潮間帶濕地提取； Zhao等[12]通過計算長時序影像中每個像元的分位數，分別選取第5和第95分位數像素合成低、高潮位基準面影像，然后使用自適應二分類器對基準面影像進行水體和陸地的二分類，根據潮灘高潮時淹沒、低潮時露出水面的特性，繪制了中國南方的潮灘地圖； 李振等[13]建立各類土地利用類型的疊合光譜圖，通過建立決策樹法進行地物分類，繪制了膠州灣海岸帶土地利用圖； 張康永[14]基于GEE平臺，通過分析灘涂光譜信息的周期性變化規律，選取分類特征進行影像合成，綜合隨機森林分類法和面向對象的圖像分割技術，開展了杭州灣以北濱海灘涂30 a來的遙感監測和驅動力分析研究。然而目前，尚缺乏針對萊州灣地區灘涂變化的相關研究。

基于以上分析，考慮到數據的獲取性、一致性及空間分辨率，本研究借助GEE平臺構建Landsat高質量密集時序影像堆棧，基于MSIC和OTSU算法實現1990年、2000年、2010年和2020年萊州灣地區灘涂信息的高精度提取，結合灘涂周邊土地覆被變化數據集從自然和人為因素方面探討30 a間萊州灣灘涂資源的數量變化以及時空分布特征，分析其演變驅動機制，以期為該地區科學開發和合理保護灘涂資源提供技術參考和數據基礎。

1 研究區概況及其數據源

1.1 研究區概況

萊州灣位于山東省北部，地理坐標為E118°39′～120°35′，N36°55′～38°00′，主要包括東營市、濰坊市和煙臺市部分地區，內有以黃河為主的多條河流以及龍口港等重要港口，坐落著黃河三角洲國家自然保護區，擁有豐富的灘涂資源，區內氣候為溫帶季風氣候。本研究依據海岸帶環境地質調查規范(1∶25萬)，以萊州灣地區海岸線為中心，向陸海兩側做10 km緩沖區作為研究區邊界(圖1)。研究區內受不規則半日潮影響，平均潮差為1.01 m，最大潮差為2.97 m。

圖1 研究區地理位置Fig.1 Geographical location of the study area

1.2 數據源及處理

本研究以1990年、2000年、2010年和2020年為研究年份，基于GEE云計算平臺，便捷調用覆蓋研究區域目標年份及其前后1 a(時間窗口為3 a)所有符合云量要求的Landsat TM和OLI地表反射率影像集合，構建初始密集時間序列影像堆棧。利用對應的去云函數進行云掩模以獲取有效觀測像素，構成研究區域目標年份的高質量密集時間序列影像堆棧。經統計，研究區域目標年份每個像素的有效觀測數量空間分布及整體分布比例如圖2所示。2000年、2010年和2020年濱海灘涂以及其周邊土地覆蓋類型(內陸地物)驗證樣本均在Google Earth高空間分辨率影像中直接獲取。其中，采集的灘涂驗證樣本均選取低潮時期位于灘涂邊緣的點。1990年驗證樣本通過文獻檢索，查找歷史地圖等方式獲取。最終，1990年驗證樣本共1 022個(其中灘涂232個，其他土地覆被790個)； 2000年驗證樣本共1 136個(其中灘涂251個，其他土地覆被885個)； 2010年驗證樣本共1 127個(其中灘涂281個，其他土地覆被846個)； 2020年驗證樣本共1 068個(其中灘涂231個，其他土地覆被837個)。

(a) 1990年 (b) 2000年圖2-1 研究區域目標年份有效觀測數量空間分布及比例Fig.2-1 Interannual distribution and proportion of each image stack

(c) 2010年 (d) 2020年圖2-2 研究區域目標年份有效觀測數量空間分布及比例Fig.2-2 Interannual distribution and proportion of each image stack

2 研究方法

2.1 技術路線

由于潮汐的不確定性，為了提取完整的灘涂斑塊信息，本研究結合前人的相關研究[1,10]。首先，通過構建Landsat高質量密集時序影像堆棧代替單景影像獲取遙感信息的方法； 其次，基于MSIC算法獲取研究區目標年份最高、最低潮影像； 最后，利用OTSU算法進行圖像二值化分類，實現1990年、2000年、2010年和2020年萊州灣地區灘涂信息的高精度提取。結合模糊邏輯分析方法(fuzzy-based segmentation parameter，FbSP)和面向對象分析技術獲取灘涂的周邊土地覆被變化數據集，在此基礎上探討近30 a間萊州灣灘涂資源的數量變化以及空間分布特征，分析該地區灘涂演變驅動機制，具體流程如圖3所示。

圖3 技術流程Fig.3 Technical flowchart

GEE是谷歌公司推出的一個基于云計算的全球尺度遙感數據處理平臺，其API庫中提供了MSIC方法。該方法是對影像堆棧中同一位置的所有像元，逐像元計算指定光譜指數值，并利用內置的排序函數進行排序，以最大值像元作為合成影像中對應位置的像元。對所有位置的像元重復該過程即可得到指定光譜指數最大值合成影像。

OTSU算法又名大津算法，是Otsu[15]提出的一種自動選取最優閾值進行圖像分割的方法。該方法根據圖像灰度值，自動選取使類間方差最大、類內方差最小的最優分割閾值，對圖像進行二值化分割，實現圖像自動化的二分類。

面向對象的分類是通過對影像分割，基于分割得到的對象為基本元素進行分類的方法。在分割得到的單個對象中，所有的像元均可認為是同質像元[16]。該方法不僅考慮像元的光譜特征，還考慮到影像的紋理和質地等因素。因此，面向對象的分類結果可以很好地避免噪聲的影響，具有良好的整體性[17]。在本研究中，為了提高分割精度和效率，采用FbSP自動確定多尺度分割的最佳分割參數。

依據CCUS技術的研發進展，預計2035年，第一代捕集技術的成本及能耗與目前相比降低15%～20%；2045年前后，第二代捕集技術實現商業應用，成本及能耗將比第一代技術降低10%～15%；2050年，CCUS技術實現廣泛部署，建成多個CCUS產業集群。

2.2 灘涂及周邊土地覆蓋類型提取

2.2.1 灘涂提取

本文將灘涂定義為漲潮時被海水浸沒，退潮時露出海面的砂質、淤泥或軟泥質的無植被海岸淺灘。由于水體像素的改進型歸一化差異水體指數(modified normalized difference water index，mNDWI)值明顯高于非水體像素，因此選取mNDWI指數進行MSIC合成，一方面能夠突出水體與非水體像素的特征差異，同時能夠更好地獲取目標年份水體信息的最大覆蓋區域，即最高潮影像； 另外，水體像素的歸一化差分植被指數(normalized difference vegetation index，NDVI)值低于非水體像素(植被、灘涂)，使用NDVI指數獲取最低潮影像可顯著增強植被、灘涂與水體間的特征差異，還可有效避免因潮汐的不確定性對灘涂提取產生的干擾[1]。因此，本研究采用MSIC算法，分別選用mNDWI指數、NDVI指數從Landsat高質量密集時序影像堆棧中合成目標年份的最高、最低潮影像。利用OTSU算法對最高潮影像進行二值分割，實現水體和非水體的分離，通過保留最大面積斑塊的方法獲取目標年份的最大海水面，對最低潮影像進行掩模提取得到最大潮水淹沒區影像。為了去除淹沒區植被的影響，再次采用OTSU算法進行二值分類，實現灘涂信息的高精度提取。

2.2.2 灘涂周邊土地覆蓋類型提取

針對萊州灣地區潮間帶的景觀特點，從研究的目的和遙感制圖的角度出發，將研究區中除灘涂以外的地物類型分為： 內陸植被、濱海植被、養殖池/鹽田、人工表面、農田、裸地、內陸水體、海水。

本研究基于eCognition9.2平臺，通過FbSP方法自動確定目標年份單景最優影像多尺度分割的最佳分割參數，在獲取最佳分割尺度下基于面向對象技術目視解譯獲取灘涂周邊土地覆被圖。并以2020年的灘涂周邊土地覆被圖作為底圖參考制作1990年、2000年、2010年的灘涂周邊土地覆被圖。

2.3 精度驗證

本研究基于混淆矩陣，計算總體分類精度和Kappa系數，用以評價灘涂和周邊地物的分類精度。混淆矩陣是根據像元的地表真實情況與該像元的分類情況進行比較計算得出的。總體分類精度是被正確分類的像元數與總像元數的比值，是分類結果和檢驗數據一致性的評價指標。Kappa系數采用一種多元技術，綜合考慮混淆矩陣所有因素的一種精度評價指標。通常認為介于0.61～0.80則具有高度一致性； 高于0.80則幾乎完全一致[18]。

3 結果與分析

3.1 灘涂及周邊地物分類精度

利用選取的驗證點對4期分類結果分別進行精度驗證，計算得到了各研究年份的總體分類精度和Kappa系數。表1是根據選取的驗證樣本生成的研究區2020年灘涂及周邊地物遙感解譯結果的混淆矩陣，結果表明2020年灘涂及周邊土地覆被分類結果總體精度為94%，Kappa系數為0.92。1990年、2000年和2010年分類結果的總體分類精度分別為88%，91%和90%，Kappa系數分別為0.85，0.89和0.88，表明分類結果與地表真實情況具有高度的一致性，各年份的分類精度均能滿足進一步分析的要求。

表1 2020年灘涂及周邊地物土地覆被分類結果的混淆矩陣Tab.1 Confusion matrix of land cover classification results of tidal flats and surrounding land covers in 2020

3.2 萊州灣濱海灘涂時空變化分布特征

萊州灣濱海地區1990—2020年土地覆蓋變化情況如圖4所示，疊加分析可知，2020年萊州灣地區的灘涂總面積為822.38 km2，較之1990年的1 365.84 km2減少了543.46 km2，縮減幅度約為40%。近30 a期間，該地區的灘涂面積呈現持續減少的趨勢。1990—2000年灘涂面積減少了150.94 km2； 2000—2010年間的減少幅度最大，共減少304.78 km2； 2010—2020年灘涂的縮減速率明顯減緩，10 a間僅減少了87.74 km2。由于萊州灣地區涵蓋范圍較廣，包含濰坊和煙臺等多個市，各市對于灘涂濕地的保護力度和開發程度均不同，從而導致了萊州灣地區各岸段的灘涂變化狀況之間也有所差異。根據圍墾程度和岸線侵占情況，本研究將研究區分為東營市岸段、濰坊市岸段和煙臺市岸段。從整體上來看，東營市岸段灘涂面積占比最大，濰坊市岸段次之，煙臺市岸段最小； 縮減幅度最大的是濰坊市岸段，為323.56 km2，最小的是煙臺市岸段，僅82.01 km2。東營市岸段的灘涂在3個岸段中是最廣闊的，占萊州灣灘涂總面積的一半以上，該岸段在1990—2000年的減少幅度較小，為51.82 km2，在2000—2010年的減幅最大，達到了131.40 km2，在2010—2020年有小幅度的增加，增幅為45.32 km2。該岸段除黃河口地區附近灘涂海一側邊界整體上不斷向海推進，且有較大的推進幅度外，其余地區灘涂海一側邊界均不斷后移，陸側邊界持續向海推進。濰坊市岸段的灘涂整體上寬幅變化不大，呈現較規則的條狀分布，在1990—2000年、2000—2010年、2010—2020年灘涂的面積減幅先增大后減小，同樣在2000—2010年達到峰值，為144.31 km2，灘涂面積呈持續減少的趨勢，灘涂海一側邊界緩慢的向陸推進的同時，陸一側邊界不斷向海推移。煙臺市岸段的灘涂寬幅變化較大，呈現西寬北窄的形式，是灘涂面積最小的一個岸段，在1990—2000年、2000—2010年和2010—2020年灘涂面積的減少幅度逐漸升高，分別為10.65 km2，29.09 km2和42.27 km2，灘涂面積同樣呈持續減少的趨勢，整體邊界變化與濰坊市岸段相同。

(a) 1990年 (b) 2000年

(c) 2010年 (d) 2020年圖4 1990—2020年萊州灣地區潮間帶土地覆被Fig.4 Land cover of intertidal zone in Laizho Bay from 1990 to 2020

3.3 灘涂與周邊土地覆被的相互轉化

1990—2020年灘涂與周邊地物相互轉化的情況如圖5所示，灘涂面積呈凈流失狀態，灘涂持續被侵占，最主要的侵占類型為養殖池/鹽田，萊州灣沿岸除黃河口地區以外，侵占灘涂的情況均十分嚴重。其次是濱海植被，雖總體上不斷減少，但仍有相當一部分灘涂轉化為了濱海植被。灘涂與海水、內陸水體(河流)的相互轉化量雖然較大，但面積總體變化上幾乎持平，一方面是海平面上升使灘涂邊界縮退，另一方面是河流泥沙輸入使灘涂邊界向海擴張，二者對灘涂的影響呈均勢，整體上對灘涂面積變化影響不大。人工表面對灘涂的侵占雖然對比養殖池/鹽田類型少，在1990—2000年間幾乎沒有侵占情況，然而，在2000—2010年和2010—2020年2個時段中侵占面積劇烈增加，有持續增長的趨勢，空間上主要發生在濰坊市岸段，主要表現為圍海圈地。其他地物如裸地和內陸植被等，雖然與灘涂也存在相互轉化，但轉化量較小。

圖5 1990—2020年灘涂與周邊地物相互轉化Fig.5 Mutual transformation tidal flats and surrounding features from 1990 to 2020

1990—2020年間灘涂與周邊地物的轉化量如表2所示。可見，養殖池/鹽田的面積急劇擴張，在30 a間增長了多達672.27 km2，增長率到達了45%以上。其中侵占灘涂面積為414.20 km2，灘涂面積的減少量對于養殖池/鹽田面積增長的貢獻率達到了62%，是最主要的貢獻者，在1990—2000年、2000—2010年和2010—2020年間養殖池/鹽田分別侵占灘涂135.42 km2，163.16 km2和95.04 km2，在2010年侵占量達到峰值之后減少。人工表面直接侵占灘涂的面積較少，30 a間共侵占24.68 km2。萊州灣地區互花米草等濱海濕地植被擴張侵占灘涂的面積較大[19]，近30 a間共有100.30 km2的灘涂轉化為濱海植被。灘涂分別向海水和內陸水體轉化了129.08 km2和46.06 km2，同時亦有155.58 km2海水和10.17 km2內陸水轉化為灘涂，整體上對灘涂面積影響較小。

表2 1990—2020年土地利用轉移矩陣Tab.2 Land use transition matrix from 1990 to 2020 (km2)

3.4 驅動力分析

灘涂的時空變化受到自然因素和人為因素的雙重影響。自然因素包括海平面上升、入海河流泥沙輸入以及潮流變化等； 人為因素則為圍海造田和養殖池開發等一系列開發活動。

海平面上升是灘涂面積變化重要的自然原因之一。根據2020年中國海平面公報，山東沿海海平面較常年高65 mm，比2019年高11 mm，預計未來30 a山東沿海海平面還將上升55～165 mm[20]。自然狀態下，海平面上升勢必會直接導致灘涂海側邊界被海水侵蝕，在無大量泥沙輸入的情況下，引起灘涂面積減少，這也是萊州灣除黃河口地區以外灘涂海側邊界縮退，灘涂面積減少的重要誘因。此外，泥沙輸入也是灘涂生長、變化的重要原因。萊州灣地區河流資源豐富，擁有以黃河為主的多條河流，這些河流帶來了巨量的泥沙輸入，其中黃河的泥沙輸入量占總量的99%以上[21]。黃河口地區附近落潮歷時大于漲潮歷時的潮流特性[22]，也使得入海的泥沙在黃河口附近淤積成陸，灘涂海側邊界向海延伸。正是在上述因素的共同影響下，黃河口地區灘涂面積總體呈現出與萊州灣整體上相反的增長趨勢，僅在2000—2010年期間有短暫下降，后恢復增加。2020年相較于1990年灘涂面積凈增長了8.58 km2，在一定程度上減緩了整個萊州灣地區灘涂被侵占的形勢。而萊州灣其他地區缺少足夠的泥沙輸入，并且落潮歷時小于漲潮歷時的潮流特性不利于泥沙淤積成陸[23]，在海平面上升的驅動下，導致灘涂的海側邊界不斷被侵蝕，呈現出逐漸向陸移動的特征。

人類活動也是灘涂面積變化的主要因素，萊州灣灘涂整體上在海一側邊界向陸推進的同時，人類活動不斷使陸一側邊界向海推進。在這種前追后趕的形式逼迫下，灘涂面積必然持續性的減少。根據第七次全國人口普查公報，到2020年山東省人口為1億[24]，在人口增加的壓力下，向海要地成為了沿海居民的主要途徑[25]。自1991年山東省政府提出建設海上山東以來，實施耕海牧魚戰略，重視開展灘涂養殖，沿海建立大量鹽田，灘涂開發力度逐年增加[26]。受其影響，僅1990—2000年10 a間，萊州灣地區灘涂面積就減少了150.94 km2。2000—2010年更是達到了304.78 km2的峰值，到2010年山東省的鹽田面積在沿海各省市中位居榜首，沿海養殖面積也在劇烈增加[27]。在灘涂保護刻不容緩的形勢下，2015年山東省海洋局發布《中共中央國務院關于加快推進生態文明建設的意見》，明確濕地保護目標，嚴控圍填海總量，控制海洋開發強度； 2018年國務院發布《關于加強濱海濕地保護，嚴格管控圍填海的通知》，原則上不再接受省級政府申請的萊州灣地區的圍填海項目； 2019年山東省政府發布《關于加強濱海濕地保護嚴格管控圍填海的實施方案》，嚴格管理圍填海項目，控制灘涂資源開發。基于這一系列保護措施，使得2010—2020年的灘涂侵占面積較之上一個10 a減少了約71%，侵占面積銳減為87.74 km2，侵占量降低至低谷值。

4 結論

本研究基于GEE平臺，采用MSIC-OTSU和面向對象方法進行萊州灣灘涂高精度提取，研究萊州灣地區灘涂資源數量和時空分布特征，構建灘涂及其周邊土地覆被數據集，分析1990—2020年萊州灣濱海灘涂時空演變與驅動機制，得出了如下結論：

1)萊州灣地區1990—2020年30 a間灘涂面積整體呈持續減少的趨勢，2020年灘涂面積總量為822.38 km2，比1990年減少了543.46 km2，縮減幅度約40%。在2000—2010年間變化最為劇烈，灘涂面積急劇減少了304.78 km2，在2010—2020年減小幅度變緩，僅減少了87.74 km2。濰坊市岸段灘涂面積縮減量最大，30 a間減少了323.56 km2，煙臺市岸段灘涂面積縮減量最小，僅為82.01 km2。

2)萊州灣地區灘涂變化最主要的因素是養殖池/鹽田的擴張。2020年養殖池/鹽田的總面積達到了1 493.18 km2，30 a間擴展了約45%，占灘涂減少量的62%。其次是由于植被在灘涂中的擴張，使灘涂轉化為濱海植被，轉化量為100.30 km2。灘涂與海水、內陸水體的相互轉化較少。此外，雖然當前人工表面對于濱海灘涂的侵占面積較小，但未來有繼續增加的趨勢。

3)萊州灣地區灘涂變化同時受到自然和人為驅動力影響。萊州灣黃河口以外其他地區由于不利于泥沙淤積成陸使灘涂面積增加，而黃河口地區因局部落潮歷時大于漲潮歷時的潮流特性使輸入的大量泥沙淤積成陸，灘涂海側邊界向海推進，一定程度上減緩了萊州灣整體上灘涂面積縮減的形勢。