蘇志明,孫永福,2,3*,宋玉鵬,繆棟杰,4,宋丙輝,杜 星
(1.自然資源部 第一海洋研究所,山東 青島 266061;2.國家深海基地管理中心,山東 青島 266237;3.青島海洋科學(xué)與技術(shù)試點國家實驗室 海洋地質(zhì)過程與環(huán)境功能實驗室,山東 青島 266235;4.南京大學(xué) 地理與海洋科學(xué)學(xué)院,江蘇 南京 210023)
海岸帶是海洋與陸地系統(tǒng)交叉作用、相互影響的地理單元,是研究濱海城市發(fā)展的重要平臺[1]。河口三角洲地區(qū)是典型的河流動力和海洋動力交互控制地帶,同時又是多尺度動力結(jié)構(gòu)的耦合區(qū)[2]。三角洲岸線長度、形態(tài)和所圍陸地面積的變化會直觀地反映這種河流和海洋對河口三角洲的動力作用,海岸帶的合理利用與開發(fā)對生態(tài)環(huán)境的可持續(xù)發(fā)展有重要的影響,一直以來是研究的熱點區(qū)域[1-8]。改革開放以來,我國的大陸岸線在自然和人為等因素的影響下發(fā)生了巨大的變化。黃河三角洲是中國乃至世界各大河三角洲中海陸變遷最活躍的地區(qū),黃河尾閭不斷變遷、不斷向海延伸,海岸線蝕退淤進(jìn)交替演變,三角洲面積也在逐年發(fā)生變化。通過對黃河三角洲大陸岸線的變化檢測與成因分析,可以為大陸岸線的變遷規(guī)律研究提供借鑒和思考,對于沿海生態(tài)環(huán)境保護、資源開發(fā)和可持續(xù)發(fā)展意義重大。
遙感具有宏觀、高效、可重復(fù)觀測的優(yōu)點,結(jié)合GIS技術(shù)能夠快速準(zhǔn)確地分析大陸岸線長時間序列的變化信息,從而實現(xiàn)海岸線的動態(tài)監(jiān)測[9]。王奎峰等[10]以黃河三角洲1976—2014年的MSS、TM、ETM+影像為數(shù)據(jù)源,采用一般高潮線法提取不同時期的海岸線,結(jié)合RS遙感與GIS技術(shù),對黃河三角洲的岸線變化進(jìn)行監(jiān)測。李賀等[11]基于1976—2018年的LandSat長時間序列衛(wèi)星影像,對黃河三角洲進(jìn)行了沖淤演變的研究。王一鳴等[12]選取1973-2016年41景影像資料,借助ArcGIS軟件得到研究區(qū)域面積,并結(jié)合利津站水沙資料,利用SPSS軟件分析河口地區(qū)年造陸面積與入海年輸沙量的關(guān)系。近年來新興的遙感大數(shù)據(jù)處理云平臺Google Earth Engine,憑借內(nèi)置的大量免費數(shù)據(jù)集和超強的數(shù)據(jù)處理能力成為遙感分析有力的工具[13-16]。Kilian等[17]基于Python利用圖像的監(jiān)督分類和亞像元分類方法,結(jié)合GEE平臺編寫了CoastSat開源軟件,使用戶能夠從30年來公開的衛(wèi)星圖像中獲取全球任何沙質(zhì)海岸線的位置。在2019年Deng等[18]利用GEE平臺研究了長江流域開放水體的動態(tài)監(jiān)測,并且提出了一種新的水體提取規(guī)則MIWDR進(jìn)行長江流域的水體提取。2018年Wang等[19]利用現(xiàn)有的成熟地表水指數(shù)方法,提出一種快速估算年際地表水面積的方法,并將GEE平臺將所提出的方法應(yīng)用于1990—2017年長江中游流域地表水的變化分析中。
綜上所述,利用軟件處理遙感影像已經(jīng)能夠提取并分析黃河三角洲長時間序列的水邊線,但容易產(chǎn)生數(shù)據(jù)質(zhì)量和主觀選取造成的數(shù)據(jù)差異性,而在大數(shù)據(jù)平臺(GEE)支撐下的遙感水邊線提取技術(shù)能夠通過調(diào)用該區(qū)域所有的遙感影像進(jìn)行批量處理,在不減少影像幅數(shù)的情況下剔除掉影響水邊線提取的噪音像素,在所有的良好觀測像素中采用統(tǒng)一的水邊線提取規(guī)則進(jìn)行水體提取,可以有效減少這種差異性并且提高處理的速度。因此,本文基于遙感大數(shù)據(jù)平臺Google Earth Engine對1984—2018年35期遙感影像數(shù)據(jù)進(jìn)行處理,分析近35 a黃河三角洲面積的變化及其原因。
黃河三角洲位于黃河入??谔?北臨渤海,東靠萊州灣,面積約為5 400 km2,主要位于東營境內(nèi),地勢西南高、東北低;年均降水量530~630 mm,雨量主要集中在夏季。歷史上黃河三角洲分為古代、近代和現(xiàn)代黃河三角洲,古黃河三角洲形成于1855年以前。1855年6月黃河回流入渤海后,尾閭決口50多次,大的改道11次,分別于1855—1934年和1934年至今形成以寧海和漁洼為頂點的近代和現(xiàn)代三角洲[20]。本文中的研究區(qū)主要為現(xiàn)代黃河三角洲陸上部分,地理范圍為(118°30′~119°24′E,37°30′~38°18′N)。除了黃河河道曾經(jīng)多次變道,水沙通量也發(fā)生顯著的變化。歷史上黃河年均入海沙量約占全球入海沙量的6%,居世界第二位;含沙量達(dá)25 kg/m3,居世界第一位[21]。但是受自然因素和人類活動的共同影響,黃河入海徑流量和輸沙量自20世紀(jì)50年代以來表現(xiàn)出明顯減少的趨勢,尤其近年來急劇減少[3]。從2002年7月開始,黃河水利委員會實施了黃河調(diào)水調(diào)沙工程,利用上游的水庫調(diào)節(jié)水沙,從而減輕下游河道淤積和河床升高。黃河三角洲是黃河從上游攜帶泥沙而沉積形成,所以入??趨^(qū)域多為粉砂泥質(zhì)岸線,遠(yuǎn)離入海口區(qū)域多為人工岸線。綜合前人研究[22-24]和歷史資料,將研究區(qū)分為4個岸段:刁口河岸段、人工岸段、河口岸段和萊州灣岸段(圖1)。
圖1 研究區(qū)域 Fig.1 Study area
遙感影像的數(shù)據(jù)來源于GEE數(shù)據(jù)庫,這些數(shù)據(jù)分別是LandSat 5的TM傳感器、LandSat 7的ETM+傳感器和LandSat 8的OLI傳感器經(jīng)過大氣校正后的地表反射率數(shù)據(jù)集,覆蓋周期為16 d,空間分辨率為30 m。
徑流和泥沙數(shù)據(jù)采用1984—2018年黃河利津水文站實測的年徑流量和年輸沙量數(shù)據(jù),其中1984—2002年數(shù)據(jù)來自《東營市水利志》[25],2003—2018年數(shù)據(jù)來自于《中國河流泥沙公報》[26]。黃河三角洲降水?dāng)?shù)據(jù)來源于國家氣象科學(xué)數(shù)據(jù)中心的中國地面氣候資料日志數(shù)據(jù)集(V 3.0)(http:∥data.cma.cn),利用Matlab軟件剔除氣象數(shù)據(jù)集中的異常值,按照年份提取出山東省內(nèi)黃河流經(jīng)地區(qū)的氣象站降水?dāng)?shù)據(jù),包括章丘站、惠民站、濟南站和墾利站。
1.3.1 岸線提取
目前提取海岸線的方法主要有人機互譯法[27-29]、同月同潮位法、平均低潮位法、一般高潮線法[9]、平均高潮線法[5-6,30]和最大潮高潮線法[7]等,在實際應(yīng)用中都取得了不錯的效果。本文采用平均高潮線法進(jìn)行岸線提取,取瞬時水邊線與一般高潮線的中線作為岸線,即兩線所圍面積的一半為最終面積。岸線提取流程見圖2。
圖2 岸線提取流程Fig.2 The process of shoreline extraction
1)瞬時水邊線的提取
本文通過修正歸一化水體指數(shù)(MNDWI)提取瞬時水邊線,MNDWI是徐涵秋[31]針對McFeeters[32]提出的歸一化水體指數(shù)(NDWI)改進(jìn)的。MNDWI的數(shù)學(xué)表達(dá)式為:
式中,LGreen和LMIR分別為綠波段和中紅外波段的像元亮度值。
基于GEE利用MNDWI法提取水體的方法為:①將研究區(qū)在每一年內(nèi)所有Landsat影像進(jìn)行去云處理,只保留每景影像的有效觀測像元;②經(jīng)過裁剪和鑲嵌合成該年份研究區(qū)的無云影像;③計算MNDWI水體指數(shù);④利用閾值分割法提取水體,提取面積最大的要素轉(zhuǎn)化為矢量;⑤將數(shù)據(jù)導(dǎo)入至Arc GIS中,參考Google Earth歷史影像進(jìn)行目視解譯。
2)一般高潮線的提取
一般高潮線定義為高潮灘和中潮灘的界線,在提取一般高潮線時參考樊彥國等[5]的提取方法,區(qū)別在于處理過程的前半部分在GEE中完成,相比于ENVI等遙感圖像處理平臺,GEE的圖像處理速度顯著提升。提取過程為:①在GEE中將經(jīng)過預(yù)處理的影像進(jìn)行纓帽變換,得到亮度、綠度和濕度三個波段合成的影像;②非監(jiān)督分類,根據(jù)分類的結(jié)果在Arc Map中選取訓(xùn)練樣本,通過直方圖和散點圖評估訓(xùn)練樣本,使不同分類的顏色盡量不重疊;③監(jiān)督分類,將影像分為高潮灘、中低潮灘、水體和陸地四類;④分類后在Arc Map中目視解譯提取出一般高潮線。
1.3.2 Mann-Kendall趨勢及突變檢驗
在時間序列趨勢分析中,最初由Mann和Kendall提出了Mann-Kendall(M-K)非參數(shù)檢驗方法,此方法不需要樣本遵循一定的分布,也不受少數(shù)異常值的影響,被許多學(xué)者用于徑流、降水和氣溫等非正態(tài)分布的數(shù)據(jù)。
對一個序列通過計算得到檢驗統(tǒng)計量Z,并將其與α顯著水平下的Z1-α/2相比較,如果|Z|>Z1-α/2,就可以認(rèn)為序列中存在相應(yīng)的變化趨勢。文中取95%顯著度,Z1-α/2=1.96。在突變檢測中,統(tǒng)計量UFk為時間序列X按照順序計算的秩序列標(biāo)準(zhǔn)化參數(shù),UFk為標(biāo)準(zhǔn)正態(tài)分布,對于給定的顯著性水平α,若|UFk|>U,表明序列存在明顯的趨勢變化。UBk為按照時間序列X逆序計算的秩序列標(biāo)準(zhǔn)化參數(shù)。若統(tǒng)計量UFk大于0,則表明序列呈上升趨勢,若小于0則為下降趨勢。若UFk和UBk存在交點,且交點位于臨界點之間,那么交點對應(yīng)的時間為發(fā)生突變的時刻[33-34]。
1.3.3 交叉小波變換(Cross Wave Transform,XWT)
與傳統(tǒng)的傅里葉變換相比,小波變換能夠有效地提取信號中的時頻特性,有利于分析長時間序列的數(shù)據(jù)[35]。本文通過小波變換分析黃河三角洲面積和泥沙量徑流量以及降水量的周期性變化特征,然后通過交叉小波和小波相干進(jìn)一步分析數(shù)據(jù)之間的響應(yīng)關(guān)系。
交叉小波可以對2個時間序列在不同時頻域中的相互關(guān)系進(jìn)行分析研究[36-37],設(shè)WX(s)、WY(s)分別為給定的兩個時間序列X和Y的交叉小波變換,其交叉小波譜為,其中的復(fù)共軛,交叉小波功率譜密度為,其值越大表示兩個時間序列的相關(guān)程度越高[38]。
1.3.4 小波相干(Wavelet Transform Coherence,WTC)
小波相干譜可以反映2個時間序列在時頻空間上的局部相關(guān)程度[39],定義為:
式中:S為平滑器;s為伸縮尺度;和分別為X、Y的小波變換;為交叉小波譜。對上述小波分析的結(jié)果顯著性都采用理論紅噪聲譜進(jìn)行檢驗[39]。
為了分析黃河三角洲面積和黃河入海水沙量(包含輸沙量、徑流量和降水量)的關(guān)系,首先分析水沙量趨勢:對1984—2018年水沙時間序列數(shù)據(jù)做五年滑動平均處理,并用M-K非參數(shù)檢驗法對水沙量進(jìn)行趨勢檢驗和突變檢測。
1)采用滑動平均法,取k=5繪制出年輸沙量、年徑流量和年降水量的趨勢圖:做五年滑動平均后可看出年輸沙量(圖3)呈顯著的下降趨勢;做五年滑動平均后可看出年徑流量(圖4)在1984—2001年呈顯著下降趨勢、2001—2008年呈上升趨勢、2009—2018年總體呈下降趨勢;做五年滑動平均后可看出年降水量(圖5)總體呈波動上升趨勢。黃河三角洲入???984—2018年的年輸沙量、年徑流量和年降水量的M-K趨勢檢驗統(tǒng)計值Z值和斜率見表1。年輸沙量和年徑流量的Z值與斜率均小于0,二者呈增大的趨勢,年徑流量的減小速率大于年輸沙量,但年徑流量|Z|<1.96,趨勢不顯著;年降水量的Z值與斜率均>0,呈增加的趨勢,|Z|<1.96,趨勢不顯著。
表1 趨勢性檢驗Table 1 Trend test
2)用M-K檢驗法對年輸沙量、年徑流量和年降水量進(jìn)行趨勢分析和突變檢驗、年輸沙量的統(tǒng)計量UF(圖3b)和年徑流量的UF(圖4b)整體均為負(fù)值,表明二者一直是減少的趨勢,年輸沙量突變發(fā)生在1996年和1998前后、年徑流量突變發(fā)生在1985和2012年前后,兩者突變點的時間存在偏差,這與黃河三角洲的“水沙異源”特點以及一些其他非自然因素的影響有關(guān)[40];年降水量的統(tǒng)計量UF和UB(圖5b)雖然有交點,但均在±1.96之間,變化趨勢不顯著,故無突變點。
圖3 年輸沙量滑動平均圖和M-K檢驗Fig.3 Moving average chart and M-K test of annual sediment discharge
圖4 年徑流量滑動平均圖和M-K檢驗Fig.4 Moving average chart and M-K test of annual runoff
圖5 年降水量滑動平均和M-K檢驗Fig.5 Moving average chart and M-K test of annual precipitation discharge
首先提取出瞬時水邊線(圖6)和一般高潮線(圖7),并基于平均高潮線法計算得出黃河三角洲在1984—2018年的陸地面積,然后對計算結(jié)果開展分析。
圖6 1984—2018年瞬時水邊線變化Fig.6 Changes of instantaneous waterline during 1984—2018
圖7 1984—2018年一般高潮線變化Fig.7 Changes of mean high tide line during 1984—2018
2.2.1 擬合結(jié)果分析
繪制了黃河三角洲總體(圖8a)和各岸段(圖8b)近35 a的面積時間序列變化散點圖,然后通過非參數(shù)方法局部加權(quán)回歸(LOWESS)[41]擬合。從黃河三角洲總面積(圖8a)可知:陸地面積表現(xiàn)出了明顯的先增加、后緩慢減少的趨勢。1984—1997年為黃河三角洲總面積的增長期,在1993年達(dá)到峰值3 443.31 km2,比1984年增加了144.09 km2,河口岸段面積的增長為總面積的增長做出了主要貢獻(xiàn);1997—2018年為黃河三角洲總面積的緩慢減少期,2018年比1993年減少了66.01 km2,減少速率隨時間逐漸降低趨于平緩,主要是由河口岸段面積增長速度小于北部刁口河岸段面積的侵蝕速度所致。
圖8b明顯地反映出黃河三角洲各分區(qū)面積的變化趨勢,其中刁口河岸段面積呈顯著的減少趨勢,河口岸段呈顯著的增加趨勢,1984—1996年大于1996—2018年面積的增長速率;萊州灣岸段和人工岸段沒有顯著的變化趨勢,萊州灣岸段近35 a面積極差為23.11 km2,人工岸段為25.51 km2。
圖8 面積的時間序列變化Fig.8 Time series changes of area
2.2.2 河口岸段面積與水沙量交叉小波和小波相干結(jié)果分析
利用交叉小波和小波相干方法分析黃河三角洲河口岸段面積和水沙量的時間序列在不同時間尺度的變化和相互關(guān)系。交叉小波變換利用2個時間序列的連續(xù)小波變換分析共同高能量區(qū)和相位關(guān)系,小波相干分析2個時間序列在時頻空間中局部的相關(guān)關(guān)系,可以進(jìn)一步解釋在交叉小波能量譜中低能量區(qū)的相關(guān)性。
圖9給出了河口岸段面積與降水量、泥沙量以及徑流量的交叉小波和小波相干功率譜圖。由河口岸段面積和降水量的交叉小波功率譜(圖9a)可知,兩者在1989—1991年存在一年的顯著的共振周期,通過小波相干圖(圖9b)可知在2003—2010年存在周期為5~7 a的顯著的負(fù)相位關(guān)系;由河口岸段面積和泥沙量與徑流量的交叉小波和小波相干功率譜(圖9c~圖9f)可知,面積和泥沙量與徑流量不存在顯著的相關(guān)關(guān)系。
圖9 各因子和河口岸段面積的交叉小波和小波相干功率譜Fig.9 CWT and WTC power spectrum between factors and estuary coastal area
黃河三角洲主要是由黃河挾帶泥沙沖積而形成的粉砂淤泥質(zhì)海岸,其面積穩(wěn)定性較差,受入海水沙、河流改道以及海洋動力條件等因素的影響。從1978—1998年黃河三角洲的面積總體呈現(xiàn)增長的趨勢,隨后直到2016年起黃河三角洲面積開始呈現(xiàn)減少的趨勢[12]。本研究也得到相似的結(jié)果:1984—1997年黃河三角洲總面積呈增長趨勢,且增長速度逐漸變緩;1998—2018年黃河三角洲面積呈緩慢減少趨勢。1976年黃河入??谟稍瓉淼谋辈康罂诤雍涌诟臑槟喜壳逅疁虾涌?這導(dǎo)致了自1976年后刁口河岸段面積急劇減少、現(xiàn)行河口岸段面積增加,是影響黃河三角洲總面積的主要原因。
本文根據(jù)岸段的特點,將整體岸段按照發(fā)展特點分為刁口河岸段、河口岸段、萊州灣段和人工岸段四個區(qū)域。與牛明香[22]劃分的刁口段、東營港及附近岸段、入??诙渭叭R州灣岸段四個區(qū)域相比,雖然岸線提取方法有差異,兩者的數(shù)據(jù)不完全一致,但整體趨勢大致相同。
萊州灣岸段在研究時段內(nèi)沒有顯著的變化趨勢,1984—2007年面積大致穩(wěn)定,2007年后有輕微的減少趨勢,這種變化一方面是由于黃河入海泥沙減少,渤海冷流南下[42]輸運到萊州灣岸段的泥沙量相應(yīng)減少;另一方面由于防潮堤的建設(shè)減緩了岸線的擴張。河口岸段面積受黃河入海泥沙影響最顯著,在研究時間段內(nèi)呈顯著的增長趨勢,按照增長速率可分為1984—1996年和1996—2018年兩個時間段。由于勝利油田孤東圍堤之間有穩(wěn)固的防潮堤,所以該人工岸段相對穩(wěn)定。刁口河岸段從1976年黃河改道后一直呈蝕退趨勢,雖然人工固化海岸逐漸增多,但仍然不及海水動力作用的侵蝕速度,故總體面積呈減少趨勢。
因為水沙量對現(xiàn)行河口岸段的影響最大,故單獨選出河口岸段分析其面積的變化與利津站年輸沙量、年徑流量和黃河三角洲年降水量的關(guān)系。通過河口岸段面積的時間序列圖(圖8b)可知,1984—1996年河口岸段面積呈顯著增長的趨勢,同時期的輸沙量和徑流量通過5 a滑動平均(圖3a、圖4a)亦可得到為1984—2018年的最大值時間段,年降水量的5 a滑動平均也為增長的趨勢,說明1984—1996年水沙量是河口面積增長的重要驅(qū)動力;1997—2003年河口岸段面積呈緩慢的減少趨勢,這也與水沙量有密切的關(guān)系。由圖3b,圖4b和圖5b可知年輸沙量突變發(fā)生在1997年,年輸沙量、年徑流量和年降水量1997—2002年均為下降的趨勢。根據(jù)利津水文站的觀測數(shù)據(jù)可知,1996年黃河改道清8汊后,原清水溝流路失去水沙的供應(yīng),1997年黃河斷流202 d,岸線不斷被侵蝕,一直到2004年黃河處于嚴(yán)重的斷流狀態(tài),導(dǎo)致清8汊流路的造陸面積變小,河口岸段面積整體呈減少趨勢。2004—2011年由于調(diào)水調(diào)沙工程河口岸段面積呈增加趨勢;2011—2018年北部沙嘴朝向(圖6c和圖7c)逐漸由東向變?yōu)楸毕?河口岸段面積呈減少趨勢。此時間段輸沙量、徑流量和降水量與面積的變化關(guān)系不顯著,可能是因為黃河2002年調(diào)水調(diào)沙之后,入海水沙模式與之前不同,具體原因尚需進(jìn)一步探討。
通過交叉小波和小波相干分析了河口岸段面積與水沙量在不同時間尺度的變化和相互關(guān)系,由河口岸段面積與泥沙徑流量交叉小波譜和小波相干譜(圖9c~圖9f)可知,面積與泥沙徑流不存在顯著的相關(guān)關(guān)系,這進(jìn)一步佐證了河口岸段的面積不僅受入海水沙的影響,還受潮汐、海流、泥沙粒徑及調(diào)水調(diào)沙工程等因素的影響,未通過顯著性檢驗的低能量區(qū)均以負(fù)相位為主,也可以說明河口面積的變化滯后于水沙的變化;由小波交叉譜圖(圖9a)可知,河口岸段面積與降水量只在1989—1991年存在1 a的顯著共振周期,其他未通過顯著性檢驗的低能量區(qū)均以負(fù)相位為主,由小波相干譜(圖9b)可知在2003—2010年存在周期為5~7 a的顯著的負(fù)相位關(guān)系,說明河口岸段面積變化滯后于黃河三角洲年降水量,這與黃河入海泥沙對黃河三角洲造陸存在明顯的“后效”影響[22]一致。
基于Landsat衛(wèi)星數(shù)據(jù),利用遙感大數(shù)據(jù)平臺Google Earth Engine和GIS技術(shù)分析了黃河三角洲地區(qū)岸線變遷趨勢,并研究了河口岸段面積與黃河三角洲年輸沙量、年徑流量、年降水量的關(guān)系。得到主要結(jié)論:
1)1984—2018年黃河三角洲總面積先增加后緩慢減少,變化主要發(fā)生在北部刁口河岸段和現(xiàn)行河口岸段,因人工固化和防潮堤的修建,人工岸段和萊州灣岸段面積保持穩(wěn)定,對黃河三角洲總面積的影響較小。
2)輸沙量和徑流量總體趨勢基本一致,豐枯交替,整體呈下降趨勢。河口岸段面積與水沙量存在密切的關(guān)系,并且輸沙量、徑流量和降水量對面積的變化存在明顯的“后效”影響,比面積變化提前一年。此外,海洋中的潮汐和海流等海水動力作用以及人工固岸和調(diào)水調(diào)沙等非自然因素也對河口岸段的變化起著至關(guān)重要的影響。
本文使用的岸線提取方法相對方便快捷,具有一定的普適性,平均高潮線法提取岸線受潮汐影響較小,能夠滿足宏觀分析的精度。LandSat系列衛(wèi)星的TM、ETM+和OLI數(shù)據(jù)量大,可以篩選出少云清晰的影像,而GEE平臺的高性能計算能力使分析處理數(shù)據(jù)更高效,特別是在對全球各地岸線的大尺度和長時間序列分析。所以,此方法未來可在其他河口地區(qū)進(jìn)行試驗。