2.6.2　响应矩阵法（response matrix method）

仍以前面的例题为例，来对比介绍。

【例2.22】二单元潜水含水层水量分配问题（嵌入法[例2.17]）。

例题描述见[例2.17]，采用嵌入法建立的LP模型为

（1）无抽水情况下，二单元含水层天然地下水水位。由式（2.42），代入Q1=Q2=0、L=10000m、B=10000m、补给强度N=0.36m/a、导水系数T=3.6×106m2/a，将h1、h2用h01和h02代替，可以解出天然水位h01和h02（水位值从河水位算起）：

（2）抽水情况下，二单元含水层地下水水位。设抽水情况下，二单元含水层的水位降深为s1、s2，则实际含水层地下水位h1、h2可以表示为

（3）抽水响应矩阵确定。设βi，j表示j 井以单位抽水量单独抽水在i 点产生的降深，则按线性叠加原理，有

式中 m——抽水井数。

利用上式对节点i展开，有

将式（2.44）代入式（2.43）并考虑到水位约束限制条件，可有

这样改造后的LP模型如下：

这里的βi，j矩阵称为响应矩阵，下面来看如何确定βi，j矩阵中各元素值。

设单元1以单位抽水量抽水，单元2不抽水，注意到式（2.43）、式（2.44），这时式（2.45）有

代入h01=NLB/T、h02=2NLB/T，可从该式中可以解出β1，1、β2，1：

类似，设单元2以单位抽水量抽水，单元1不抽水，可以得到

这样响应矩阵为

上式中的响应矩阵由于数值较小，实际确定响应矩阵时，可对单位抽水量取某一数值，这里取I=106m3/a。此时的响应矩阵为

不过这么处理后，后面在形成式（2.47）时，别忘了式中Q 值应该除以I。

（4）应用响应矩阵对水位变量进行处理。响应矩阵代入水位或降深公式（2.45），有

或者写成：

进一步整理，有

如果采用I=106m3/a后得到的响应矩阵，则上式应写成：

该式整理后与式（2.47）是一样的。式中的0.1389 为I 脉冲抽水量形成的降深，0.1389/I 为单位脉冲抽水量形成的降深，0.1389Q1/I 为Q1抽水量单独抽水在点1形成的降深。

这里也可看出，脉冲抽水量I 的引入，对最终的方程没有任何实质影响，只是当时提出这一说法的人考虑到：这么处理下，响应矩阵的数值更像实际的抽水引起的降深值了（比如0.1389m 是可以实际观测到的。不然0.1389×10-6m 降深就不好理解了，这么小的降深实际上根本就观测不到）。但脉冲抽水量I 只是中间的一个替换概念，脉冲抽水量I的引入对最终方程的建立多少有些画蛇添足之举。

下面把通过响应矩阵处理后的规划模型与原嵌入法模型对照如下：

从上面对照中可以看出，通过响应矩阵βi，j，将嵌入法中与水位变量有关的式子通过响应矩阵变换成了另外的形式，去掉了水位变量，去掉了数值方程。

本例题中，去掉的数值方程虽然只有2个，但对于实际问题，去掉的是整个嵌入的数值方程。

（5）响应矩阵法模型求解。

【例2.23】用Excel求解结果见表2.11。

表2.11　[例2.23]Excel规划求解结果

表2.11中的水位h1、h2是利用式（2.45）回代求出的。以方案1为例：

或者采用脉冲抽水量I，回代计算为

通过上例，可归纳采用响应矩阵法建模的步骤为：

（1）确定建模中需要的节点初始水位（一般需利用数值方程求解）。

（2）利用数值方程，让每口抽水井单独抽水，分别求解响应矩阵各元素值。

（3）利用响应矩阵，根据水位与降深间的关系，将原嵌入法中含水位的所有变量从模型中消除（包括约束条件和目标函数中的水位变量，不包括非负的水位约束项）。

（4）求解不含水位（或降深）变量的LP模型。

（5）需要计算节点水位时，可根据水位与降深间关系，利用响应矩阵回代计算。

线性系统与叠加原理：

系统的作用是通过储存、传输、延时和平滑等作用，将系统的输入转换为输出。假设输入与输出均是时间t的函数，以x（t）表示输入、y（t）表示输出，则系统作用可表示为

式中h（t）或→代表系统的作用，它是反映系统特性的函数。地下水资源系统中的许多过程，如降水入渗、地下水位变化、抽水或回灌、地下水向地表水排泄等均是时间过程。

若系统输入x1（t）产生输出y1（t），输入x2（t）产生输出y2（t），而当输入为x1（t）+x2（t），输出为产生输出y1（t）+y2（t），即

则称这个系统可以应用叠加原理。线性系统必须满足叠加原理。但满足叠加原理的系统不一定是线性系统，它还需要满足倍比原理：如x（t）→y（t），有nx（t）→ny（t）。

满足叠加原理和倍比原理的系统才是线性系统，也称为齐次线性系统。否则称为非齐次线性系统。对于非齐次线性系统，可以从中分解出一个线性系统，这就是通常将地下水系统中某些变量的关系处理为线性关系的依据。

如y（t）=a x（t）+b 中，y（t）与x（t）显然可看成有线性关系，给出不同输入，有不同线性输出：

但当输入为[x1（t）+x2（t）]时，输出为a[x1（t）+x2（t）]+b≠y1（t）+y2（t），不满足叠加性。

叠加原理在地下水资源系统管理中有重大意义的原因在于，在求解若干项输入引起系统效应时，可以先单独分析各个输入产生的输出，然后将这些输出加以叠加以获得多项输入的总体效应。这就是响应矩阵法应用的理论基础。

【例2.24】二维稳定问题的响应矩阵法（嵌入法[例2.18]）。

[例2.18]采用响应矩阵法嵌入法给出的LP模型为：

按前述响应矩阵法建模的步骤来求解本题：

（1）确定建模中需要的节点初始水位。本题为承压含水层，抽水前初始水位水平，各节点均为20m（高于隔水顶板），不用求解。

（2）利用数值方程，让每口抽水井单独抽水，分别求解响应矩阵各元素值。本题中，共有6、7、10、11、14、15、18、19 计算节点，其中6、7、10 按要求各布置1 口抽水井，其余各点优化1 口抽水井。这样，这8 个点都是可能的抽水井，需要每口井单独计算。

设点6单独抽水，其余井不抽水，取脉冲抽水强度I=0.01m/d，即Q6=0.01m/d，其余Q=0，代入相关参数到上面数值方程中，用Excel求解结果如图2.26所示。

图2.26　Q6=0.01m/d，其余Q=0时Excel规划求解结果

图2.26的表中，Excel直接求解的是脉冲抽水下各点的水位值，初始水位（20m）与相应单元水位之差是剩余降深，即响应矩阵。按βi，j的定义，上边降深s 行剩余降深对应的βi，j依序为β6，6、β7，6、β10，6、β11，6、β14，6、β15，6、β18，6、β19，6，即单元6单位脉冲抽水时对点6、7、10、11、14、15、18、19的降深影响值。

类似，可以分别求出每口井脉冲抽水下的降深值，结果汇总到图2.27中。

按βi，j的定义，图2.27 的表中元素行列转置后，就是我们要求的响应矩阵表（图2.28）。

图2.28中，逐列求和列于最后一行，便于后面计算使用。

（3）利用响应矩阵，根据水位与降深间的关系，将原嵌入法中含水位的所有变量从模型中消除。从方程中可见，只有目标函数含有水位变量，需要利用响应矩阵及水位与降深的关系，计算消除水位变量。

图2.27　分别采用脉冲抽水后的降深Excel计算结果

图2.28　[例2.24]脉冲抽水形成的响应矩阵表

注意到图2.28中表末行数据，有

加上前面的水量约束和非负约束：

（4）求解不含水位（或降深）变量的LP 模型。取P=0.01m/d，用Excel求解上面的规划模型，解如图2.29所示。

(https://www.daowen.com)

图2.29　[例2.24]Excel规划求解结果（P=0.01m/d）

（5）需要计算节点水位时，可根据水位与降深间关系，利用响应矩阵回代计算。水位回代计算，如：

回代计算结果与嵌入法[例2.18]相同。

上面回代计算式中，如下加点说明，就会理解得更清晰了：

顺便说明一点，Excel计算中，可以在单元格格式中设置显示小数位数，结果显示较少的小数位数会看起来更清晰。而实际的精度不会因为显示位数的减少而降低。数据计算时，可以利用“=”方式引用单元格内的数据，这样不会因为舍入误差的影响，导致最后结果出现明显的比对误差。

【例2.25】潜水含水层三单元非稳定抽水最优分配（嵌入法[例2.19]）。

对于非稳定问题，响应矩阵法计算仍然采用前面的5个步骤：

（1）确定建模中需要的节点初始水位。

（2）利用数值方程，让每口抽水井单独抽水，分别求解响应矩阵各元素值。非稳定情况下，第i节点n 时段末的实际水位为

式中 m——抽水井节点数；

n——时段数，取具体某一值（1≤n≤最大时段数）。

如果采用脉冲抽水得到的响应矩阵，流量部分应该除以脉冲流量I，即

上式可理解为稳定流情况下水位计算公式的进一步扩展。

单位脉冲响应矩阵可用图2.30来表示，图2.30是j 点以单位抽水量从时段初开始连续抽水至n 时段，在i 点引起的水位降深及仅第1时段抽水，以后停抽在i 产生的水位恢复曲线。

图2.30　响应矩阵β（i，j，n-k+1）示意图

从恢复曲线看，当k=1时候，响应矩阵为β（i，j，n），表示第1时段j 点单位抽水量抽水在相隔n 时段对i点所产生的降深。而当k=n 时候，响应矩阵为β（i，j，1），表示第n 时段j 点单位抽水量抽水，在第n 时段末对i点所产生的降深。由此可见，k 越小，说明观测时间间隔越长，单位脉冲引起的剩余降深越小。k 越大，时间间隔越短，单位脉冲引起的剩余降深越大。因此，响应矩阵β（i，j，n-k+1）表示j 点第k 时段以单位抽水量抽水，在相隔（n-k+1）时段对i点产生的水位降深。从图2.30中可以看出，β（i，j，n-k+1）含义是剩余降深。特别应该注意的是：第k 时段抽水与间隔（n-k+1）时段后的降深的对应关系（图2.30 k=1与n-k+1=n、k=n 与n-k+1=1等的对应关系）。

图2.30中，假设非稳定过程分3个时段，则图中有3个β 值。最大降深β 值代表第1时段抽水，第1时段末的降深，也代表第2时段抽水，第2时段末的降深，还代表第3时段抽水，第3时段末的降深；图中停抽后的第二个降深β 值，代表第1时段抽水，第2时段末的降深，也代表第2时段抽水，第3时段末的降深；最后一个最小降深β 值，代表第1时段抽水，第3时段末的降深。

数值方程式（2.22）为

上式中，单元1为脉冲抽水量，其余Q 取零值，可以求出相应时段末各单元的水位h，然后利用s=h0-h 求出的s 就是相应单元抽水时的响应矩阵元素，Excel结果如图2.31所示。

图2.31　单元1取I=1000m3/d时计算水位与降深值

图2.31中最后一行就是响应矩阵各元素值，只是对应β（i，j，n-k+1）矩阵时，顺序要逆序进行调整。

类似取，其余Q 取零值，可以得到图2.31最后一行6个计算数值。取，其余Q 取零值，还可以得到图2.31最后一行6个计算数值。这3行与β（i，j，n-k+1）矩阵对应关系见表2.12。注意到上表中第2行，单元1第1时段抽水第2时段停抽时，单元3第2时段末水位还是继续下降的，这可能是如图2.32所示的情况。停抽后，井周水位回升快，远处水位继续向井流动，则表现出继续下降的情况。