由于全球变暖的持续发展和社会经济的快速发展，世界很多国家或地区的水资源处于脆弱状态，水资源紧缺形势非常严重。因此，水资源脆弱性是国际社会普遍关心的全球性问题，也是中国可持续发展面对的重大战略问题，水资源脆弱性评价是世界水文水资源研究领域的重要前沿课题^[1–5]。

评价水资源脆弱性的常用方法是综合指数法，即首先建立相应的指标体系，其次对各指标进行赋权，然后利用一些数学模型综合成脆弱性指数^[6]，其中加权综合法和模糊综合评价法是最常用的两种数学模型^[1,7–8]。Wu等^[9]利用加权综合法评价地下水的脆弱性；Sullivan^[10]从供需角度建立社会经济和环境因子并确定指标权重，然后利用加权综合法计算水资源脆弱性指数；Dixon^[11]采用模糊综合评价法对地下水脆弱性进行了评价，并作出脆弱性图。综合指数法具有很好的研究基础，考虑脆弱性影响因素比较全面，但该方法在应用时也存在一些不足，如在赋权过程中主观性较强^[12]；除此之外，加权综合法和模糊综合评价法均是一种线性的加权评价方法，而评价的本质是非线性的^[13]。为弥补线性加权法的不足，一些非线性评估方法逐渐被应用到脆弱性或风险评价中，如分形理论^[14]、集对分析理论等^[15]、模糊物元模型^[16]及非线性模糊综合评价^[17]及函数法等^[18]。非线性评估模型的提出在一定程度上推动了水资源脆弱性的研究，但是仍然存在一定的局限，如需要对指标进行赋权，仍然无法避免主观性的影响；另外大部分模型需要确定指标与脆弱性间复杂的数学关系式，如非线性模糊综合评价法需要确定隶属函数的表达式等，但这些关系式随研究地区或研究内容的不同而需作相应的改变，不便于推广^[19]。

基于上述讨论，作者在对水资源脆弱性指标进行降维处理的基础上，分析水资源脆弱性与指标变量之间的变化关系，根据微分方程建模技术构建一种水资源脆弱性评估的通用模型，该模型不需要对指标进行赋权，可以避免人为的干扰。

1 水资源脆弱性非线性评估模型构建

水资源脆弱性评价模型的建模与计算步骤如图1所示，详细建模过程见1.1～1.4节。

1.1 水资源脆弱性指标及标准化处理

关于水资源脆弱性的定义，学术界还没有统一的定论，本文采用文献[3]中水资源脆弱性的定义：水资源脆弱性是指水资源系统在面对潜在危险时所表现出的特殊的性质或状态，主要包括降水量P、人均水资源量W_p、水资源满足程度S_r、水资源开发利用率U_r和污水处理率D_r。其中水资源开发利用率^[20]和水资源满足程度的计算公式分别为：

${U_{\rm{r}}} = \frac{{W_{\rm{u}}}}{W}$

(1)

式中，W_u为流域或区域已开发利用的水资源量，W为水资源总量。

${S_{\rm r}} = \frac{{{W_{\rm as}}}}{{{W_{\rm td}}}}$

(2)

其中，W_as为可供水量，W_td为总需水量。所谓指标标准化处理，是指保持指标的同趋势化，以保证指标间的可比性。而指标往往分成以下类型：成本型、效益型、适度型和区间型，所谓成本型指标是指指标数值越小越好的指标，效益型指标是指数值越大越好的指标，适度型指标是指数值越接近某个常数越好的指标，区间型指标是指数值越接近某个区间（包括落在该区间）越好的指标。因此，降水量P、人均水资源量W_p、污水处理率D_r和水资源满足程度S_r是成本型指标，水资源开发利用率U_r是效益型指标。由于成本型指标较多，作者拟将所有指标处理成越小越好的指标，即V随所有指标变量的增大而减小。标准化方法如下：

$\left\{ \begin{array}{l}{y_{ij}} = \frac{\displaystyle{\mathop {\min }\limits_i \left\{\displaystyle {{a_{ij}}} \right\}}}{{{a_{ij}}}} \,\,\,\,\,{\text{(效益型)}}\text{，}\\[5pt]{y_{ij}} = \frac{\displaystyle{{a_{ij}}}}{\displaystyle{\mathop {\max }\limits_i \left\{ {{a_{ij}}} \right\}}} \,\,\,\,\, {\text{(成本型)}}\end{array} \right. $

(3)

其中，a_ij为指标的原始值，y_ij为指标标准化后的值。因此，经过标准化处理后，V随所有指标变量的增大而减小。

1.2 指标的降维处理

由于影响水资源脆弱性的变量有很多，定量分析它们之间的函数关系非常困难。可以利用降维的思想将多维指标将转换成1维指标，从中筛选和提取特征信息，从而使进一步研究脆弱性与变量之间的函数关系变得简单。投影寻踪是一种对高维数据进行降维的方法，其基本思想和建模步骤参见文献[21–22]，主要包括构造投影函数和估计最优投影方向两个步骤。投影函数一般为线性函数，表达式如下：

$x(i) = \sum\limits_{j = 1}^n {{\textit{ω}_j}{y_{ij}}} $

(4)

式中， $\left\{ {{y_{ij}}|i = 1, 2, \cdots, m;j = 1, 2, \cdots, n} \right\}$ 为标准化后的数据序列，m、n分别为样本容量和指标个数， ${\mathit{\boldsymbol{\omega}}} = \left( {{\textit{ω} _1},{\textit{ω} _2},} \right.$ $\left. { \cdots ,{\textit{ω} _n}} \right)$ 为投影方向。

在文献[21–22]中，要求投影方向为单位向量，但没有给出具体理由。本文以投影方向和为1作为优化投影方向的约束条件，于是可构造如下的最小化问题来估计最优投影方向。

$\begin{array}{l}\min \;\;\;\;Q\left( \textit{ω} \right) = - \frac{\displaystyle{\sum\limits_{i = 1}^m {{{\left[ {x\left( i \right) - \bar x} \right]}^2}} }}{\displaystyle{m - 1}},\\{\mathop{\rm s.t}\nolimits} \;\;\;\;\;\displaystyle\sum\limits_{j = 1}^n {{\textit{ω} _j} = 1} \end{array}$

(5)

显然，式（5）是一个条件极值问题，可以构建一个拉格朗日函数^[23]：

$L\left( {\textit{ω} ,\lambda } \right) = Q\left(\textit{ω} \right) + \lambda \left( {\sum\limits_{j = 1}^n {{\textit{ω}_j}} - 1} \right)$

(6)

$\begin{aligned}[b]\frac{\displaystyle{\partial L}}{\displaystyle{\partial {\textit{ω} _j}}} = & - \frac{\displaystyle{2}}{\displaystyle {m - 1}} \displaystyle \sum\limits_{{i} = 1}^m {\left( \displaystyle{\sum\limits_{{k} = 1}^n {{c_{ik}}{\textit{ω} _k}} } \right)} {c_{ij}} + \lambda =\\[7pt]& - \frac{\displaystyle{2}}{\displaystyle{m - 1}}\displaystyle \sum\limits_{{k} = 1}^m {\left( {\textstyle \displaystyle \sum\limits_{i = 1}^m {{c_{ij}}{c_{ik}}} } \right)} {\textit{ω}_k} + \lambda = 0\end{aligned}$

(7)

式中， $ \left( \displaystyle{{y_{ij}} \! - \! \frac{1}{n}\sum\limits_{i \! = \! 1}^m {{y_{ij}}} } \right) \! = \! {c_{ij}},{y_{ik}} - \displaystyle \frac{1}{n}\sum\limits_{k \! = \! 1}^n {{y_{ik}}} = {c_{ik}},j \! = \! 1,2,\cdots , n$ 。

因此，只要方程（7）的系数矩阵的行列式不为0，方程（7）有一个唯一解 $\textit{ω} = \left( {{\textit{ω} _1},{\textit{ω} _2}, \cdots ,{\textit{ω} _n}} \right)$ 。

1.3 水资源脆弱性微分方程模型

由于经过标准化处理后，V随所有指标变量的增大而减小，将降水量P、人均水资源量W_p、水资源满足程度S_r、水资源开发利用率U_r和污水处理率D_r投影成一维变量后x后，根据式（4）可知V也随着x的增大而减小，因此V是x的减函数。

张晓慧等^[13]认为一个好的评估函数应该具有如下性质：1）评估函数应该具有单调性，包括单调增加或单调减小，如果是单调增加，则评价指标值越高，评估值越大；2）评估函数应该是有界且平稳的，平稳性是指：评价指标值的变化应使得评估结果连续平稳的变化，没有跳跃。钱龙霞等^[3]认为评估函数的变化除了是连续平稳的，还应该满足性质3）：先是越来越快，到达某一拐点后越来越慢。

性质1）和2）说明评估函数应该是单调连续的；根据导数理论和拐点理论^[23]可以将性质3）用数学的语言来描述，即导函数开始是增函数，到达某一拐点后是减函数，且导函数是连续函数。根据上面的分析，可以提出如下假设：

1）V是连续函数，且自变量x的定义域为[a,b]；

2）V的导函数是连续的，即V是光滑函数；

3）V的导函数在[a,c]是增函数，在[c,d]上是减函数，即c为V的拐点；

为了保证水资源脆弱性的可比性，提出假设4）：V的取值范围为[0,1]。由于V是x的减函数，则V（a）=1，V（b）=0。为了简化问题，假设导函数是对称的，即c为a和b之间的中点。

根据以上假设可以构建如下的微分方程模型

$\frac{{{\rm d}V}}{{{\rm d}x}} = \alpha {\left( {x - c} \right)^2} + h,a \le x \le b$

(8)

1.4 水资源脆弱性非线性评估模型

求解方程（8）可得水资源脆弱性非线性评估模型为：

$V = \frac{\alpha }{3}{\left( {x - c} \right)^3} + hx + \beta ,a \le x \le b$

(9)

根据假设1）～3）及水资源脆弱性的特点可知方程（9）满足：

$\left\{ \begin{array}{l}h < 0,\\\alpha {\left( {a - c} \right)^2} + h < 0,\\\alpha {\left( {b - c} \right)^2} + h < 0\end{array} \right.$

(10)

$\left\{ \begin{array}{l}\frac{\displaystyle{\alpha }}{\displaystyle{3}}{\left( {a - c} \right)^3} + hx + \beta = 1,\\[7pt]\frac{\displaystyle{\alpha }}{\displaystyle{3}}{\left( {b - c} \right)^3} + hx + \beta = 0\end{array} \right.$

(11)

方程组（11）是关于 $ \alpha $ 、h和β的非齐次线性方程组，根据线性代数中线性方程组解的理论、不等式组（10）及指标的数据样本资料即可获得方程（9）中的未知参数，进而获得水资源脆弱性V与投影变量x之间的非线性评估模型。

2 实例研究 2.1 研究区概况

北京位于华北平原西部，主城区分布在永定河、北运河冲击平原上，偏远区县主要分布在大清河、永定河、北运河、潮白河山区。北京属暖温带半干旱半湿润性季风气候，雨量年际和季节分配极不均匀，主要集中在夏季，全市多年平均降水量仅约为570 mm（1956—2012年）。北京水资源严重紧缺，人均水资源量约200 m³（以2008年为例），地表水资源开发率高达86%，同时地下水处于严重超采状态。另一方面，由于人口的快速增长，北京生活用水需求逐年增加，使得水资源供水压力进一步增大。水资源的高度开发利用引发一系列生态问题，如河流断流、水质恶化及地面沉降等，水资源系统受到很大影响。

2.2 水资源脆弱性模型构建与验证 2.2.1 水资源脆弱性函数模型的建立

首先对1979—2012年的水资源脆弱性指标进行标准化处理（式（3）），将处理后的数据代入方程（7）可得系数矩阵的行列式为–2.256 3，所以方程（7）有唯一解， $ \textit{ω} $ =(0.737 467, −0.001 079, 0.897 492, 0.009 677, −0.645 72)。

因为

$\begin{aligned}x\left( i \right) =& \sum\limits_{j = 1}^5 {{\textit{ω} _j}y(i,j)} =\\& 0.737\,\,\,467 \times y\left( {i,1} \right){\rm{ - }}0.001\,\,\,079 \times y\left( {i,2} \right)+\\& 0.897\,\,\,492 \times y\left( {i,3} \right) + 0.009\,\,\,677 \times y\left( {i,4} \right)-\\& {\rm{ }}0.645\,\,\,72 \times y\left( {i,5} \right)\text{。}\end{aligned}$

而对 $\forall i,j$ ，有 $0 \le y(i,j) \le 1$ ，

因此

$- 0.645\,\,\,72 \le \sum\limits_{j = 1}^5 {{\textit{ω} _j}y(i,j)} \le 1.645\,\,\,716\text{。}$

再根据1.2节中假设1）可以确定a和b分别为-0.645 72和1.645 716，由此可得c为0.5。解线性方程组（11）可得：

$\left[ \begin{array}{l}\alpha \\h\\\beta \end{array} \right] = \left[ \begin{array}{l} - 3.822\,\,\,8\\\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,1\\\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,0\end{array} \right] + \left[ \begin{array}{l} - 4.570\,\,\,9\\\,\,\,\,\,\,\,\,\,- 2\\\,\,\,\,\,\,\,\,\,\,\,\,\,\,1\end{array} \right]k$

(12)

其中，k为一切实数，联立式（10）和（12）求解可得 $0.718\,\,\,196 \!<\! k \!<\! 0.827\,\,\,292$ ，因此k为(0.718 196, 0.827 292)之间的一切实数，于是可取k为0.82，进而可求得 $\beta = 0.82$ ， $\alpha = 0.419\,\,\,629$ ， $h = - 0.64$ ，于是可得水资源脆弱性评估函数为：

$V = 0.139\,\,\,9 {\left( {x - 0.5} \right)^3} - 0.64x + 0.82$

(13)

2.2.2 北京市1979—2012年水资源脆弱性计算与分析

将1979—2012年的水资源脆弱性投影值数据代入方程（13）即可得到北京市1979—2012年的水资源脆弱性，计算结果如图3所示。其次利用Hierarchical Cluster对1979—2012年北京市的水资源脆弱性进行聚类，聚类结果显示1979—2012年北京水资源均处于强脆弱或极度脆弱状态，如图4所示。利用判别分析筛选影响北京水资源脆弱性的关键因子，筛选过程和结果见表1、2。

由表1和2可知，水资源开发利用率和污水处理率是影响水资源脆弱性的敏感因子。以水资源开发利用率和污水处理率的选择为例说明敏感因子的筛选过程，由表2第1步可知，水资源开发利用率的F值最大，大于spss软件默认的3.84，且Wilk’s Lamda值最小；再由表1第1步可知，水资源开发利用率移出的F值为29.434，大于spss软件默认的2.71，因此水资源开发利用率第1个进入模型。由表2第2步可知，水资源开发利用率进入模型后，污水处理率进入的F值最大，大于spss软件默认的3.84，且Wilk’s Lamda值最小；再由表1第2步可知，污水处理率移出的F值大于模型默认的2.71，因此第2个进入模型的是污水处理率。同理可以分析得出其他3个变量没有被选择，即它们只是影响水资源脆弱性的一般因子。

由图3和4可知，1979—2012年北京市的水资源脆弱性值介于0.77和0.85之间，均为强脆弱或极度脆弱。不仅如此，1979—1987年水资源脆弱性有增大的趋势，2004—2012年水资源脆弱性有减小的趋势，水资源脆弱性在1988—2004年间的波动比较小。由于水资源开发利用率和污水处理率是影响北京市水资源脆弱性的敏感因子（表2），而1979—1987年水资源开发利用率均较大，且污水处理率较低（均为10%左右），而且呈现下降的趋势，因此1979—1987年水资源脆弱性有增大的趋势。而2004年以来水资源开发利用率有降低的趋势，且污水处理率逐年增加（其中2012年已达到83%），这使得再生水利用量逐年增加，这在一定程度上扩大了北京市的水资源总量，而用水量也有下降的趋势，这使得水资源开发利用率又有所下降，因此2004年以后水资源脆弱性有下降的趋势。另一方面，水资源脆弱性是多种因素对水资源系统影响的综合表现^[2]，除了水资源开发利用率和污水处理率之外，其他因素如降水量、人均水资源量及水资源满足程度也会造成各年份水资源脆弱性的差异。需要强调的是，以上分析只能说明1979—2012年水资源脆弱性的趋势，并不能说明其它年份水资源脆弱性也遵循同样的规律。

以上分析说明本文建立的模型计算结果具有一定的合理性，可以用来预测未来的水资源脆弱性。

2.3 2020年北京市水资源脆弱性评价 2.3.1 水资源供需平衡分析

水资源供需平衡分析供需分析的结果为2020年在1956—2012年57种来水条件下的可供水量序列和需水量序列，来源于文献[3]，如图5所示。将供水、需水序列代到式（1）、（2）即可得到在57种来水条件下的水资源开发利用率和水资源满足程度等水资源脆弱性指标的数值。

2.3.2 水资源脆弱性计算与分析

当不考虑利用外调水和再生水时，将2020年在57种来水条件下的脆弱性指标数据序列进行标准化处理，代入式（7）求得系数矩阵的行列式为–0.001 7，因此方程（7）有唯一解 $\textit{ω} = \left( {0.994\,\,\,514,{\rm{ - }}0.001\,\,\,43,} \right.$ $\left. {0.370\,\,\,181,{\rm{ - }}0.368\,\,\,95,0.005\,\,\,688} \right)$ 。根据式（4）可以确定方程（3）里的参数a和b分别为–0.370 38和1.370 383，解线性方程组（11）可得：

$\left[ \begin{array}{l}\alpha \\h\\\beta \end{array} \right] = \left[ \begin{array}{l} - 6.234\,\,\,9\\\,\,\,\,\,\,\,\,\,\,\,\,1\\\,\,\,\,\,\,\,\,\,\,\,\,0\end{array} \right] + \left[ \begin{array}{l}7.920\,\,\,1\\\,\,\,\,\,\,\,- 2\\\,\,\,\,\,\,\,\,\,\,\,1\end{array} \right]k$

(14)

联立式（10）和（14）求解可得 $0.787\,\,\,225 < k < 0.930\,\,\,837$ ，于是可取k为0.93，进而可求得 $\beta = 0.93$ ， $\alpha $ = 1.130 793， $h = - 0.86$ ，于是可得水资源脆弱性函数为：

$V = 0.376\,\,\,931 {\left( {x - 0.5} \right)^3} - 0.86x + 0.93$

(15)

将2020年在57种来水条件下的脆弱性指标的标准化值序列代入式（4）可得到2020年在57种来水条件下的水资源脆弱性投影向量数据序列，并将其代入式（15）可得2020年北京市在不同来水条件下的水资源脆弱性，结果如图6（a）所示。

由图6（a）可知，当不考虑外调水和再生水时，2020年北京市在1956—2012年来水条件下的水资源脆弱性值介于0.882和0.883之间，水资源均处于极度脆弱状态。而1979—2012年的水资源脆弱性介于0.77和0.85之间，均为强脆弱或极度脆弱，说明2020年的水资源脆弱性结果仍然保持了历年水资源的高度脆弱状态，与过去的计算结果相吻合。另一方面，随着北京人口和经济的快速增长，2020年北京市用水量将会增加^[2]；而且北京面临水库储水量锐减、地下水位持续快速下降等问题^[2]，水资源的战略储备濒临枯竭，这些必然会导致2020年水资源的极度脆弱状态。

如果考虑南水北调水和再生水，相当于增加了水资源量，2020年人均水自流水资源量、水资源开发利用率和水资源满足程度均会发生变量，2020年外调水和再生水利用量来源于文献[3]。将相应数据代入式（4）、（7）、（9）、（10）和（11），同理可得水资源脆弱性函数为：

$V = - 0.005\,\,\,84 {\left( {x - 0.5} \right)^3} - 0.46x + 0.73$

(16)

将2020年利用外调水和再生水后的水资源脆弱性投影向量数据序列代入式（16）可得2020年利用外调水和再生水后在不同来水条件下的水资源脆弱性，如图6（b）所示。

由图6（b）可知，利用再生水和南水北调水后，北京市在不同来水条件下的水资源脆弱性均在0.73左右，平均下降幅度仅有17.2%，这是由于水资源脆弱性是水资源系统（承险体）在面对潜在危险时所表现出的特殊的性质或状态，是多种因素对水资源系统影响的综合表现^[2]，包括降水量、人均水资源量、水资源满足程度、水资源开发利用率和污水处理率等，利用再生水和外调水后，虽然某些指标如水资源满足程度和人均水资源量会增大、水资源开发利用率会减小，但2020年的降水量和污水处理率并没有变化。

3 结　论

本文建立了基于降维思想的水资源脆弱性非线性评估函数模型，该函数模型可以避免人为的干扰，具备单调性、有界性和变化的连续平稳性等性质，能定量模拟和刻画水资源脆弱性与指标之间的变化关系，较易推广。研究结果表明：1）水资源开发利用率和污水处理率是影响北京市水资源脆弱性的敏感因子；2）在1956—2012年的来水条件下，2020年北京市水资源均处于极度脆弱状态；利用再生水和南水北调水后，2020年北京市在不同来水条件下的水资源脆弱性均有所下降，平均下降幅度为17%。

由于本文研究建立在对指标进行降维处理的基础上，而在数据处理过程中难免会丢失一些信息，另外Uitto^[24]指出在风险分析中人们经常忽视了脆弱性的时变特性，即在脆弱性微分方程模型中要加入时间变量，这些问题都有待于今后进一步研究。