中华白海豚Sousa chinensis又名印度太平洋驼背豚，为沿岸河口定栖性小型齿鲸类，属海豚科Delphinidae驼背豚属Sousa。主要分布于西太平洋和印度洋的热带和温带水域，常见于我国东南沿海(Ross et al.，1994；Jefferson & Karczmarski，2001)，其中以珠江口水域数量最多，约有1500头(Jefferson & Hung，2004)。早在1988年，中华白海豚就被列为国家Ⅰ级重点保护濒危野生动物，也是中国海洋哺乳动物中唯一的国家Ⅰ级重点保护物种。近年来，由于珠江口区域经济的高速发展，所带来的环境污染问题也日趋严重。据最近发布的《2013年广东省海洋环境质量公报》，珠江口水域在广东省各水域中污染最为严重。研究显示，珠江口水域中华白海豚种群正在走向衰落，其数量正以每年2.46%的速率下降(Huang et al.，2012)。因此，了解珠江口水域中华白海豚种群分布及其与环境因子的关系对于我国中华白海豚的保育工作至关重要。

栖息地选择是指动物为了达到生存目的(如觅食、休息、迁徙、繁殖、躲避敌害等)而寻找相对适宜的栖息地的过程，是动物对环境的一种行为适应(刘宣，2007)。动物栖息地选择的研究一直以来是动物生态学的研究热点，亦是开展珍稀濒危物种研究及生物多样性保护的基础(Johnson，1980；Caughley，1994)。作为适应海栖生活的海洋哺乳动物，其栖息地的选择受到生态学家极大关注。研究表明，栖息地地形、食物、盐度、深度、离岸距离、叶绿素含量等因子均会对海洋哺乳动物的分布格局产生影响(Hastie et al.，2003，2004；Fazioli et al.，2006；陈炳耀等，2007；Vázquez Castán et al.，2007)。中华白海豚具有数量少、活动性强、遇见率低的特点，因此对其栖息地选择的研究尤为困难。目前国内外关于中华白海豚栖息地选择的研究主要集中在其出现频次与水深、水温、盐度、潮汐变化、离岸距离、风力、海水透明度的关系(Karczmarski et al.，2000；Atkins et al.，2004；翟飞飞，2006；陈炳耀等，2007；刘宣，2007；陈涛等，2010)，但大多局限于传统的经验性描述、单因子分析及简单的线性回归分析。

物种的空间分布是多种环境因子共同作用的结果，简单的单因子分析和线性回归分析并不能很好地解释其栖息地选择机制。在众多统计方法中，广义可加模型(generalized additive models，GAM)作为近年来生态学空间格局分析的重要方法之一，它能应用非参数的方法检测数据结构，找出其中规律，进而得到更好的预测结果(Hastie & Tibshirani，1986)。由于具有能直接处理响应变量与多个解释变量之间非线性关系的优点，因此GAM被广泛应用于探究种群分布的时空特征及其与环境因子之间的关系的研究(牛明香等，2010)。目前国内外主要利用GAM分析渔业资源与环境因子的关系，对于鲸豚类的研究还比较少。Forney(2000)通过对短吻真海豚Delphinus delphis和鼠海豚Phocoena phocoena的栖息地选择进行GAM建模，发现GAM能有效降低环境变异引起的不确定性，是研究鲸类分布与环境因子关系的良好工具。本文通过截线调查法监测珠江口水域中华白海豚密度分布状况，并通过监测其栖息地范围内的10类海水水质因子(pH、深度、表层水温、盐度、溶解氧、叶绿素a、化学需氧量、活性磷酸盐、无机氮、活性硅酸盐)，首次运用GAM分析了珠江口中华白海豚的栖息地选择同海水因子的关系。研究结果有助于进一步揭示中华白海豚的栖息地选择机制，并为更全面了解中华白海豚的分布状况、栖息地评价及保育工作提供科学依据。

本研究区域为中华白海豚分布最密集和拥有资源数量最多的区域——珠江口水域，属南亚热带海洋性季风气候，光照充足，雨量充沛(图 1)。研究区域包括珠江口中华白海豚国家级自然保护区，北起中山与深圳对开水域，南至横琴、三灶和大澳岛，西至淇澳岛，东至香港大屿山(21.948 871°~22.463 839°N，113.450 360°~113.674 215°E)，总面积达700 km²。

图 1 研究区域及截线调查站位设置分布图 Fig. 1 Research area and line-transect survey sampling sites in Pearl River Estuary

图选项

珠江口中华白海豚的密度分布数据来自本实验室长期船基截线调查和标记重捕的监测结果，调查样线的设置综合考虑了海深、海岸、岛礁、避免重复等因素，并尽可能遵循分布的均匀性和随机性等概率原则。在适宜的监测条件下(蒲福海况0~3级，天气晴朗，能见度大于1200 m)，研究船只以13~15 km·h^-1的航速沿调查路线航行。由2名研究人员采用2台长焦相机(佳能5D mark Ⅱ)同时观察，一旦发现中华白海豚，使船只尽量靠近海豚，并进行拍照和摄像。调查过程中记录观测的时间、位置、航速、航向、能见度和航程，并记录中华白海豚目击时间、位置、头数等。位置、航速和航程信息由GPS(Garmin 78S)获得，调查时间为2013年2月—2014年9月，共进行104次航行，平均每月出海调查5.2次。根据船基截线调查和照相识别所记录的中华白海豚目击数据，采用GIS ArcView 3.1绘制中华白海豚目击位置分布图，并通过核密度分析，将每个目击点的中华白海豚目击数据散布到整个研究区域，以获得珠江口水域中华白海豚的密度分布数据。

本研究在整个珠江口水域设置了20个站位(图 1)，基于条件限制，本研究并未对整个船基截线调查路线进行连续取样监测，而是将站位的设置均匀分布于船基截线调查的路线中。分别于2013年12月24—25日、2014年4月26—27日和2014年9月4—5日在各站位采集水样。其中水样用2.5 L采水器采集，统一采集深度1 m处的海水。每一站位水样均取2份，保存于聚乙烯瓶中，并及时带回实验室测定。共测定珠江口水域的10类海水水质因子：pH、深度、表层水温、盐度、溶解氧、叶绿素a、化学需氧量、活性磷酸盐、无机氮、活性硅酸盐。其中，深度由水深计(SD5000)测定，pH、表层水温和盐度由YSI多参数水质分析仪(YSI professional plus，美国维赛)测定。叶绿素a、溶解氧、化学需氧量、活性磷酸盐、无机氮和活性硅酸盐的分析测定均按照《海洋监测规范》(GB 17378—2007)执行。水质评价标准依据《海水水质标准》(GB 3097—1997)。

采用GAM分析中华白海豚密度分布同水质因子的关系。GAM的表达式为：

其中，muY为响应变量的期望值；g()是连接函数；a₀是截距；ε为服从正态分布的随机误差项；f_n()是解释变量X_n的单变量函数，假设每个f_n均有光滑函数，然后利用散点图进行逐步平滑估计。模型中的误差分布为正态分布，所使用的散点图平滑为薄板样条(TPRS)平滑。本研究以中华白海豚密度分布为响应变量，水质因子为解释变量。利用正态Q-Q图和正态P-P图判断中华白海豚密度数据的分布型以大致确定其连接函数。在GAM的非线性分析中，需要克服解释变量之间的共曲线性关系导致的系数标准误的偏差，对于水质因子之间可能存在的共线性关系，采用Pearson相关系数来判断两者的相关性，当两者存在共线性关系时，只选取一个变量为解释变量(李增光等，2012；欧阳芳，戈峰，2013)。本研究采用逐步回归法检验模型中每个因子的重要性及遴选最优模型。首先，对单一水质因子做GAM，以确定中华白海豚密度分布和单一水质因子的关系，然后基于AIC准则(Akaike’s information criterion)检验逐步加入水质因子后模型的拟合 程度，AIC值越小，表明模型的拟合效果越好(Akaike，1974；Hastie & Tibshirani，1990)。 采用F检验评估各解释变量的显著性。根据估计自由度EDF(estimated degrees of freedom)和参考自由度RDF(referred degrees of freedom)判断解释变量与响应变量是否属于非线性关系。正态Q-Q图和正态P-P图的绘制通过SPSS 17.0实现，Pearson相关分析、GAM的构建及检验通过R-project 3.1.2中的mgcv软件包实现(Wood，2006)。

在2013年2月—2014年9月期间共记录中华白海豚327群1649头，各站位中华白海豚密度分布范围在0.02~5.26头/次之间。中华白海豚的分布主要集中在站位6、7、9、10、11、12、13、14、15、18(图 1虚线内部分)，呈现由河口向外海先增加后降低的趋势。其中，中华白海豚密度最高的为站位13，最低的为站位20。

基于2014年4月水质数据绘制的珠江口水域各站位10类水质因子分布态势图(图 2)。各站位深度、盐度均表现为由西北至东南呈递增趋势；各站位化学需氧量、活性磷酸盐、无机氮和活性硅酸盐均表现为由西北至东南呈下降趋势；珠江口水域各站 位pH、表层水温从北至南呈缓慢递增趋势；各站位溶解氧含量由北至南呈缓慢下降趋势；叶绿素a含量在各站位变化较大，其中叶绿素a含量较高的站位均为近岸或近岛礁分布，尤以珠江口西侧水域含量较高。其他月份水质因子的分布态势与4月类似。

图 2 珠江口水域各站位10类水质因子分布态势 Fig. 2 Distribution of ten water quality parameters in all sampling sites in Pearl River Estuary 图中每个站位圆的大小分别表示各水质因子含量的高低。 The size of the circle represents the content of water quality parameters in each site.

图选项

根据中华白海豚密度分布绘制的正态P-P图和正态Q-Q图(图 3)，散点基本落在对角线附近，认为中华白海豚的密度分布特征为正态分布类型，因此GAM以Identity link作为连接函数。Pearson相关分析表明，pH和深度，无机氮和叶绿素a，溶解氧和温度，盐度和硅酸盐之间存在显著相关性(R>0.5)，因此，本研究只选取其中一个作为解释变量。

图 3 中华白海豚密度分布的正态P-P图(左)和正态Q-Q图(右) Fig. 3 Normal P-P plot and normal Q-Q plot of the dolphins density distribution

图选项

本研究综合了2013年12月、2014年4月和2014年9月的水质分析结果，对于一组具有多个解释变量的数据集，构建模型时可构造出一个或多个可行或可能的模型。例如：数据集中有i个解释变量，则以不同解释变量组合构建的模型组包括可能模型的数量为：n=C¹i+C²i+…+C^i-1_i+Ci_i(欧阳芳，戈峰，2013)。因此，本研究先对单一水质因子建模，然后基于AIC准则检验逐步加入水质因子以获取最优模型。当对单一水质因子建模时，深度因子较其他因子具有更低的AIC值和更高的偏差解释率，因此在此后的模型构建中在深度因子的基础上逐步加 入其他水质因子。在模型构建中每生成一个GAM，均使用卡方检验与上一最优模型进行比较，从而获得新的最优模型，在R中的语句为“anova(模型n-1，模型n，test="Chi")”。

具体构建过程如下：

模型1 g(密度分布)~s(深度)，其中AIC=204.8429；

模型2 g(密度分布)~s(深度)+s(化学需氧量)，其中AIC=195.0984；

……

模型n g(密度分布)~s(深度)+s(化学需氧量)+s(活性磷酸盐)+s(盐度)+s(表层水温)+s(无机氮)，其中AIC=161.1911；

表 1为GAM统计结果。表中的每一行显示了逐步加入每一解释变量的拟合过程和解释变量加入后的相关信息。逐步回归分析表明，该模型对中华白海豚密度分布的总偏差解释率达到86.3%。其中贡献最大的为深度，贡献率为42.0%。其次为化学需氧量、活性磷酸盐、盐度，贡献率分别为12.7%、9.2%、18.2%。表层水温和无机氮的影响较低，对偏差的解释率仅为2.2%和1.7%；F检验表明，深度、化学需氧量、活性磷酸盐、盐度这4个海水因子对中华白海豚的密度分布影响显著(P<0.05)，表层水温和无机氮对其分布影响不显著(P>0.05)；深度、化学需氧量、活性磷酸盐、盐度、温度的估计自由度EDF均大于1，表明它们与中华白海豚密度分布之间存在非线性关系。无机氮的估计自由度等于1，表明其与中华白海豚密度分布之间为线性关系；随着海水因子的逐步加入，模型中的AIC值继续变小，表明模型具有很好的拟合度。

由表 1可知，深度对中华白海豚密度分布的影响极其显著，是影响中华白海豚密度分布的最主要因子，对整个模型的偏差解释率达到42.0%。随着深 度的增加，中华白海豚的密度分布呈先增加后降低的趋势，其主要分布在15 m以内，以5~10 m之间分布最多(图 4：a)，这与国内外其他研究结果相一致(刘宣，2007；陈涛等，2010)。调查表明珠江西部河口中华白海豚主要分布于深度小于10 m的水域，并认为20 m等深线是其分布的临界深度。厦门水域中华白海豚出现频次与水温、风向、风力、潮汐变化没有显著相关性，偏爱于1 km以内、5~20 m之间的近岸浅水生活(翟飞飞，2006；陈炳耀等，2007)。在南非阿尔戈亚海湾，中华白海豚主要分布于近岸2 km以内、深度小于15 m的水域，并认为其对潮水涨落、表层水温、海水透明度均没有明显选择性(Karczmarski et al.，2000；Atkins et al.，2004)。一方面合适的栖息地深度能为其提供合适的生存活动空间，另一方面珠江口近岸浅海区，大量营养物质随径流入海，使得众多浮游生物繁殖生长，成为鱼类活动的密集中心，从而为中华白海豚提供了丰富的饵料。

图 4 海水因子与中华白海豚密度分布关系的广义可加模型分析图 Fig. 4 The relationship between impact factors and the dolphin’s density distribution based on the generalized additive models

图选项

化学需氧量是监测水中还原性有机污染物的重要指标，能够在一定程度上反映出水体的污染程度。依据国家《海水水质标准GB 3097—1997》，一类水质的化学需氧量须小于2.0 mg·L^-1。由表 1可知，化学需氧量对中华白海豚密度分布的影响极其显著，对整个模型的偏差解释率为12.7%。中华白海豚主要分布在化学需氧量小于2.0 mg·L^-1的一类水质区域(图 4：b)。当化学需氧量小于1.0 mg·L^-1时，中华白海豚的密度与化学需氧量呈正相关；当化学需氧量大于1.0 mg·L^-1时，随着化学需氧量的增加，中华白海豚的密度则逐渐下降。中华白海豚分布与化学需氧量的相关性提示环境污染对其密度分布造成了一定的影响。

中华白海豚的密度分布与水中无机氮和活性磷酸盐含量均呈负相关(图 4：c，f)。其主要分布区的磷酸盐含量在0.02~0.04 mg·L^-1之间，无机氮的含量在1~2 mg·L^-1之间，均超过国家一类水质标准。无机氮和活性磷酸盐能为浮游生物和藻类的生长繁殖提供必不可少的微量元素。但高浓度的氮磷含量会导致水体的富营养化，由此带来的毒赤潮会严重威胁海洋生态的安全。研究表明，甲藻产生的甲藻毒素能够对座头鲸Megaptera novaeangliae造成致命的威胁(Geraci et al.，1989)。 根据化学需氧量、无机氮和活性磷酸盐的含量，采用改良的营养状态指数法(张景平等，2009)得出珠江口水域中华白海豚栖息地的富营养状态柱状图(图 5)。当营养状态指数E>1时，该水域呈富营养化状态，珠江口水域20个站位中有10个站位的E大于1，富营养化率为50%。珠江口水域中华白海豚栖息地面临富营养化的威胁。

图 5 珠江口水域各站位营养状态指数 Fig. 5 珠江口水域各站位营养状态指数

图选项

盐度对中华白海豚密度分布的影响也呈显著水平。当盐度小于10‰时，中华白海豚的密度随着盐度的增加而增加；在盐度为10‰~30‰之间的水域，随着盐度的增加，其密度呈缓慢下降的趋势；当盐度大于30‰时，其密度则大幅降低(图 4：d)。Jefferson(2000)通过生态调查发现，在盐度大于22‰的水域，珠江口中华白海豚的目击率与盐度的增长呈负相关。盐度是影响水生生物生长的重要环境因素之一，每种水生生物都会生存在适当的渗透压和离子浓度下。在等渗环境下因盐度差而额外付出的代谢能最少，更有利于水生生物的生存，而长期处于非等渗环境中可能导致其出现水肿、脱水、皮肤溃烂甚至死亡的情况。由GAM分析可知，中华白海豚主要分布于盐度为10‰~30‰之间的水域，对盐度变化有很大的适应性，属广盐性生物。珠江水系的径流量是仅次于长江的中国第二大河流，每年大量的冲淡水及营养物质经珠江口汇入南海，使得咸淡水交汇处具有丰富的鱼类资源。这种对盐度的广适性大大增加了中华白海豚对饵料获取的几率，更有利于其生存。

本研究首次利用GAM探究了中华白海豚密度分布与各种海水因子的关系，结果表明，GAM能够很好地揭示中华白海豚的密度分布与海水因子间的非线性关系，并具有很高的偏差解释率。与传统回归方法相比，GAM更具有灵活性和普适性。GAM研究表明，去除具有共线性关系的海水因子后，与中华白海豚密度分布具有显著相关性的海水因子为深度、化学需氧量、活性磷酸盐和盐度。鉴于海水因子对中华白海豚栖息地选择的重要性，在进行栖息地规划及评价时，充分考虑其栖息地环境特征尤其是海水因子特征十分必要。另外，通过监测中华白海豚栖息地海水因子特征，了解其栖息地生存环境状况，能够为制定相应的保育措施提供科学指导。

中华白海豚对栖息地的选择是多因子共同作用的结果，本研究主要探究其分布与海水因子的关系，除此之外，食物、地形和船只噪声也被认为是影响鲸豚分布的重要因素，因此在之后的研究，将考虑加入更多相关因子来完善该模型，以便为中华白海豚栖息地的选择和分布的研究提供更多依据。此外，由图 4可知本研究构建的GAM仍具有较大的离差，可能是由于不同时间的监测结果具有较大的离散度，也可能存在数据的不足，因此在之后的研究中会严格控制不同批次的监测流程，并充分积累数据。另外，由于条件限制，本研究的站点设置覆盖范围较窄，在之后的研究中有必要适当增加站点，以期更全面地探究影响中华白海豚分布的机制，同时也为其他海域的同类研究提供科学借鉴。