
1 系统模型的建立
浅水湖泊垂向浓度分布近似均匀一致,从而使水质模型简化为二维问题.根据水平方向的浓度分布,将湖泊划分为n个单元,如图1所示的滇池被分为31个单元。对于某一单元,其水体浓度将由中心点的浓度来代表。系统滤波模型包括2个矩阵方程——状态方程和观测方程。
1.1 状态方程
对某一单元(i,j),如图2所示,根据质量平衡原理,时段dt间的水质变化应有如下关系:



图1 滇池测点分布及单元划分、水量交换(湿季)图
1—西山中;2—灰湾西;3—灰湾中;4—灰湾东;5—西华街中;6—观音山西;7—观音山中;8—观音山东;9—白鱼口;10—海口西;11—海宝口;12—大河尾;13—海口东;14—太史村;15—化肥厂
式中:各变量的整数下标i-1,i,i+1,j-1,j等,为一个单元的中心在坐标x,y上的位置,分数下标i-1/2,i+1/2,j-1/2等,为两个相邻单元的交界面的坐标x,y上的位置,x为水流的主方向,y为与x正交的流向;Vi,j为(i,j)单元的水体体积,它与四周相邻单元的界面面积为Fi-1/2,j,Fi+1/2,j,Fi,j-1/2,Fi,j+1/2;Q为通过某一界面的流量,C为某一污染物浓度,例如Ci,j为某种污染物在Vi,j中心的浓度,Ci-1/2,j为界面为Fi-1/2,j上的浓度,Dx、Dy分别为某种污染物沿x,y向的扩散系数,Δx,Δy为两个相邻单元中心在坐标x、y上的间距,Si,j为(i,j)单元的源强,Wi,j为(i,j)单元的模型计算浓度随机误差,称系统噪声,gi,j为单元系统噪声权重系数.

图2 计算单元水体示意
两相邻单元面上的浓度与该相邻单元中心点的浓度有关,可近似线性组合表示:

式中:α、β为界面污染物浓度权重系数。
生化需氧量BOD和化学需氧量COD降解,用一级反应动力学关系计算,则(i,j)单元的源强为

式中对于BOD,k=k1为BOD的降解速率系数;对于COD,k=k2为COD的降解速率系数;Ui,j为单元(i,j)的污染负荷强度,即河流等向单元水体的总排污率,排入为正,排出为负。将式(2)、式(3)代入式(1),整理后得(i,j)单元水质状态变化微分方程:

式中:

对于n个单元,每个单元都可写出上面的质量平衡方程,n个方程联立,化为矩阵形式,得整个湖泊以矩阵方程描述的系统状态微分方程:

式中:C(t)=[c1(t),c2(t),…cn(t)]T,为t时n个单元的污染浓度列向量;U(t)=[u1(t),u2(t)…un(t)]T,为t时n个单元污染物质输入的负荷强度列向量,W(t)=[w1(t),w2(t)…wn(t)]T,为t时n个单元的浓度模型计算随机误差列向量,称系统噪声列向量;A=[ɑi,j]为单元污染物浓度状态转移矩阵;B=[bi,j],输入变量分配矩阵;G=[gi,j]为系统噪声分配矩阵,其中i=1,2,…,n,j=1,2,…,n。
1.2 观测方程
t时对系统中设有监测点的单元观测的水质浓度值,以向量Z(t)表示,它与系统状态变量C(t)的关系可表达成如下的观测方程:

式中:Z(t)=[z1(t),z2(t)…zL(t)]T,为t时L个单元观测的浓度列向量;H(t)=[hi,j(t)],i=1,2,…L,j=1,2,…n,为水质浓度观测分配矩阵,由观测点在系统单元中的位置确定;M(t)=[m1(t),m2(t)…mL(t)]T,为观测随机误差列向量,称观测噪声向量。
式(5)、式(6)是建立的浅水湖泊水质系统滤波模型,应用滤波方法求解,从而做出水质实时过程预报。