馮永亮

(唐山學(xué)院 基礎(chǔ)教學(xué)部，河北 唐山 063000)

0 引言

物種敏感度分布(Species Sensitivity Distribution，SSD)是一種用于描述不同物種對某個脅迫因素敏感度差異的統(tǒng)計分布模型。SSD的概念最早于1978年由美國環(huán)保局(U.S.EPA)在推導(dǎo)全國環(huán)境水質(zhì)標(biāo)準(zhǔn)時提出，隨后于1985年確定了以SSD的5%分位數(shù)所對應(yīng)的濃度(HC5)推導(dǎo)環(huán)境基準(zhǔn)值，并以此代替之前的專家判斷法，沿用至今[1]。我國在2017年頒布的《淡水水生生物水質(zhì)基準(zhǔn)制定技術(shù)指南》中也推薦使用SSD模型進(jìn)行淡水水生生物水質(zhì)基準(zhǔn)推導(dǎo)。此外，隨著SSD理論的發(fā)展，許多學(xué)者開始采用SSD與環(huán)境監(jiān)測數(shù)據(jù)相結(jié)合的概率方法評價污染物的生態(tài)風(fēng)險[2-5]。Chen[3]分別采用急性和慢性毒性數(shù)據(jù)，構(gòu)建了9種污染物的SSDs曲線，并結(jié)合環(huán)境監(jiān)測濃度構(gòu)建聯(lián)合概率曲線(JPCs)評價了中國臺灣北部基隆河中污染物的生態(tài)風(fēng)險，發(fā)現(xiàn)Zn和Cu的風(fēng)險最高，有機(jī)物的風(fēng)險可忽略。目前SSD模型在許多國家和地區(qū)得到了很好的應(yīng)用[6-7]，已經(jīng)成為污染物環(huán)境質(zhì)量標(biāo)準(zhǔn)制定和生態(tài)風(fēng)險評價的核心方法。

隨著SSD理論的發(fā)展和廣泛應(yīng)用，一些學(xué)者或機(jī)構(gòu)為SSD模型的構(gòu)建和生態(tài)風(fēng)險評價開發(fā)了相應(yīng)的軟件。這些軟件大體上可分為兩類：一類是基于Excel表格開發(fā)的，如USEPA開發(fā)的一款用于因果評估的在線應(yīng)用程序(CADDIS)中的SSD Generator，荷蘭國家公共衛(wèi)生和環(huán)境研究所(RIVM)開發(fā)的ETX 2.0[8]；另一類是基于專業(yè)的統(tǒng)計軟件R開發(fā)的，如澳大利亞聯(lián)邦科學(xué)和工業(yè)研究組織(CSIRO)基于Burr III分布開發(fā)的軟件BurrlizO[9]，Kon等開發(fā)的網(wǎng)頁版工具M(jìn)OSAIC_SSD[10]。這些軟件主要集中于SSD的構(gòu)建和環(huán)境標(biāo)準(zhǔn)(HC5)的推導(dǎo)，存在著可選統(tǒng)計模型較少、未能與環(huán)境數(shù)據(jù)相結(jié)合等缺點。本研究結(jié)合環(huán)境監(jiān)測數(shù)據(jù)，首次給出了基于SSD的概率生態(tài)風(fēng)險評價的Matlab函數(shù)“PERA”，為污染物的生態(tài)風(fēng)險評價提供技術(shù)支持。

1 SSD的原理和假設(shè)

不同生物生活史、生理、形態(tài)、行為和地域分布的差異構(gòu)成了生態(tài)系統(tǒng)的生物多樣性。對于生態(tài)毒理學(xué)，這些生物多樣性意味著不同生物對特定污染物具有不同的敏感性。物種間的這種敏感性差異可以利用特定的概率分布曲線即SSD曲線來描述。

物種敏感度分布的基本假設(shè)是不同物種對特定污染物的敏感性能夠用一些統(tǒng)計分布模型(如Log-normal和Log-logistic分布等)來描述。生態(tài)毒理學(xué)的數(shù)據(jù)可以看成是來自于這些分布模型的樣本，并用于SSD參數(shù)的估計。SSD模型的基本形式見圖1。

圖1 物種敏感度分布曲線示意圖

圖1中的圓圈為原始毒性數(shù)據(jù)點，可以是急性(如半數(shù)致死濃度，LC50)或慢性(如無觀測效應(yīng)濃度，NOECs)效應(yīng)濃度。擬合的曲線為毒性數(shù)據(jù)的累積概率曲線(即SSD)。圖中的箭頭表示SSD具有正向和反向兩個方面的應(yīng)用。其中反向(箭頭Y→X)可應(yīng)用于污染物的環(huán)境質(zhì)量標(biāo)準(zhǔn)(EQC)的獲得，主要是通過曲線計算特定物種損害比例所對應(yīng)的污染物濃度，通常取HC5作為EQC，從而保護(hù)生態(tài)系統(tǒng)中95%的生物不受污染物的影響；正向(箭頭X→Y)可以用于生態(tài)風(fēng)險評價，主要由SSD曲線計算污染物在特定濃度下(一般為環(huán)境監(jiān)測濃度)導(dǎo)致有害生物效應(yīng)的累積概率，即受影響物種的比例(PAF)，進(jìn)而定量表征污染物的生態(tài)風(fēng)險。

2 基于SSD的概率生態(tài)風(fēng)險評價

化學(xué)物質(zhì)的生態(tài)風(fēng)險主要是通過環(huán)境監(jiān)測數(shù)據(jù)與毒性數(shù)據(jù)進(jìn)行表征的，表征的方法包括簡單的風(fēng)險商HQs法和相對復(fù)雜的概率生態(tài)風(fēng)險評價法。風(fēng)險商HQs法(即環(huán)境監(jiān)測濃度除以毒性閾值)是簡單的單點估計法，具有簡單、透明和對數(shù)據(jù)要求低等優(yōu)點，但其只適用于污染物篩選的初級階段[11-12]。為了更加精確詳細(xì)地表征生態(tài)風(fēng)險，很多學(xué)者推薦采用更為復(fù)雜的概率方法，即聯(lián)合概率曲線(JPC)法。JPC方法通過比較環(huán)境監(jiān)測數(shù)據(jù)和毒性數(shù)據(jù)的分布特征可以定性和定量地表征生態(tài)風(fēng)險，常被用作細(xì)化各種污染物生態(tài)風(fēng)險的高水平方法[13]。下面給出JPC的構(gòu)建過程和其生態(tài)學(xué)意義。

設(shè)EC為污染物的環(huán)境監(jiān)測濃度，SS為污染物的毒性效應(yīng)濃度，隨機(jī)變量X1和X2分別來自于EC和SS的概率分布。

Pr(logEC>logSS)=Pr(X1>X2)。

(1)

式(1)為污染物的環(huán)境監(jiān)測濃度大于其毒性效應(yīng)濃度的概率，有兩種等價的積分形式表達(dá)式：

(2)

或

(3)

其中PDFXi，CDFXi(i=1，2)分別表示隨機(jī)變量Xi的概率密度函數(shù)和累積分布函數(shù)；式(3)中1-CDFX1(x)表示隨機(jī)變量X1超過給定數(shù)值x的概率，即傳統(tǒng)概率論中的生存函數(shù)，在環(huán)境毒理學(xué)中通常將其定義為超越函數(shù)EXFX1(x)=1-CDFX1(x)。此時式(3)可等價轉(zhuǎn)化為下式

(4)

由式(2)和式(4)右邊的積分表達(dá)式可知，Pr(X1>X2)即為污染物環(huán)境風(fēng)險在統(tǒng)計學(xué)上的期望值，其在數(shù)值上等于Van Straalen[14]給出的生態(tài)風(fēng)險δ和Cardwell[15]給出的期望總風(fēng)險(ETR)。

基于式(4)，以毒性數(shù)據(jù)的累積概率(即dCDFX2(x))為自變量，環(huán)境監(jiān)測數(shù)據(jù)的超越概率(即EXFX1(x))為因變量，構(gòu)造的曲線即為聯(lián)合概率曲線(具體構(gòu)建過程見圖2)。該曲線上的點表示導(dǎo)致不同物種損害水平的概率[2]，曲線與兩坐標(biāo)軸圍成的面積即為污染物環(huán)境風(fēng)險的期望值Pr(X1>X2)，也有學(xué)者稱該值為發(fā)生有害生物效應(yīng)的總體風(fēng)險概率(ORP)[16-17]。

類似地，JPC也可以根據(jù)式(2)進(jìn)行構(gòu)建，即暴露數(shù)據(jù)的累積概率(dCDFX1(x))為自變量，毒性數(shù)據(jù)的累積概率(CDFX2(x))為因變量，曲線與兩坐標(biāo)軸圍成的面積為Pr(X1>X2)。

3 概率生態(tài)風(fēng)險評價的Matlab程序

基于上述分析過程，構(gòu)建用于污染物概率生態(tài)風(fēng)險評價的Matlab函數(shù)“PERA”，其程序代碼如下所示。

function [HCq，ORP]=PERA(Tdata，q，option，ECdata)

% Tdata為污染物的毒性數(shù)據(jù)，需為向量形式。

% q為SSD的q分位數(shù)HCq，0

% option=1和option=2分別表示基于Log-normal和Log-logistic分布構(gòu)造SSD曲線。

% ECdata為污染物的環(huán)境監(jiān)測數(shù)據(jù)，需為向量形式。

% [HCq]=PERA(Tdata，q，option)，返回SSD曲線和HCq值。

% [HCq，ORP]=PERA(Tdata，q，option，ECdata)，除返回SSD曲線和HCq值外，還返回JPC曲線和ORP值。

if nargin<3

error('輸入?yún)?shù)必須為3個或4個')；

end

data=sort(Tdata)；

[n，m]=size(data)；

n=max(n，m)；

p=zeros(n，1)；

for i=1：n

p(i)=i/(n+1)；

end

data1=log10(data)；

if option==1

dis='normal'；

else if option==2

dis='logistic'；

end

end

[log，logci]=mle(data1，'distribution'，dis)；

p2=icdf(dis，q，log(1)，log(2))；

HCq=10^p2；

N=10000；

mm=0.2；

mi=min(data1)；

mi=min(mi)-mm；

mx=max(data1)；

mx=max(mx)+mm；

x1=linspace(mi-0.1，mx+0.1，N)；

p11=cdf(dis，x1，log(1)，log(2))；

p=p*100；

p11=p11*100；

MS=6；

figure

plot(data1，p，'ro'，'MarkerSize'，MS，'MarkerFaceColor'，'r')

hold on

LW=1.5； % 圖像中線的寬度。

plot(x1，p11，'k-'，'LineWidth'，LW)

axis([mi-0.1，mx+0.1，0，101])；

FS=12；% 圖像中字號的大小。

xlabel('log_{10}(concentration)')

ylabel('Potentially Affected Percentage (%)')

set(gca，'FontSize'，F(xiàn)S，'Fontname'，'times new Roman')

set(get(gca，'XLabel')，'FontSize'，F(xiàn)S，'Fontname'，'times new Roman')

set(get(gca，'YLabel')，'FontSize'，F(xiàn)S，'Fontname'，'times new Roman')

%%以下程序用于構(gòu)造聯(lián)合概率曲線JPC。

if nargin==4

data2=log10(ECdata)；

[phat，pci]=mle(data2，'distribution'，'normal')； %擬合環(huán)境數(shù)據(jù)的參數(shù)。

y=0：0.00001：1；

p22=icdf(dis，y，log(1)，log(2))； %求SSD的反累積對數(shù)濃度數(shù)值。

p33=cdf('normal'，p22，phat(1)，phat(2))； %求相對于SSD反累積對數(shù)濃度數(shù)值對應(yīng)的環(huán)境監(jiān)測數(shù)據(jù)的累積函數(shù)值。

p44=1-p33； %求相對于SSD反累積數(shù)值對應(yīng)的環(huán)境監(jiān)測數(shù)據(jù)的超越概率。

ORP=trapz(y，p44)； %計算JPC曲線與兩坐標(biāo)軸所圍區(qū)域的面積。

p44=p44*100；

y=y*100；

figure

plot(y，p44，'r'，'LineWidth'，LW)

title('JPC'，'Fontname'，'times new Roman'，'FontSize'，F(xiàn)S)

xlabel('Percent of species affected (%)')；

ylabel('Probability exceedance (%)')；

set(gca，'FontSize'，F(xiàn)S，'Fontname'，'times new Roman')

set(get(gca，'XLabel')，'FontSize'，F(xiàn)S，'Fontname'，'times new Roman')

set(get(gca，'YLabel')，'FontSize'，F(xiàn)S，'Fontname'，'times new Roman')

end

4 Matlab程序的應(yīng)用實例

采用Li[5]對萊州灣西部表層海水重金屬Zn的環(huán)境監(jiān)測數(shù)據(jù)(2016年9月監(jiān)測)和Zn對海洋生物的慢性毒性數(shù)據(jù)對本文Matlab程序進(jìn)行實際運行。

給定Zn的毒性數(shù)據(jù)Tdata=[20，1000，1144.962，97，1000，2154.395，10000，142.5，313，100，10，5000，400，380，4500，1000，970，200，139，480，230，600，520，124.0718，200，178，2500，630，88.5，15000，170.2，605，55，160，64，10，500，616.4414，80，170.9976，100，270，674，2000，2768.212，960，421.7163，1000，250，160，5900，1944.402]；

Zn的環(huán)境監(jiān)測數(shù)據(jù)ECdata=[39，40.6，39.1，38.6，43.6，40.4，41.2，41.9，41.1，37.8，38.4，36.8，38.4，38.7，39.1，38.5，41.2，40.7，42.4，39.5]；

在Matlab 2017b下運行：

>> [HC5，ORP]=PERA(Tdata，0.05，2，ECdata)

HC5=29.0852

ORP=0.0698

所得SSD曲線和JPC曲線分別見圖3和圖4。該結(jié)果與Li[5]的研究成果一致。

圖4 萊州灣西部表層海水中Zn的聯(lián)合概率曲線(JPC)

5 討論

專業(yè)的統(tǒng)計軟件R是學(xué)者們在做關(guān)于污染物生態(tài)風(fēng)險評價時使用最多的一個軟件平臺[7，18-19]。R軟件強(qiáng)大的統(tǒng)計軟件包(如fitdistrplus)能夠為SSD的建立、HC5的估計等提供更多自由全面的選擇，但該軟件的熟練運用需要使用者具有一定的概率統(tǒng)計學(xué)功底。本文給出的Matlab函數(shù)“PERA”只需要使用者具有污染物的毒性數(shù)據(jù)就能夠方便地得到SSD曲線和其重要參數(shù)HC5，為概率統(tǒng)計學(xué)功底較弱的環(huán)境工作者提供了技術(shù)支持。

目前關(guān)于生態(tài)風(fēng)險評價的軟件主要集中在SSD的構(gòu)建和環(huán)境標(biāo)準(zhǔn)(HC5)的推導(dǎo)，如CSIRO開發(fā)的BurrlizO[9]，U.S.EPA開發(fā)的SSD Generator和Kon等開發(fā)的MOSAIC_SSD[10]等。這些軟件雖然可以讓使用者方便地得到SSD和HC5，但并沒有與污染物的環(huán)境監(jiān)測數(shù)據(jù)相結(jié)合，無法進(jìn)行污染物的概率生態(tài)風(fēng)險評價。本文給出的Matlab函數(shù)“PERA”只需要使用者具有污染物的毒性數(shù)據(jù)和環(huán)境監(jiān)測數(shù)據(jù)，就可以方便地得到聯(lián)合概率曲線JPC和污染物的總體風(fēng)險概率ORP。雖然RIVM開發(fā)的軟件ETX 2.0也可以給出JPC和ORP，但該軟件中的SSD模型只能選Log-normal[8]，而本文所給函數(shù)中的SSD模型有Log-normal和Log-logistic兩種選擇。另外，本文所給Matlab函數(shù)的開源代碼，使用者可以對其作適當(dāng)增加或修改(如增加或修改SSD模型的概率分布類型等)，較ETX2.0更具靈活性。目前關(guān)于SSD的軟件主要是基于Excel表格和專業(yè)統(tǒng)計軟件R開發(fā)的，尚未有基于Matlab的程序或函數(shù)，本文首次給出了基于SSD的概率生態(tài)風(fēng)險評價的Matlab函數(shù)“PERA”，為污染物的生態(tài)風(fēng)險評價提供進(jìn)一步的技術(shù)支持。