膠州灣水體氮磷營養物基準研究?

徐夢琪，婁 琦，崔曉瑩，李正炎,2??

(1.中國海洋大學環境科學與工程學院, 山東 青島 266100; 2.中國海洋大學海洋環境與生態教育部重點實驗室, 山東 青島 266100)

海灣處于陸地和海洋相互作用強烈的區域，其生態系統較為復雜，極易受到人類活動和陸源污染的影響[1]。由于海灣地區人口分布密集，且快速的城市化和工業化進程使世界上大多數沿海水域都面臨著巨大的環境壓力[2-3]。過量的人為營養物，如氮(N)和磷(P)，被排放到沿海水域，導致其出現水質惡化和富營養化問題，從而引發有害藻華(Harmful algal blooms, HABs)等生態災害[4]。《海水水質標準》(GB 3097—1997)是為中國海域的環境監測和風險評估制定的，然而，考慮到不同地區的自然特征差異，目前制定的海水營養物標準還不能為所有海域提供充分有效的保護。為了保護各海灣生態系統免受富營養化和HABs的有害影響，應針對不同海灣等沿海水域制定特定的營養物基準，以開展相應的營養物控制工作。海灣營養物基準定義為對海灣地區沒有不利或有害影響的環境營養物參數的最大劑量或濃度[5]。營養物基準的制定和實施對海灣生態系統的保護、恢復和可持續發展具有重要意義。

在營養物基準研究方面，美國在過去幾十年中發布了一系列針對不同水體包括河流、湖泊和水庫、河口和沿海海域以及濕地[6-9]的技術指導文件。參照條件法、機械建模和壓力-響應模型是制定營養物基準的三種典型方法，由美國環境保護局(the United States Environmental Protection Agency, US EPA)推薦并廣泛運用[10-12]。除此之外，從水生態保護的角度，歐盟也制定了相應的技術指南使各成員國相對統一地建立營養閾值控制富營養化[13-15]。相比之下，在過去十幾年中，中國在營養物基準研究方面也取得一些重要進展。但已有研究主要集中在湖泊和河口區域[16-19]，針對海灣等沿海水域的營養物基準研究仍然十分匱乏。

膠州灣被稱為青島的“母親灣”，面臨來自工業、農業、旅游、港口等各方面的環境壓力。隨著環灣地區高強度的人類開發活動，近幾十年來N、P等營養物不斷輸入，并呈現不同程度的增加趨勢，同時，由于灣內相對較弱的水交換能力，造成了膠州灣水體污染和富營養化等問題[20-21]。近年來，膠州灣水體呈現明顯的“高N、低P、低硅(Si)”的營養鹽結構，總氮(Total nitrogen, TN)和總磷(Total phosphorus, TP)的年輸入通量分別為6 945和160.6 t/a[22]，且目前膠州灣N/P值呈現快速上升趨勢，Si/N值則一直遠低于Redfield比值[23]。由于這些變化，海灣的葉綠素a(Chlorophylla, Chla)濃度以及浮游植物的生物量都在增加，嚴重威脅著膠州灣的生態系統安全[22, 24]。目前，膠州灣的營養物基準尚未建立，這在一定程度上影響了膠州灣及其陸域營養物控制的科學性和富營養化風險評估。本文根據2014—2017年的膠州灣現場監測數據，采用頻數分布法、分類回歸樹模型以及非參數拐點分析等方法推導出營養物基準候選值，并基于赤潮敏感藻生態響應，通過室內多藻種實驗進一步開展營養物基準校驗，最終提出膠州灣水體的營養物基準推薦值。

1 研究區域和數據來源

膠州灣(35°97′N—36°18′N, 120°04′E—120°23′E，見圖1)是一個典型的半封閉型海灣，位于山東半島南部，平均水深7 m，最大水深64 m，出現在灣口附近，面積約360 km2。膠州灣三面被青島、紅島和黃島環繞，且有6條主要河流匯入膠州灣，為其帶來大量的點源和非點源營養物，顯著影響了膠州灣的生態環境[21]。自1990年以來，膠州灣內赤潮頻發，記錄次數超過30次，這對青島市社會經濟的可持續發展造成了威脅[25-26]。

本研究中使用的膠州灣環境監測數據來自中國海洋大學未公開的現場調查結果和國家海洋環境監測中心的水質監測報告。現場調查在2014—2017年的3—10月進行，采樣站點共22個(見圖1)，監測指標包括鹽度(S)、pH、溶解無機氮(Dissolved inorganic nitrogen, DIN)、活性磷酸鹽(Soluble reactive phosphorus, SRP)、TN、TP、氨氮(NH3-N)、溶解氧(Dissolved oxygen, DO)和Chla。所有監測指標的樣品采集、保存、運輸、預處理方法、化學分析和記錄規范均按照《海洋監測規范》(GB17378.4—2007)等相關標準進行。

圖1 膠州灣研究區域及采樣點位置

2 研究方法

2.1 營養物基準指標篩選

一般來說，營養狀態變量涉及營養物指標和生物指標。營養物指標，如DIN、SRP、TN和TP，屬于壓力變量，能夠表征水體的富營養化狀況；藻類對營養物質的變化很敏感并能做出直接的反應，藻類生物量(例如Chla)和DO等生物指標則屬于響應變量，反映了營養物與生物初級生產力的關系[16]。

本研究結合膠州灣2014—2017年3—10月的調查監測數據，運用相關性分析(Correlation analysis, CA)和主成分分析(Principal component analysis, PCA)，綜合考慮我國海水水質評價的主要監測指標類別，進行膠州灣營養物基準指標的篩選。其中，CA是分析連續變量之間線性相關程度，根據變量數據的正態性檢驗，計算皮爾森相關系數(Pearson correlation coefficient)評估相關性。運用CA確定變量之間的相關性后進行PCA。PCA在盡量保證原信息量不丟失的基礎上，按照方差最大化的原則，通過某種線性投影使所投影的維度上數據信息量最大，轉換后的變量為主成分。CA和PCA均運用SPSS 26軟件進行數據處理和統計分析。

2.2 營養物基準推導方法

2.2.1 頻數分布法 頻數分布法是確定水體參照狀態時常用的統計分析方法，并在國內外研究中廣泛應用于制定特定營養狀態的營養物基準。頻數分布圖用于分析DIN、SRP、TN、TP、Chla和DO監測指標數據的正態性，如果數據不呈現正態分布，則采用ln進行數據轉換。分析營養狀態變量的頻數分布曲線時，如果沿海水域受到富營養化的損害程度較小，通常提出將所有采樣點數據頻數分布曲線的下25%或參照點數據頻數分布曲線的上25%作為營養物基準值。如果大多數采樣點的水體環境因富營養化而嚴重退化，則應選擇將監測數據頻數分布曲線的上5%或下5%作為營養物基準值[6, 8]。頻數分布法的分析及作圖均使用SPSS 26軟件進行處理。

2.2.2 壓力響應模型 當無法確定受富營養化影響最小的地點及其歷史數據時，可以根據現場調查的監測指標數據來分析養分濃度和初級生產力之間的壓力響應關系，通過建立壓力-響應模型來推導營養物基準值[10]。該模型能夠定量描述營養物壓力變量與生物響應變量之間的關系，尤其適用于受人為活動影響明顯水域的營養物基準。為了探索壓力變量和響應變量之間復雜的非線性關系，分類回歸樹模型(Classified and regression tree, CART)和非參數拐點分析法(Nonparametric changepoint analysis, nCPA)是兩種常見的的壓力-響應模型方法[19, 27]，本研究基于膠州灣2014—2017年3—10月的調查監測數據，應用CART模型和nCPA兩種方法分析并進行膠州灣營養物基準候選值的推導。

CART模型可以揭示一系列響應和壓力變量之間可能存在的非線性關系[28]。CART模型通過基于基尼系數的二叉樹結構對壓力變量進行分類，壓力變量的每一次分類都是根據其最佳分類值來完成，分組后按照響應變量的最大降低目標變量的分布或變異性將其分成兩個獨立的節點。重復這個分類過程，直到不能進行節點分類并形成終節點為止[30]。通常，分類回歸樹的第一個根節點是最重要的壓力指標，對響應變量有顯著影響。

nCPA法是一種估計二元關系中變化點位置的方法，為營養物基準提供了自然候選值。從操作上來講，將響應變量沿著壓力梯度排序，并將拐點定義為兩組響應變量之間差異達到最大值的點[10, 30]。此外，使用隨機自舉抽樣和90%的置信區間來定義拐點位置的不確定性。CART模型和nCPA兩種方法能夠消除設置基準值的主觀偏差，一般來說，CART模型和nCPA這兩種方法得到的基準值可以相互驗證，以提高確定營養物基準的準確性[19]。CART模型和nCPA法均使用R語言(R64 3.5.1)軟件進行數據處理及作圖。

2.3 營養物基準值驗證方法

不同于單種浮游藻的驗證實驗，本研究考慮了浮游藻的種間競爭，更好地模擬膠州灣的水生態環境，從種群的角度推導生態響應值進行營養物基準值的驗證。根據營養物水平與藻類生長之間的壓力響應關系，通過混合藻種實驗驗證上述營養物基準候選值。驗證過程包括實驗數據的獲取、壓力響應模型的建立和生態響應值的計算。

2.3.1 混合藻種實驗 近年來，中肋骨條藻(Skeletonemacostatum)和旋鏈角毛藻(Chaetoceroscurvisetus)一直是膠州灣的常見優勢藻種，同樣也是主要的赤潮藻種。膠州灣最常暴發雙相型赤潮，S.costatum也最易與其他藻種同發赤潮[31-32]。S.costatum、C.curvisetus是常見的廣溫廣鹽藻種，在夏季易引發赤潮，對于這些廣泛分布的赤潮藻種，應該重視其形成赤潮的潛在風險。除此之外，有研究表明[22-33]，引發膠州灣赤潮的浮游藻中，甲藻的比例明顯增高，并且值得注意的是，塔瑪亞歷山大藻(Alexandriumtamarense)曾于夏季在膠州灣出現過[34]。由于A.tamarense為有毒赤潮藻，需要對其加強監測。因此，結合室內可培養情況，混合藻種實驗選擇S.costatum、C.curvisetus和A.tamarense為研究對象。

混合藻種實驗進行前，將指數生長期的S.costatum、C.curvisetus和A.tamarense分別于1 700、2 000和1 000 r·min-1轉速下離心10、10和5 min，去除上清液并保留下層藻液，然后用鹽度為30±1的滅菌人工海水進行沖洗再離心2次，從而減少藻液中殘存的營養鹽。根據膠州灣水體中DIN和SRP的環境濃度水平，加入NaNO3和NaH2PO4以設置8個不同N、P濃度的實驗組(見表1)，其它營養物按照f/2培養基的配方加入滅菌后的人工海水。將離心后的三種藻液以初始密度約為1×104cells/mL置于培養基中共培養。實驗分為兩組，即固定P濃度的N組和固定N濃度的P組，并為每個濃度水平設置三個平行樣。為了使藻在生長過程中有充足的N(P)源，固定N組中N濃度及固定P組中P濃度高于實際監測濃度。不同N(P)濃度梯度則主要依據膠州灣現場監測數據中DIN和SRP的濃度范圍，同時考慮到曲線擬合的相關要求進行設置。將錐形瓶置于恒溫光照培養箱(GXZ-280B，寧波江南儀器廠)中培養，溫度為20 ℃，光照強度為3 000～5 000 Lx，光/暗循環時間為12 h∶12 h。為防止藻類聚集，錐形瓶每天定時搖動三次。并于每天的同一時間，取1 mL藻液于取樣管中，在顯微鏡下用血細胞計數板進行中肋骨條藻和旋鏈角毛藻的細胞密度計數，用浮游動物計數框進行塔瑪亞歷山大藻的計數。

表1 實驗組中的N和P濃度

2.3.2 數據分析 藻類生長曲線通常呈S形，并可以用各種模型進行擬合，擬合模型包括Logistic、Bertalanffy、Mitscherlich和Brody[35]。在這些模型中，S型生長曲線被最廣泛地用來描述生物生長過程。因此，本研究采用Logistic模型來擬合三種藻的指數生長階段。S-logistic2(1)和S-logistic3(2)分別適用于固定P組和固定N組的生長曲線擬合，具體公式如下：

(1)

(2)

式中：Bt是t時的浮游植物生物量(cells/mL)；B0是初始浮游植物生物量(cells/mL)；Bf是終止生物量(cells/mL)；代表浮游植物種群的最大密度；R是浮游植物的種群增長率(cells·mL-1·d-1)；t是培養時間(d)；b和k是模型參數。

隨著營養物濃度的變化，Bf的形狀大致與S型曲線相同，并且可以用Boltzmann模型來擬合。曲線的拐點x0是Bf的突變點，可用于確定生態響應值。Boltzmann方程(3)表示如下：

(3)

式中：a,b和c為模型參數；x為營養物濃度(mg/L)。

在此混合藻種實驗中，固定N(P)實驗組設置為8個濃度梯度，每個濃度水平設有三個平行樣。培養7天后，將每個濃度梯度的每個平行樣的藻密度值作為一組。根據各組的實驗測量數據(t,Bt)，通過Logistic擬合生長曲線得出Bf值。然后，進行Bootstrap隨機抽樣，從每個濃度梯度中隨機抽取Bf值，8個濃度梯度提取的Bf值形成一組數據，總計38組數據。此外，從38組數據中隨機抽取200組(x,Bf)數據，并利用Boltzmann模型進行擬合，得到200個x0值。將200個x0值的中值作為生態響應值，并使用95%的置信區間來表示平均響應值的不確定性。在獲得各藻種的響應值后，選擇最低值作為混合藻種實驗的生態響應值，以代表最敏感藻種的響應。使用Matlab 2019b和OriginPro2021軟件進行Bootstrap隨機抽樣和模型擬合。

3 結果

3.1 營養物基準指標的選擇

從各監測指標的相關性系數矩陣和顯著性檢驗結果(見表2)中可以直觀地看到，各監測指標之間是具有明顯相關性的。其次，表3是PCA結果，包括各成分的特征值、成分貢獻率和累積貢獻率。成分1、2、3的特征值分別為2.997、1.748、1.124，均大于1且其累積貢獻率達到73.361%，因此提取前三個成分為主成分。此外，通過各因子旋轉形成旋轉后的成分矩陣(見表4)來分析與各主成分密切相關的監測指標信息。

表2 各指標的相關系數矩陣與顯著性檢驗結果

表3 特征值、成分貢獻率與累積貢獻率

表4 旋轉后的成分矩陣

PCA結果顯示，DIN、SRP、TN和TP為第一主成分相關的指標，第二主成分相關指標為Chla和S，pH和DO則為第三主成分相關指標。DIN和SRP通常用于評估富營養化狀況，考慮我國《海水水質標準》(GB 3097—1997)中的主要監測指標類別，將DIN和SRP選為主要的營養物指標，TN、TP為輔助性營養物指標。除此之外，由于S與DIN、SRP、TN和TP的相關性較弱，且pH容易受藻類生長、海水擾動等作用影響而產生較大波動，因此不考慮將S和pH作為營養物基準的指標。而CA結果表明DO與DIN和SRP密切相關，且DO是生態健康的重要指標，同時也是水體營養狀況的早期預警指標。除此之外，Chla含量表明浮游植物的生物量，是HABs最具代表性的響應指標[36-37]。隨著營養物濃度的升高，浮游植物密度增加，DO含量減少。因此，DIN、SRP、TN和TP被選為壓力變量，Chla和DO被選為響應變量。

3.2 頻數分布法

我國《海水水質標準》(GB3097—1997)將水質分為四類，通過對膠州灣的監測數據進行統計分析，發現76.7%采樣點的DIN和98.3%采樣點的SRP屬于前三類水質，同時考慮到膠州灣整體的富營養化狀況，本研究基于膠州灣2014—2017年3—10月的調查監測數據，選擇將所有采樣點數據的頻數分布曲線的下25%來推導膠州灣營養物基準的候選值。頻數分布曲線及其統計值如圖2和表5所示，DIN、TN、SRP和TP的營養物基準候選值分別為0.152、0.394、0.003 5、0.008 9 mg/L(見圖2和表5)。

圖2 膠州灣營養狀態指標的頻數分布曲線

表5 膠州灣營養狀態指標的頻數分布統計結果

3.3 分類回歸樹模型

以Chla和DO為響應變量，DIN、SRP、TN和TP為壓力變量，使用CART模型建立壓力響應關系，選擇CART根節點的養分濃度作為膠州灣營養物基準候選值，并將DIN和TN的結果保留三位有效數字，以Chla和DO為響應變量的CART分析結果如圖3、圖4和圖5所示。從圖3中可以看出，SRP是影響Chla和DO濃度變化的最首要且最重要的壓力變量，這與目前膠州灣存在的P限制現象相符合。以Chla為響應變量，得出DIN、TN、SRP和TP的營養物基準候選值分別為0.257、0.410、0.003 2和0.009 8 mg/L(見圖4)；將DO作為響應變量，計算出DIN、TN、SRP和TP的營養物基準候選值分別為0.235、0.885、0.010 9和0.035 5 mg/L(見圖5)。

(箱圖是根據分類數據中Chl a和DO的平均濃度升序排列，表示每個終端節點上Chl a和DO的濃度分布。The boxplots represent the concentration distributions of Chl a and DO at each terminal nodes, arranged according to the average Chl a and DO concentrations in the classified data in an ascending order.)

(箱圖是根據分類數據中Chl a的平均濃度升序排列，表示每個終端節點上Chl a的濃度分布。The boxplots represent the concentration distributions of Chl a at each terminal nodes, arranged according to the average Chl a concentrations in the classified data in an ascending order.)

(箱圖是根據分類數據中DO的平均濃度升序排列，表示每個終端節點上DO的濃度分布。The boxplots represent the concentration distributions of DO at each terminal nodes, arranged according to the average DO concentrations in the classified data in an ascending order.)

3.4 非參數拐點分析法

以DO為響應變量，DIN、SRP、TN和TP為壓力變量，使用nCPA方法建立壓力響應關系。將DIN、SRP、TN和TP的拐點作為這些變量對DO的突變點，并將DIN和TN的結果保留三位有效數字。nCPA的結果顯示(見圖6)，DIN、TN、SRP和TP的營養物基準候選值分別為0.235、0.885、0.010 9和0.035 5 mg/L。由于Chla數據不足，未對該響應變量進行nCPA分析。

(折線表示響應變量在拐點處的階梯變化，虛線表示拐點90%的置信區間，圓圈表示樣本數據。The broken line shows the stepped variation of the response variable at the changepoint, the vertical dashed lines show the 90% confidence intervals of the changepoint, and the circles represent sample data.)

3.5 混合藻種實驗

C.curvisetus和A.tamarense的生長曲線與總生物量有很大不同，但S.costatum的生長曲線與總生物量相似(見圖7、8)。在固定N和P組中，S.costatum為優勢種，占總生物量的大部分，且由于營養物質豐富，S.costatum和總生物量在0～3 d和0～4 d階段的生長曲線基本呈S形，因此，選擇固N組0～3 d和固P組0～4 d的階段進行Bf值的擬合。S.costatum的DIN和SRP生態響應值分別為0.521和0.044 6 mg/L，總生物量的DIN和SRP生態響應值分別為0.642和0.047 4 mg/L(見表6)。因此，對DIN和SRP最敏感藻種的生態響應值分別為0.521和0.044 6 mg/L。根據膠州灣的監測數據，SRP∶TP和DIN∶TN比值的平均值分別為0.315和0.519，TN和TP的生態響應值則分別為1.004和0.142 mg/L。

圖7 固N組中肋骨條藻、旋鏈角毛藻、塔瑪亞歷山大藻和總生物量的生長曲線

圖8 固P組中肋骨條藻、旋鏈角毛藻、塔瑪亞歷山大藻和總生物量的生長曲線

表6 DIN和SRP的生態響應值

3.6 營養物基準推薦值

通過頻數分布法、分類回歸樹模型和非參數拐點分析法得出的營養物基準候選值和通過混合藻種實驗計算出的生態響應值如表7所示。鑒于對膠州灣現場監測數據的統計分析，76.7%采樣點的DIN和98.3%采樣點的SRP屬于前三類水質，綜合考慮膠州灣目前的水質狀況和富營養化水平，選擇所有方法得出的基準最低值作為膠州灣的營養物基準推薦值，以此對膠州灣的營養物控制和評價提出更高的標準和要求。因此，本研究建議膠州灣的DIN、TN、SRP和TP營養物基準值分別為0.152、0.394、0.003 2和0.008 9 mg/L。

表7 不同方法推導的營養物基準候選值

4 討論

由于野外調查和室內實驗的環境條件不同，營養物基準候選值與生態響應值相差較大(見表7)。沿海地區N、P的濃度會受陸源輸入、S、水停留時間、海水交換、沉積特征、海洋環流和大氣沉降等的影響[38-39]。因此，基于現場監測數據的營養物基準候選值也會受到上述因素的影響。基于營養物壓力指標和生物響應指標之間的關系，混合藻種驗證實驗模擬了浮游藻最佳生長條件，為藻類生長提供了固定的溫度、S和光照條件。然而，實驗室無法完全模擬出海灣自然環境中的水文狀況和生物環境，使得營養物質的濃度和組成結構是影響生態響應值的主要因素。膠州灣目前呈現高N、低P的營養鹽結構，因此，本研究提出的DIN和SRP營養物基準推薦值的N/P值相應較高。為了避免過高的N/P值可能會導致生態系統失衡的結果，在后續的研究中，需要進一步進行現場微宇宙實驗，綜合考慮水文條件、浮游動物等膠州灣自然生態環境條件的影響，來驗證理論推導方法得出的營養物基準候選值，同時結合較為適合的N/P值進行營養物基準候選值的校正，在實現有效控制膠州灣N、P濃度的前提下，進一步保證膠州灣生態系統的平衡和穩定。

從表8可以看出，DIN和SRP的營養物基準值在不同的研究中差異較大，這可能是由于各個國家和地區沿海海域富營養化狀況的差異以及對營養物質的敏感性不同造成的。此外，推導方法選擇的不同，也是導致營養物基準值存在一定差異的原因。因早期遼河口水質狀況較為良好，除遼河口沿海地區外，膠州灣的DIN基準值低于其他沿海地區。而不同沿海地區的SRP基準值大約是膠州灣的2～9倍，這可能是與膠州灣在2010年以來存在的P限制現象有關[23]。本研究主要針對DIN和SRP指標，而TN和TP僅作為輔助性指標。在現行的《海水水質標準》(GB3097—1997)中，水質監測和評價的營養物指標也主要針對DIN和SRP，而TN和TP主要作為淡水中的營養物基準指標。不同研究中TN和TP的營養物基準值如表9所示，由于不同國家和地區藻類對營養物質的敏感性差異，除長江口外及舟山鄰海、佛羅里達群島外，膠州灣推薦的TN基準值低于表9中其它地區中營養物基準的范圍。而膠州灣推薦的TP基準值低于或在其他研究的范圍內，但比斯坎灣除外。

表8 不同沿海地區DIN和SRP的參照條件或營養物基準

表9 不同國家和地區中TN和TP的營養物基準

此外，有研究表明膠州灣在20世紀70年代以前，受人為影響很小，N、P等營養物基本呈現背景值且較為穩定，其中，SRP濃度在20世紀60年代約為0.005 4 mg/L。自20世紀80年代起，膠州灣DIN和SRP濃度呈現波動變化趨勢，目前膠州灣DIN濃度略高于國家Ⅰ類海水水質標準，而SRP濃度基本符合國家Ⅰ類海水水質標準[22-23]。本研究提出的DIN和SRP營養物基準推薦值均低于我國《海水水質標準》(GB 3097—1997)的Ⅰ類。TN的營養物基準推薦值在我國《地表水環境質量標準》(GB 3838—2002)的Ⅰ類和Ⅱ類范圍內，TP的營養物基準推薦值低于地表水水質的Ⅰ類。總體而言，本研究中建議的膠州灣DIN、TN、SRP和TP營養物基準值是相對合理可靠的，能對膠州灣實現有效保護。

5 結論

(1)根據2014—2017年的膠州灣調查監測數據，運用相關性分析和主成分分析選擇DIN、SRP、TN和TP作為營養物指標即壓力變量，Chla和DO作為生物指標即響應變量。

(2)運用頻數分布法、分類回歸樹模型以及非參數拐點分析三種方法推導出DIN、TN、SRP和TP營養物基準候選值，然后通過混合藻種實驗計算出生態響應值開展進一步驗證，最終提出膠州灣DIN、TN、SRP和TP的營養物基準推薦值分別為0.152、0.394、0.003 2和0.008 9 mg/L。

(3)通過比較不同國家和地區的營養物基準研究，本研究中建議的DIN、TN、SRP和TP營養物基準值合理可靠，能夠為膠州灣的水質監測及評估等環境管理工作提供科學支撐。