易雨君，謝泓毅，宋 劼，楊志峰，4

（1.北京師范大學(xué)教育部水沙科學(xué)重點實驗室，北京 100875；2.北京師范大學(xué)水環(huán)境模擬國家重點實驗室，北京 100875；3.交通運輸部規(guī)劃研究院環(huán)境資源所，北京 100028；4.廣東工業(yè)大學(xué)環(huán)境生態(tài)工程研究院，廣東廣州 510006）

1 研究背景

鹽沼植被適宜生境保護是維持黃河三角洲生物多樣性的重要環(huán)節(jié)之一［1-2］。許多研究表明，鹽沼植被的演替與鹽沼的海拔高度和水文過程有關(guān)，尤其與地下水過程關(guān)系緊密［3］。不同地下水流動模式下，非生物因子梯度，如土壤鹽度與含水率梯度在空間上存在明顯差異，這種差異直接導(dǎo)致了鹽沼植被適宜生境的變化［4］。對于鹽沼植被適宜生境的模擬始于偏好曲線法，該方法基于生物對生境因子的偏好，構(gòu)建環(huán)境因子與生物偏好之間的關(guān)系，得到物種生境的適宜程度［5］。但該方法沒有考慮環(huán)境因子之間的相互作用，不能適用于環(huán)境因子相互作用關(guān)系復(fù)雜的情況。

為了解決偏好曲線法的局限，一些研究開始采用多元統(tǒng)計方法，通過建立目標物種分布與環(huán)境因子之間的統(tǒng)計關(guān)系來模擬適宜生境，包括多元線性回歸法、嶺回歸法、主成分回歸法、邏輯回歸法、廣義線性模型和廣義可加模型等［6-7］。偏好曲線法和多元統(tǒng)計方法能較好地描述物種對多環(huán)境因子的響應(yīng)問題，并和水動力水質(zhì)模型相結(jié)合，得到適宜生境范圍的時空演變過程。但上述方法主要基于環(huán)境因子和生物選擇的統(tǒng)計關(guān)系［8］，無法考慮生物的動力過程。有研究提出了基于生物生長動力過程的適宜生境模擬方法［9］，該方法能夠模擬環(huán)境因子和生物過程以及兩者之間的動態(tài)作用關(guān)系［10］。目前基于生物生長動力過程的適宜生境模型主要有元胞自動機法（Cellular Automata，CA），該方法能模擬物種的生長和擴散過程［11］；同時，有研究從生態(tài)學(xué)角度出發(fā)，提出了模擬種群增長和種間競爭的模型［12］。然而，模擬種群擴散和增長的元胞自動機模型，和模擬種間競爭的模型，往往基于靜態(tài)的環(huán)境因子空間分布，或?qū)⒖臻g視為均質(zhì)統(tǒng)一體，無法體現(xiàn)環(huán)境因子以及生物分布在空間和時間上的異質(zhì)性。因此，針對目前生境模擬中面臨的環(huán)境因子時空異質(zhì)性與生物生長、擴散以及種間競爭之間無法兼顧這一問題，需提出新的生境適宜度模擬方法。

鹽沼濕地地下水埋深較淺，地下水過程對鹽沼植被根系所在的淺層土壤影響顯著而迅速，地下水過程時空變化特征對鹽沼植被生境適宜度的分布狀況會有顯著影響［13］。目前已有一些研究表明地下水動力三維數(shù)值模擬應(yīng)用廣泛，對于濱海鹽沼濕地地下水模擬也較為適用［14］。

本文基于MODFLOW 模型構(gòu)建潮汐和徑流共同作用下地下水動力模型，并構(gòu)建土壤鹽度模型模擬鹽度的時空變化；將元胞自動機模型和改進的Logistic模型耦合，構(gòu)建綜合考慮鹽沼植被生長、擴散及種間競爭的鹽沼植被種群增長-競爭動力學(xué)模型，其中，元胞自動機模型模擬植被的空間擴散，改進的Logistic 模型模擬鹽沼植被的生物量積累。將地下水動力模型、土壤鹽度模型與植被種群增長-競爭動力學(xué)模型相結(jié)合，模擬不同徑流過程和潮汐共同作用下，鹽沼植被群落的適宜生境及生物量的時空分布。

2 研究區(qū)及監(jiān)測實驗設(shè)計

2.1 研究區(qū)及監(jiān)測樣點布設(shè)黃河三角洲國家級自然保護區(qū)位于現(xiàn)行黃河流路和黃河故道的入?？谔?，總面積達到15.3×104hm2，分為實驗區(qū)、過渡區(qū)和核心區(qū)。黃河自西向東穿過研究區(qū)，最終匯入渤海，是研究區(qū)最主要的淡水源，也是區(qū)域地下水的主要補給源。淺層地下水的埋深非常淺，約為0～3 m，以微咸水、咸水和鹵水為主。地下水位較高，一般在-2.5～7.5 m。自海向陸環(huán)境梯度上，日本鰻草-互花米草-翅堿蓬-檉柳-蘆葦?shù)塞}沼植被呈典型帶狀分布。

在黃河口自然保護區(qū)內(nèi)，植被具有典型帶狀分布格局，依據(jù)鹽沼植被分布分成5個區(qū)，依次為裸灘區(qū)（W1）、翅堿蓬區(qū)（W3）、交錯區(qū)（W5）、檉柳區(qū)（W7）和蘆葦區(qū)（W10），如圖1所示。選取受人為干擾較小的區(qū)域，設(shè)置長約為4 km的樣帶，沿樣帶設(shè)置5個監(jiān)測樣點。2019年3—10月定期（每月25日）監(jiān)測植被生長狀態(tài)、土壤及地下水?dāng)?shù)據(jù)，并采集表層土壤樣品進行理化性質(zhì)分析，于2019年7月2—29日黃河調(diào)水調(diào)沙期間開展了為期28 d的持續(xù)原位監(jiān)測，每個監(jiān)測點設(shè)置兩個不同開篩深度的地下水監(jiān)測井，井深3 m，開篩深度分別為0.5～1 m 和2～2.5 m（圖1（c））。同步監(jiān)測土壤鹽度、含水率、地下水埋深、黃河水位和渤海潮高等指標，持續(xù)監(jiān)測時長約700 h。

2.2 原位監(jiān)測實驗于2019年7月2—29日期間，每隔3日采集1次淺層土壤（0～40 cm）樣品，并測定其鹽度和含水率。每個監(jiān)測點取3個平行樣，以平行樣均值作為測定值。淺層土壤含水率采用烘干稱重法測定，鹽度采用土壤浸出液法［15］測定。

淺層土壤鹽度和含水率的監(jiān)測結(jié)果（圖2）表明，裸灘區(qū)、翅堿蓬區(qū)和交錯區(qū)（W1、W3、W5）淺層土壤鹽度和含水率較為接近，與檉柳區(qū)（W7）相比，上述區(qū)域淺層土壤鹽度均明顯偏低，而含水率均明顯偏高。由海向陸至黃河梯度上，淺層土壤鹽度和含水率的變化范圍增大。近海的翅堿蓬區(qū)土壤鹽度較低而含水率較高，變化范圍均較??；與之相反，靠近黃河的檉柳區(qū)土壤鹽度較高，含水率較低，變化范圍均較大。位于兩者之間的交錯區(qū)與翅堿蓬區(qū)的鹽度和含水率較為一致，但變化范圍顯著變大。交錯區(qū)作為過渡區(qū)，其含水率和鹽度為翅堿蓬適宜范圍，又因含水率和鹽度的變化范圍大，為檉柳提供了適宜生境［16］。蘆葦區(qū)（W10）鹽度明顯偏低，含水率的均值和變化范圍最大。

圖2 淺層土壤含水率和鹽度的監(jiān)測結(jié)果

開展原位監(jiān)測實驗的同時，在潮灘上和黃河岸邊分別設(shè)定了潮位和黃河水位監(jiān)測點，使用水位記錄儀分別記錄實時潮高和黃河水位。潮汐和黃河水位的監(jiān)測結(jié)果（見圖3）表明，樣帶北側(cè)存在大小潮持續(xù)影響，潮汐周期約為14 d，大潮出現(xiàn)在7月6日和18日前后，小潮出現(xiàn)在7月12日和26日前后；樣帶南側(cè)黃河水位受到人為擾動較大，水位波動劇烈。在調(diào)水調(diào)沙前（7月2—9日）、中（7月10—16日）和后期（7月17—29日）變化顯著，調(diào)水調(diào)沙中期黃河水位激增，水位波動較大，調(diào)水調(diào)沙前期和后期水位較低。

圖3 潮高和黃河水位變化監(jiān)測結(jié)果

3 模型構(gòu)建

3.1 地下水動力模型采用MODFLOW模型模擬地下水動力條件?；谶_西定律和滲流連續(xù)方程和實際水文地質(zhì)條件，建立淺層地下水三維非穩(wěn)定流數(shù)學(xué)模型［17］，公式如下：

式中：Ω為滲流模擬區(qū)域；h為含水層水頭，m；Kx、Ky、Kz分別為x、y、z軸的滲透系數(shù)，m/d；μd為給水度，1/m；ε為源匯項，1/d；h0為初始水頭，m；h1為第一類邊界的水頭，m；t為時間，d；S1和S2代表第一類和第二類邊界。

數(shù)學(xué)模型的求解采用有限差分法，強隱式迭代法解算器（Strongly Implicit Procedure，SIP）求解，最大迭代次數(shù)為50，收斂的判別標準和殘差判別標準為0.001，網(wǎng)格考慮干濕轉(zhuǎn)換。

3.2 土壤鹽度模型基于原位監(jiān)測數(shù)據(jù)，通過冗余分析（Redundancy Analysis，RDA）研究地下水相關(guān)指標、多種環(huán)境因子與土壤含水率、土壤鹽度間的關(guān)系。由于解釋變量和環(huán)境變量的單位均不同，在進行RDA分析前，采用Z分數(shù)（Z-Score）標準化方法對所有變量進行標準化，其公式為：

式中：x為個案值；μ為總體平均值；σ為總體的標準偏差。

基于瞬時值的淺層土壤含水率（SM）和鹽度（SS）與環(huán)境因子的RDA 分析結(jié)果見圖4。前2 個排序軸特征值總和為98.14%，表明環(huán)境因子對淺層土壤含水率和鹽度均有較好的解釋。圖4顯示地表高程和潮高與淺層土壤鹽度顯著正相關(guān)，其他環(huán)境因子影響較小。其原因可能為高程越高，地表裸露時間越長，土壤中水分蒸發(fā)，而鹽分滯留在土壤中導(dǎo)致。然而，地下水鹽度與淺層土壤鹽度相關(guān)關(guān)系較弱，表明土壤鹽分變化受到地下水鹽度變化影響較小，土壤鹽分的變化是一個累積過程，與地下水鹽度瞬時值關(guān)系較弱。

采用廣義可加模型（Generalized Additive Models，GAMs）量化土壤鹽度與環(huán)境因子之間的關(guān)系。本研究中光滑函數(shù)使用薄板樣條函數(shù)，以約束性最大似然法作為光滑參數(shù)估計方法，模型擬合度檢驗使用殘差偏差法。

GAMs 模型的響應(yīng)變量為土壤含鹽度（SS），解釋變量為含水率（SM）、潮位（TIDE）、潮差（TIDEH）、地下水埋深（GWD）、地下水鹽度（GWS）、地表高程（ELE）、距離黃河河岸的距離（DIS）和黃河水位（YRWL）。采用Pearson相關(guān)系數(shù)法來檢查環(huán)境因子之間的共線性，發(fā)現(xiàn)ELE與DIS、SM與GWD、TIDE與TIDEH、DIS與GWS存在顯著相關(guān)關(guān)系（見圖5），在模型中只保留其中一個變量，因此本研究選取ELE、GWD和TIDE。

圖4 淺層土壤含水率和鹽度與環(huán)境因子RDA分析結(jié)果

圖5 環(huán)境變量的Pearson相關(guān)系數(shù)檢驗結(jié)果

通過單變量檢驗和交互因子檢驗設(shè)置光滑函數(shù)以及參數(shù)，GWD、GWS和ELE通過單變量顯著性檢驗，可以單獨作為模型的解釋變量，模型的解釋率在3.31%～50.90%。以GWD、GWS和ELE為模型參數(shù)基礎(chǔ)，考慮它們之間交互因子的影響，建立多個GAMs方程。3個因子之間的交互項均對解釋變量有顯著影響，模型解釋率在51.80%～68.10%，環(huán)境因子之間的交互作用對SS有顯著的影響。

以上述單、雙因子的GAMs模型為基礎(chǔ)，采用向前逐步回歸法對模型進行優(yōu)化。按照環(huán)境因子影響的顯著程度調(diào)整模型結(jié)構(gòu)，共建立10個含鹽度模擬GAMs模型（見表1）。以解釋率殘差偏差評價模型擬合度，解釋率越高殘差偏差越小，則模型擬合程度越好。為避免模型過度擬合，采用似然比檢驗判斷是否顯著提高模型性能，使用復(fù)雜模型和簡單模型的比例應(yīng)該近似卡方分布，故以P值進行判斷［18］。由表1可知，Model7的解釋率最高（69.50%）、殘差偏差最?。?78.50），是最優(yōu)模型，模型結(jié)構(gòu)為：

式中：s為單因素函數(shù)；te為交互作用函數(shù)。

最優(yōu)模型的SS計算值與監(jiān)測值擬合結(jié)果見圖6（a），決定系數(shù)R2=0.6959，表明模型預(yù)測值與監(jiān)測值間擬合度較高。利用2019年4、5和6月實測數(shù)據(jù)（n=12）進行模型驗證（見圖6（b）），預(yù)測值和實測值的R2為0.6505，模型可用于預(yù)測。

表1 土壤鹽度GAMs模型結(jié)構(gòu)優(yōu)化

圖6 最優(yōu)模型（Model7）的淺層土壤鹽度率定和驗證結(jié)果

3.3 植被種群增長-競爭動力學(xué)模型基于研究區(qū)典型鹽沼植被（翅堿蓬、檉柳和蘆葦）對地下水埋深和淺層土壤鹽度耐受性的實驗結(jié)果，將研究區(qū)剖分成1223×1038個30 m×30 m的正方形網(wǎng)格，耦合元胞自動機和改進的Logistic函數(shù)，建立研究區(qū)典型鹽沼植被種群增長-競爭動力學(xué)模型。

3.3.1 基于元胞自動機的種群擴散模塊 在自然條件下，鹽沼植被的種子的發(fā)芽率以及幼苗的定植率較低，地下莖的克隆繁殖是其主要的定植和繁殖擴散方式。鹽沼植被地下莖的克隆擴散過程受到植株狀態(tài)和環(huán)境因素的共同影響，其中土壤水、鹽指標是影響地下莖克隆擴散過程最為關(guān)鍵的因子。因此，本研究選用地下水埋深和淺層土壤鹽度作為根莖出芽率的函數(shù)。采用元胞自動機方法模擬［19］，傳遞規(guī)則為：（1）當(dāng)?shù)叵滤裆詈蜏\層土壤鹽度處在鹽沼植被可生長的范圍時，鹽沼植被繼續(xù)存活，并具有一定的生長力；（2）當(dāng)?shù)叵滤裆詈蜏\層土壤鹽度處在鹽沼植被最適宜生長的范圍時，鹽沼植被生長處于旺盛的狀態(tài)；（3）當(dāng)?shù)叵滤裆詈蜏\層土壤含鹽度處在鹽沼植被不可生長的范圍時，鹽沼植被在下一時刻死亡；（4）當(dāng)元胞沒有存活的鹽沼植被且地下水埋深和土壤含鹽度處于鹽沼植被可生長的范圍時，下一時刻該元胞內(nèi)有可能產(chǎn)生新的鹽沼植被，其概率P為：

式中：i為鄰域元胞編號；Ni為鄰域元胞的出芽數(shù)；r為元胞的出苗率。

在每一個元胞內(nèi)計算3種鹽沼植被的生境適宜度。假設(shè)該元胞內(nèi)的淹水深度和土壤鹽度兩個狀態(tài)中，存在任何一個狀態(tài)不適宜鹽沼植被生長，則鹽沼植被的生境不適宜，其計算公式為：

式中：SI為一種鹽沼植被的生境適宜度，其取值為0 或1，當(dāng)取值為1 時表明該元胞內(nèi)植被可以生長，取值為0時說明該元胞內(nèi)植被不可生長；SID為地下水埋深脅迫量，取1時表示該元胞內(nèi)該鹽沼植被不受地下水埋深脅迫，取0 時表示該元胞內(nèi)該鹽沼植被處于受地下水埋深脅迫無法正常生長；SIS為土壤鹽度脅迫量，取1時表示該元胞內(nèi)該鹽沼植被不受土壤鹽度脅迫，取0時表示該元胞內(nèi)該鹽沼植被受土壤鹽度脅迫而無法正常生長。

式中：Di,j為該元胞的土壤含水率；Si,j為該元胞的土壤鹽度值，g/L；Dmax和Dmin分別為地下水埋深的閾值；Smax和Smin分別為土壤鹽度的閾值。

3.3.2 基于改進Logistic函數(shù)的植被生長-競爭模塊 Logistic生長曲線適用于處于無環(huán)境干擾下的植物萌發(fā)后正常生長過程的模擬，其形狀為“S”形，分為發(fā)生、發(fā)展和成熟三個階段，斜率先增大再減小。由于Logistic函數(shù)忽略了鹽沼植被在不適宜環(huán)境下的生物量損失，因此本研究采用考慮生物量損失的Hill方程對Logistic函數(shù)進行修正［20］，同時引入鹽沼植被的種間競爭函數(shù)［10］，提出考慮植物生物量損失和種間競爭的Logistic函數(shù)的數(shù)學(xué)表達式：

式中：Δt為模擬的時間步長，本研究關(guān)注于生長季內(nèi)的變化，以月計；ΔRi為i月的鹽沼植被生物量凈增量，g/m2，本研究只考慮地上部分；ΔLi為i月的生物量增量，g/m2；ΔCi、ΔHi和ΔFi分別為i月因種間競爭、地下水埋深和土壤鹽度導(dǎo)致的生物量損失量（干重），g/m2；ri為i月的內(nèi)稟增長率；L為不適條件下的最大生物損失量，g/(m2·月)；Di為i月的地下水埋深，m；Dm為最適地下水埋深，m；Si為i月的土壤鹽度，g/L；Sm為最適土壤鹽度，g/L；p為模型形狀參數(shù)，取值為2、4 或6［20］；h為半致死地下水埋深，m；s為半致死土壤鹽度，g/L；βj,i-1為該鹽沼植被的種間競爭系數(shù)；Rj,i-1為與鹽沼植物競爭的另一種鹽沼植物j在i-1月的生物量，g/m2。

3.4 模型耦合耦合模型的結(jié)構(gòu)如圖7所示。基于地下水動力模型得到地下水水位分布條件，通過GAMs方程計算得到土壤鹽度分布條件，將土壤鹽度和地下水埋深分布作為種群動力學(xué)模型的輸入條件，利用元胞自動機模擬適宜生境的空間分布，在每個元胞內(nèi)利用改進的Logistic函數(shù)計算鹽沼植被的生物量。在此基礎(chǔ)上假設(shè)，生長季開始時，鹽沼植被在適宜生境處完成芽的擴散，在生長季內(nèi)不會產(chǎn)生新的植株，只在芽分布處進行生長。運行植被種群增長-競爭動力學(xué)模型計算鹽沼植被的適宜生境和生物量分布。

圖7 模型結(jié)構(gòu)圖

圖8 模型邊界及水文地質(zhì)參數(shù)分區(qū)

4 模型率定與驗證

4.1 地下水模型率定與驗證黃河三角洲自然保護區(qū)模擬區(qū)面積約為46 km×46 km，其中陸地面積約為687 km2，處于河口咸淡水交互區(qū)，黃河和渤海對模擬區(qū)的地下水埋深有顯著影響。將研究區(qū)劃分為10 000個460 m×460 m的矩形網(wǎng)格，地下水監(jiān)測井附近網(wǎng)格細化為230 m×230 m。根據(jù)全國重要地質(zhì)鉆孔數(shù)據(jù)庫服務(wù)平臺（http：//zkinfo.cgsi.cn/）提供的地質(zhì)鉆井圖與土壤巖性信息，將研究區(qū)土壤淺層含水層概化為非均質(zhì)、各向異性的潛水含水層（粉土黏土為主）和微承壓含水層（砂土為主）兩層，并通過空間插值的方法，將各含水層的頂板和底板概化，其中地表采用DEM衛(wèi)星數(shù)據(jù)。模型的東北側(cè)、東側(cè)和南側(cè)均為渤海，概化為已知水頭邊界；西北側(cè)為堤壩，概化為不透水邊界；西側(cè)為自然的大陸地塊，自西向東流量近似為零，概化為零流量邊界；地表存在黃河徑流，概化為河流邊界。降水補給和潛水蒸散發(fā)在研究區(qū)內(nèi)變化不大，滲透系數(shù)取值0.0167 m/d［21］。模型邊界及水文地質(zhì)參數(shù)分區(qū)見圖8，參數(shù)取值參考文獻［21］。

模型計算時段為2019年7月1—29日。使用4口地下水觀測井的實測數(shù)據(jù)進行模型參數(shù)率定和驗證，其中W1、W3和W7用于模型的參數(shù)率定，W5用于模型驗證。模擬結(jié)果與實測數(shù)據(jù)的擬合度采用標準化殘差（RMSE）進行評估，RMSE值越小擬合效果越好。模型調(diào)參后RMSE=7.79%，模型誤差在合理范圍內(nèi)，模型的地下水位時間序列圖（圖9（a））中模擬值與實測值較為接近，證明模型率定完成。以地下水觀測井W5實測數(shù)據(jù)進行模型的驗證。模型驗證中，RMSE=9.18%，模型的地下水位時間序列圖（圖9（b））中模擬值與實測值的較為接近，滿足精度要求。

圖9 模型率定和驗證結(jié)果

4.2 植被種群增長-競爭動力學(xué)模型率定與驗證在黃河三角洲地區(qū)關(guān)于蘆葦、翅堿蓬和檉柳對地下水埋深和土壤鹽度的耐受性，很多學(xué)者開展了室內(nèi)移植實驗、原位監(jiān)測實驗和模型模擬的研究。結(jié)合前人已有的研究成果和課題組的野外監(jiān)測，確定了各物種出芽數(shù)和出苗率的取值見表3?？紤]到目前相關(guān)研究中，對于地下水埋深脅迫和土壤鹽度脅迫影響的協(xié)同影響的研究較少，本研究不區(qū)分出芽數(shù)和出苗率在地下水埋深脅迫和土壤鹽度脅迫下的差異。

表2 3種典型鹽沼植被的生境閾值

表3 不同生境條件下各物種出芽數(shù)和出苗率取值

Logistic函數(shù)中，3種鹽沼植被的最大生物量的設(shè)置參照文獻［4］的取值，蘆葦為2437.24 g/m2，檉柳為3137.75 g/m2，翅堿蓬為5970.15 g/m2。模型的最適地下水埋深和最適鹽度取上下限的均值，蘆葦、檉柳和翅堿蓬的最適地下水埋深依次為-0.185、0.85和0.42 m，最適鹽度依次為6.5、16和12.71 g/L。3個物種的內(nèi)稟增長率取值參見表4。3種鹽沼植被在生長季開始時的初始生物量構(gòu)建的模型參數(shù)設(shè)為100 g/m2。

表4 內(nèi)稟增長率參數(shù)取值

表5 種間競爭系數(shù)取值

圖10 適宜生境模型模擬結(jié)果

3種鹽沼植物的種間競爭系數(shù)參考文獻［4］，在地下水埋深脅迫和土壤含鹽度脅迫下取值見表5。模型模擬了不同土壤鹽度下蘆葦、檉柳和翅堿蓬的生物量，與已有監(jiān)測結(jié)果的對比見圖10，表明模型合理。

5 結(jié)論

在徑流和潮汐的共同作用下，黃河口三角洲的水鹽條件發(fā)生動態(tài)變化，進而導(dǎo)致鹽沼植被群落呈明顯的帶狀分布。為明確生境條件與植被群落間的影響關(guān)系，探究濱海鹽沼濕地生態(tài)系統(tǒng)對自然和人類活動的響應(yīng)規(guī)律。本研究基于野外觀測數(shù)據(jù)建立了地下水水位、土壤水分和土壤鹽度的地下水動力模型；將元胞自動機模型和改進的Logistic函數(shù)相耦合，構(gòu)建了植被種群動力學(xué)生境模型，并將水動力模型與種群動力學(xué)生境模型相耦合，實現(xiàn)了包含土壤水鹽協(xié)同作用和植被生長擴散作用兩個過程的植被群落適宜生境及生物量的定量模擬和預(yù)測，突破了目前基于地下水動力學(xué)過程的適宜生境模型無法模擬植被生長擴散作用的局限。