1979—2012年萊州灣南岸海水入侵與區(qū)域海岸線變動(dòng)時(shí)空耦合分析

蒙永輝， 王集寧， 張麗霞， 羅 梅

(山東省地質(zhì)環(huán)境監(jiān)測(cè)總站,濟(jì)南 250014)

0 引言

海水入侵問題是近海岸帶地區(qū)地質(zhì)環(huán)境災(zāi)害研究領(lǐng)域中的一個(gè)熱點(diǎn)，也是國際學(xué)術(shù)界一直關(guān)心的一個(gè)焦點(diǎn)環(huán)境問題[1-6]。海水入侵是指在人為和自然因素的干擾下，濱海地區(qū)海水和陸地地下淡水的平衡狀態(tài)遭到破壞，導(dǎo)致海水不斷向內(nèi)陸方向侵入的過程和現(xiàn)象[7-9]。雖然海水的密度(礦化度)高于淡水，但是陸地淡水的水位高于海水，所以濱海地區(qū)的地下淡水和聯(lián)通海水通過壓力抗衡，在理論上能夠維持一個(gè)相對(duì)平衡穩(wěn)定的界面。但是，這種穩(wěn)定一旦受到外界的干擾，例如： 氣候變化、淡水資源的過渡開采、海平面上升等，平衡就會(huì)被打破，海水就會(huì)不斷地侵入內(nèi)陸淡水層。大量的海水向陸地倒灌，就會(huì)嚴(yán)重破壞區(qū)域生態(tài)環(huán)境，導(dǎo)致原本可以利用的淡水資源被污染，地下水水質(zhì)變咸，土壤發(fā)生鹽堿化，生產(chǎn)和生活用水困難，對(duì)區(qū)域社會(huì)發(fā)展造成諸多負(fù)面影響并引起一些環(huán)境問題，如淡水水資源匱乏、耕地資源退化、土壤次生鹽堿化、原生植被凋亡等[10-14]。

位于山東省濰坊市北部的萊州灣南岸海岸帶地區(qū)是我國最早發(fā)現(xiàn)海水入侵現(xiàn)象的地區(qū)之一[9,15]，也是我國目前海水入侵最為嚴(yán)重的地區(qū)[8,16]，已引起了當(dāng)?shù)卣块T和研究人員的高度重視。山東省人民政府早在1990年就出臺(tái)了山東省萊州灣地區(qū)海水侵染災(zāi)情分析與綜合治理方案，1994年出版了山東省萊州灣地區(qū)海水侵染綜合治理規(guī)劃[17]。相關(guān)領(lǐng)域科研工作者在海水入侵特征、成因、過程、模擬和評(píng)價(jià)方面開展了大量研究，取得了豐富的成果[3-5,18-27]。在海水入侵過程與特征研究方面，李新運(yùn)等[20]對(duì)萊州灣南岸的海水入侵過程進(jìn)行了分析和預(yù)測(cè)； 孫云華等[23]分析了1979―2008年海岸地貌過程演變與人類活動(dòng)地下海水入侵的關(guān)系； 蒙永輝等[24]研究分析了濰坊北部海咸水入侵的特征和現(xiàn)狀。在海水入侵成因與對(duì)策研究方面，劉桂儀[25]探討了萊州灣南岸的海水入侵原因和防治對(duì)策； 豐愛平等[21]分析了萊州灣南岸1980―2000年海岸侵蝕過程并分析了海水入侵原因； 孫云華等[23]以1979、1989、2001和2008年4個(gè)時(shí)相的遙感影像為數(shù)據(jù)源，采用景觀空間分析方法，分析了人類活動(dòng)對(duì)海岸帶地貌過程及海水入侵的影響。在海水入侵的模型模擬與評(píng)價(jià)研究方面，李福林[26]對(duì)萊州灣東岸海水入侵進(jìn)行了動(dòng)態(tài)監(jiān)測(cè)與數(shù)值模擬研究； 蘇喬等[27]對(duì)萊州灣南岸海水入侵現(xiàn)狀進(jìn)行評(píng)價(jià)； 陳廣泉等[22]評(píng)價(jià)分析了萊州灣地區(qū)海水入侵災(zāi)害風(fēng)險(xiǎn)。綜上所述，對(duì)萊州灣地區(qū)海水入侵研究的成果相當(dāng)豐碩，但是對(duì)海水入侵特征過程仍缺乏定量系統(tǒng)的研究，對(duì)最新的海水入侵動(dòng)態(tài)關(guān)注較少； 對(duì)海水入侵變化原因多為定性描述，缺少定量探討； 對(duì)海水入侵年際時(shí)空變化規(guī)律并不明確[8]。

本文在已有萊州灣海水入侵監(jiān)測(cè)資料的基礎(chǔ)上，結(jié)合現(xiàn)代遙感觀測(cè)手段和空間分析技術(shù)，對(duì)1979―2012年萊州灣南岸海水入侵發(fā)生與發(fā)展動(dòng)態(tài)演變過程進(jìn)行數(shù)字重建，根據(jù)相應(yīng)時(shí)間段的海岸線信息，分析其時(shí)空演化特征，重點(diǎn)探討和分析區(qū)域海水入侵過程與海岸線變遷的時(shí)空耦合關(guān)系。

1 研究區(qū)與數(shù)據(jù)源

1.1 研究區(qū)概況

本文的研究區(qū)位于萊州灣南岸濱海地區(qū)(圖1)，地處膠東半島西北部，瀕臨渤海，地理坐標(biāo)在118°40′～119°40′E， 36°40′～37°20′N之間，涉及壽光、濰坊(寒亭)和昌邑3個(gè)市。氣候?qū)倥瘻貛Ъ撅L(fēng)型大陸性氣候，干旱少雨，年內(nèi)降雨量不均。地質(zhì)構(gòu)造復(fù)雜、地貌類型多樣，主要為由河流沖洪積與海水作用形成的沖海積平原，廣泛分布著各類第四紀(jì)沉積物和粉沙質(zhì)海岸地貌[9]。其中第四系沉積層厚度較大，平均200～300 m，沉積物成因類型以海相、湖沼相、河流沖積物為主； 入海河流形成的古河道帶透水性較強(qiáng)，成為咸水入侵的良好通道[24]。該地區(qū)資源豐富、人口密集、海洋經(jīng)濟(jì)發(fā)達(dá)，擁有得天獨(dú)厚的鹵水礦資源，是我國沿海最大的鹵水礦區(qū)，但也是我國最早發(fā)現(xiàn)海水入侵現(xiàn)象的地區(qū)之一[16]。受氣候變化和人類活動(dòng)的共同影響，該區(qū)域已經(jīng)成為我國目前海水入侵最嚴(yán)重的地區(qū)，給當(dāng)?shù)毓ぁ?農(nóng)業(yè)生產(chǎn)和人民生活造成了極大危害。

圖1 研究區(qū)位置 (影像底圖為TM B4(R),B3(G),B2(B)假彩色合成影像)Fig.1 Location of study area

1.2 數(shù)據(jù)源

本文中使用的海水入侵?jǐn)?shù)據(jù)主要根據(jù)魯勘字[2011]14號(hào)項(xiàng)目《黃河三角洲高效生態(tài)經(jīng)濟(jì)區(qū)(濰坊)海(咸)水入侵調(diào)查與監(jiān)控預(yù)警系統(tǒng)建設(shè)報(bào)告》和萊水字[2011]39號(hào)文件，以及相關(guān)文獻(xiàn)中的數(shù)據(jù)綜合整理得到[17,21,23]。對(duì)海水入侵?jǐn)?shù)據(jù)整理、利用中，存在空間和尺度的誤差，需要進(jìn)行空間配準(zhǔn)和一致性檢驗(yàn)，以去除彼此矛盾的數(shù)據(jù)。根據(jù)這些數(shù)據(jù)的采集時(shí)間，同時(shí)考慮遙感數(shù)據(jù)質(zhì)量、受云層干擾情況等，去匹配相應(yīng)時(shí)間的衛(wèi)星影像。最終，本文挑選6期(1979, 1990,1995,2000,2003和2013年)遙感數(shù)據(jù)用于監(jiān)測(cè)岸線變化(表1)。對(duì)遙感數(shù)據(jù)進(jìn)行了預(yù)處理，將DN值轉(zhuǎn)換為光譜反射率[28-29]。

表1 用于海岸線監(jiān)測(cè)的遙感數(shù)據(jù)Tab.1 Remote sensing data used for shoreline monitoring

2 研究方法

本文的總體技術(shù)流程如圖2所示。

圖2 技術(shù)流程圖Fig.2 Technical flow chart

在多時(shí)相海水入侵?jǐn)?shù)據(jù)整理和相應(yīng)時(shí)段遙感影像選擇的基礎(chǔ)上，首先對(duì)海水入侵專題圖進(jìn)行配準(zhǔn)和數(shù)字化，整合1979―2012年6期海水入侵鋒線； 對(duì)遙感影像進(jìn)行預(yù)處理，提取1979―2013年的海岸線位置分布信息。然后，采用縱向剖面分析(垂直斷面建模方法)[28]，每隔200 m采樣，獲取每個(gè)空間采樣點(diǎn)的定量變化信息； 再根據(jù)端點(diǎn)速率法(end point ratio，EPR)分析模型，分別計(jì)算和統(tǒng)計(jì)分析入侵鋒線、海岸線的時(shí)間分段變化速率及多年平均變化率; 最后，通過Pearson相關(guān)分析與顯著性t檢驗(yàn)，分析海岸線變化與海水入侵在時(shí)間和空間上的耦合關(guān)系。

2.1 海岸線提取與鋒線數(shù)字化

通過遙感影像監(jiān)測(cè)海岸線的變化，在國內(nèi)外已有大量的算法研究和應(yīng)用。本文主要利用面向?qū)ο蠓诸惣夹g(shù)，根據(jù)岸線在遙感影像中的空間分布特征和水陸光譜差異，依據(jù)歸一化差值水體指數(shù)(normalized difference water index，NDWI)類間光譜特征差異最大原理，完成水陸分離； 然后通過人工互信息操作，實(shí)現(xiàn)海岸線的高精度識(shí)別，具體技術(shù)流程參考文獻(xiàn)[28]。首先采用大津算法(OTSU)完成NDWI水陸分離； 然后依據(jù)空間關(guān)系特征判斷實(shí)現(xiàn)陸地水域與海洋判別； 最后通過互信息人工操作對(duì)海岸線提取結(jié)果進(jìn)行優(yōu)化處理和微地形修正，得到研究區(qū)的各期海岸線分布圖?；バ畔⒑筇幚聿僮髦饕ǎ?河口處理需保留大型河口港灣特征和小河齊陸地基線； 瀉湖需考慮其內(nèi)環(huán)境，當(dāng)瀉湖與海洋有寬闊水域通道時(shí)，在瀉湖內(nèi)部采集海岸線。海水入侵鋒線是在幾何精確配準(zhǔn)的基礎(chǔ)上，對(duì)歷史時(shí)期的海水入侵鋒線進(jìn)行手工數(shù)字化，得到各個(gè)歷史時(shí)期海水入侵鋒線的空間分布位置。

需要說明的是，本文中的海岸線數(shù)據(jù)沒有考慮平均大潮位修正，不是真正意義的海岸線； 準(zhǔn)確地說是利用遙感影像提取的瞬時(shí)水邊線。理論上做潮位修正效果可能會(huì)更好，但因這個(gè)區(qū)域地表變化很大，SRTM高程數(shù)據(jù)在這個(gè)區(qū)域不是特別準(zhǔn)確，潮位修正后的結(jié)果和遙感影像套合不上。所以，綜合考慮到本文中分析數(shù)據(jù)是通過EPR 模型計(jì)算多個(gè)時(shí)期時(shí)間序列岸線的平均變化速率，盡管用水邊線代替海岸線時(shí)潮位會(huì)對(duì)分析結(jié)果會(huì)產(chǎn)生一定的影響(不可避免)，但基本可以近似地反映區(qū)域海岸線的動(dòng)態(tài)變化趨勢(shì)(與其他文獻(xiàn)中描述的趨勢(shì)基本一致)，具有一定的可比性，對(duì)分析海水入侵鋒線位置影響不大。

2.2 縱向斷面建模方法

垂直斷面建模方法是一種對(duì)數(shù)字岸線線性目標(biāo)變遷速率計(jì)算傳統(tǒng)的方法，其中以端點(diǎn)速率法 (EPR)計(jì)算模型最為簡單常用[28]。EPR模型統(tǒng)計(jì)計(jì)算岸線的變化速率原理如圖3所示。

圖3 EPR局部垂直(縱向)斷面建模計(jì)算示意圖Fig.3 Sketch map of EPR model

圖中P1，P2，P3，…,Pn基于基線的局部垂直采樣點(diǎn)(n的大小由采樣密度決定)；timei是獲取時(shí)間；D1，D2，D3，…，Di是垂線(縱向剖面線)與變遷線的交點(diǎn)；distance(Di-Dj)是Di點(diǎn)與Dj點(diǎn)之間的距離。

首先建立垂直剖面線，本文在提取海水入侵鋒線和岸線的基礎(chǔ)上，通過緩沖區(qū)分析建立基線，設(shè)置采樣間隔為200 m； 然后根據(jù)垂線與所有斷面上的岸線交點(diǎn)距離，計(jì)算平均變化率。那么，對(duì)于2個(gè)歷史時(shí)期的岸線變化速率，可直接用2個(gè)歷史時(shí)期的變遷距離除以時(shí)間來計(jì)算，EPR1,2的數(shù)學(xué)表達(dá)式為

(1)

式中：D1，D2分別為time1，time2兩個(gè)歷史時(shí)期岸線的位置；distance(D1-D2) 為D1，D2兩點(diǎn)之間的距離。

對(duì)于多個(gè)時(shí)期的一組時(shí)間序列岸線的平均變化速率，EPRave的數(shù)學(xué)表達(dá)式為

(2)

式中：timelatest為最新歷史時(shí)間；timeoldest為最老歷史時(shí)間；Dlatest,Doldest分別為timelatest,timeoldest兩個(gè)歷史時(shí)期岸線的位置。

3 結(jié)果分析

3.1 海水入侵鋒線動(dòng)態(tài)變化過程

1979―2012年萊州灣海岸帶區(qū)域海水入侵鋒線位置動(dòng)態(tài)變化見圖4。從咸水入侵鋒線圖可以清楚地看出區(qū)域海水入侵的動(dòng)態(tài)變化過程。1979―1989年是區(qū)域海水入侵最快的階段，隨后區(qū)域海水鋒線移動(dòng)明顯變緩。圖5為海水入侵EPR時(shí)空動(dòng)態(tài)分布圖。綜合分析海水入侵的動(dòng)態(tài)過程和EPR時(shí)空動(dòng)態(tài)分布圖，可以從總體上將區(qū)域海水入侵分為3個(gè)階段： ①1979―1989年的快速入侵階段。此階段是萊州灣南岸海水入侵發(fā)展最快的時(shí)期，每年入侵面積增加29.22 km2； 海咸水平均每年入侵速率為380 m/a，最大入侵距離在昌邑市青鄉(xiāng)南部，入侵速率達(dá)到690 m/a； ②1989―2000年的慢速發(fā)展階段。此階段萊州灣南岸海水入侵的速度明顯變緩，無論是入侵速度，還是入侵面積，都明顯下降。海水入侵鋒線基本上在1995年的位置附近波動(dòng)； ③2000―2012 年的相對(duì)穩(wěn)定與徘徊階段。此階段海咸水入侵動(dòng)態(tài)鋒線基本處于徘徊狀態(tài)，并且局部出現(xiàn)了退縮(圖5(e))。

圖4 1979―2012年萊州灣南岸海水入侵鋒線動(dòng)態(tài)過程 (影像底圖為TM B3(R)，B2(G)，B1(B)假彩色合成影像)Fig.4 Processes of seawater intrusion strike in south coast of Laizhou Bay from 1979 to 2012

(a) 1978—1989年海水入侵速率 (b) 1989—1995年海水入侵速率(c) 1995—2000年海水入侵速率

(d) 2000—2008年海水入侵速率(e) 2008—2012年海水入侵速率(f) 1979—2012年海水入侵速率

圖5EPR空間動(dòng)態(tài)分布圖

Fig.5SpatialdynamicdistributionofEPR

3.2 海岸線動(dòng)態(tài)變化遙感制圖與分析

從解譯的結(jié)果可以看出，萊州灣南部岸線的時(shí)空變化并不是一致的(圖6)。

圖6 1979―2013年區(qū)域海岸線變遷圖 (影像底圖為MSS B3(R),B2(G),B1(B)假彩色合成影像)Fig.6 Regional coastline transition from 1979 to 2013

從空間上看，區(qū)域海岸線總體上呈現(xiàn)向陸地蝕退的趨勢(shì)(局部區(qū)域受到人類活動(dòng)的干擾(如人工修建海上設(shè)施)除外)； 從時(shí)間上看，不同地點(diǎn)區(qū)域海岸線的蝕退速率也不一樣，相同地點(diǎn)不同時(shí)段的蝕退速率也不同。

3.3 海水入侵與岸線變化耦合關(guān)系分析

1979―2012年萊州灣南岸的海水入侵與岸線變化速率見圖7。從圖7可以看出2條EPR曲線的變化特征和時(shí)空耦合關(guān)系。在空間分布格局上，海水入侵的強(qiáng)度和海岸線的蝕退強(qiáng)度相關(guān)性具有很強(qiáng)的一致性； 在時(shí)間變化速率上，海水的入侵速率變化曲線和海岸線變化的變化速率曲線形狀相似。這一變化特征說明，萊州灣南岸地區(qū)海水入侵和海岸線蝕退在空間分布上和時(shí)間變化上均具有一定的相關(guān)性。

(圖中Y軸表示基準(zhǔn)岸線，正值表示向海方向，負(fù)值表示向陸地方向)圖7 1979―2012年萊州灣南岸岸線與海水入侵鋒線變化速率圖Fig.7 Change rates of coastline and seawater intrusion front line of Laizhou Bay from 1979 to 2012

表2從定量的角度統(tǒng)計(jì)了區(qū)域岸線蝕退和海水入侵的變化特征以及二者之間的相關(guān)關(guān)系。表中海岸蝕退正值表示向海方向淤進(jìn)，負(fù)值表示向陸地方向蝕退； 海水入侵負(fù)值表示向陸地方向侵入； ** 表示在0.01水平(雙側(cè))顯著相關(guān)。

表2 岸線蝕退與海水入侵變化量統(tǒng)計(jì)與耦合關(guān)系分析Tab.2 Relationship analysis between seawater intrusion and coastline changes (m/a)

1979―2012年萊州灣南岸海水入侵速率均值達(dá)到了 177.23 m/a； 最大值分布在昌邑市青鄉(xiāng)南部區(qū)域，入侵速率達(dá)到 435.28 m/a； 最小值分布在濰坊北部寒亭,為 124.02 m/a。1979―2012年區(qū)域岸線的蝕退速率約為10.44 m/a； 蝕退的最大速率為124.02 m/a，空間位置上分布在昌邑市北部下營鎮(zhèn)； 蝕退的最小速率分布在寒亭區(qū)以北的岸線段。此處由于人工建造海上游樂場和港口，區(qū)域岸線呈現(xiàn)向海洋快速淤進(jìn)，最大速度達(dá)到191.14 m/a。海岸線EPR與入侵鋒線EPR相關(guān)系數(shù)為0.407，顯著性水平P值為0.00，通過99%的信度水平(雙側(cè))檢驗(yàn)。由此可見，區(qū)域海水入侵與岸線蝕退在時(shí)空上存在強(qiáng)耦合關(guān)系。

4 結(jié)論與討論

4.1 結(jié)論

本文對(duì)1979―2012年萊州灣南岸海水入侵鋒線演化的時(shí)空過程進(jìn)行數(shù)字重建和區(qū)域海岸線的多期遙感監(jiān)測(cè)； 通過EPR定量模型探索性地分析了二者之間的時(shí)空耦合關(guān)系機(jī)制。研究得到以下結(jié)論：

1)區(qū)域海水入侵經(jīng)歷了由快到慢的變化過程，1990年以后速率明顯減緩，入侵鋒線基本穩(wěn)定在1995年鋒線附近； 2008―2012年入侵鋒線局部發(fā)生后退，目前保持相對(duì)穩(wěn)定。這一研究結(jié)論與之前相關(guān)學(xué)者的研究結(jié)果基本一致[22]。

2)除了局部人工造陸導(dǎo)致岸線向海擴(kuò)張外，1979―2012年區(qū)域海岸線以蝕退型海岸為主，但是蝕退速度在空間上差異明顯。

3)海水入侵鋒線變化與海岸線進(jìn)退二者之間在時(shí)空上存在強(qiáng)耦合關(guān)系，相關(guān)系數(shù)達(dá)到0.407，顯著性水平P<0.01(雙側(cè))，說明海岸線工程對(duì)區(qū)域海水入侵速率有顯著的影響。以上研究結(jié)論可為區(qū)域海水入侵的預(yù)防和治理提供數(shù)據(jù)支撐和科學(xué)依據(jù)。

4.2 討論

1)在變化原因上，影響咸淡水界面遷移速率的因素可分為自然因素和人為因素，自然因素主要是地質(zhì)條件、大氣降水減少和風(fēng)暴潮等； 人為因素主要是地下水開采及鹵水的開采。特殊的地質(zhì)條件(廣泛分布的第四紀(jì)海相、湖相沉積以及河流沖積物，透水性強(qiáng))為研究區(qū)內(nèi)海咸水入侵提供了有利條件； 而氣候變暖及海平面升高是該區(qū)海咸水入侵的大環(huán)境背景，人類活動(dòng)則改變了海水入侵時(shí)空變化進(jìn)程和入侵強(qiáng)度[9,24]。

2)在變化特征上，1979―1989年海咸水快速入侵階段，主要是由于地下淡水資源的過渡開采，地下水位大幅度下降； 1989―2000年海咸水慢速發(fā)展階段，由于當(dāng)?shù)刂饾u認(rèn)識(shí)到海咸水入侵的危害，對(duì)地下水開采進(jìn)行了適當(dāng)限制，同時(shí)近海鹵水的開采量在逐年增大； 2000―2012 年海咸水相對(duì)穩(wěn)定與后退階段，主要是得益于海水入侵治理工程得實(shí)施，通過協(xié)調(diào)和控制地下水與鹵水的開采，較好地控制了咸淡水界面。因此，按照研究區(qū)咸淡水界面變化情況，大量的研究分析共同表明： 地下水資源開發(fā)和鹵水資源開采起了決定性作用，是區(qū)域海水入侵的直接原因[3,19, 22-23]。

3)海岸線的進(jìn)退是影響區(qū)域海水入侵的又一疊加因素。海岸線的淤進(jìn)，尤其是人工修筑的岸線，在一定程度上可以阻止海水入侵。可能的原因是人工修筑的岸線有利于阻擋風(fēng)暴潮，堵塞海水入侵的部分通道。

以上研究結(jié)論建立在文中數(shù)據(jù)和方法的基礎(chǔ)上。由于EPR模型方法計(jì)算線性目標(biāo)的變化率是通過剖面上兩點(diǎn)距離與時(shí)間的比值得到的，對(duì)于短期的加速、緩慢或者逆轉(zhuǎn)變化不夠敏感。尤其是在計(jì)算多期平均變化率時(shí)，短期的逆轉(zhuǎn)性變化對(duì)整體變化趨勢(shì)會(huì)產(chǎn)生一定的影響。此外，多源數(shù)據(jù)的空間精度誤差問題、時(shí)間尺度問題不可回避。在空間尺度上，遙感數(shù)據(jù)源空間分辨率為30～60 m，在提取的海岸線精度上難免有些誤差； 海水入侵歷史數(shù)據(jù)資料的幾何變形，導(dǎo)致局部精度難以控制。在時(shí)間尺度上，受到觀測(cè)數(shù)據(jù)的限制，通過6期數(shù)據(jù)分析變化過程，時(shí)間分辨率精度不夠高，對(duì)演化進(jìn)程分析尤其是時(shí)間拐點(diǎn)的確定不夠準(zhǔn)確。這些可能的誤差會(huì)對(duì)數(shù)據(jù)分析以及部分研究結(jié)論產(chǎn)生一定的影響。在未來的研究中有些方面將有待于進(jìn)一步改進(jìn)和提高，例如加大海水入侵鋒線動(dòng)態(tài)監(jiān)測(cè)時(shí)間密度，改進(jìn)算法提高對(duì)短期動(dòng)態(tài)變化的敏感性等。