灰色系统理论是在20世纪80年代初被提出并发展推广应用的。实际上，灰色系统是黑箱理论的一种推广。黑箱理论认为系统内部结构、参数、特征等一无所知，只能从系统外部的输入、输出信息来研究。如果一个系统的内部特征全部确知，那么这个系统就是一个白箱，我们可以根据系统内部的信息建立各种模型进行研究。灰色系统就是介于“黑”、“白”之间的，或者说灰色系统是指部分信息已知，部分信息未知的系统。
灰色系统作为系统科学的一个分支，正处于发展时期。灰色系统理论包括了关联分析、灰色模型、灰色预测决策、灰色优化控制等等。这里就灰色理论在水文地质学中较常应用的几个方面着重进行介绍。
一、关联分析
关联分析事实上是动态过程发展态势的量化比较分析。关联分析的实质是几种曲线（代表着相应的变量）间的几何形状的比较分析，即认为几何形状越接近，则发展变化的态势越接近，关联程度就越大。这种分析对数据量没有太高的要求。利用作图的直观方法，也可看出几种曲线的形状变化，但只能是定性的，不能够定量化。
作关联分析时，先要指定参考的数据列x0，被比较的数据列常记为：x1，x2，…，xn
现代水文地质学
关联系数的计算公式为：
现代水文地质学
式中：i为比较数列的个数；k为时段；ξi（k）是第k个时刻比较曲线xi与参考曲线x0的相对差值，称这一差值为xi对x0在k时刻的关联系数。上式中ζ是分辨系数，一般取值范围为（0，1）。这样一来，就可以分别求得i个比较曲线与参考曲线的关联系数。在实际应用中，在计算关联系数之前，当各数据列所代表的变量的量纲不同时要进行无量纲化处理，而且要求所有数列有公共交点。为了解决这一问题，应对各数列进行初值化处理。具体作法是用每一个数列的第一个数去除其他数。
由于关联系数很多，信息过于分散，不便于比较，为此有必要将各个时刻关联系数集中为一个值，故引入关联度ri：
现代水文地质学
这样，通过各比较曲线与参考曲线关联度的比较，就可定量地分析变量间的关联程度。
二、优势分析
若参考序列不止一个，比较序列也不止一个，可分别计算所有序列间的关联度，构造关联度矩阵，通过对矩阵中元素大小的分析，找出优势因素和劣势因素。
设有n个参考序列，m个比较序列：
现代水文地质学
则每一个母因素对每个子因素都有关联度，用rij表示与xj的关联度，则有：
现代水文地质学
R中每一行表示同一母因素对不同子因素的影响，每一列表示不同母因素对同一子因素的影响。如果R中第j列的各个子元素大于其他各列的相应子元素，则称子元素xj相对于其他因素为最优，即xj为优势子因素；若R中第i行的各元素均大于其他各行的相应子元素，则称母元素为优势因素。
三、生成数
在水文地质系统的研究中，不少因素都存在着随机性，如降水补给、地下水的污染等等。处理随机因素一般用数理统计方法寻找概率分布或统计规律。而灰色理论则是用数据处理的方法来找数据间的规律。一般称某种数据的处理方式为生成方式。
灰色系统理论常用的数据生成方式有3类：
（一）累加生成（AGO）
原始数据进行累加得到的新数据与数列累加后的数列称为生成数列。记x0为原始数列：
现代水文地质学
生成数列为x（1）：
现代水文地质学
若x（1）（k）=x（0）（i）则称为一次累加生成，记为1-AGO。当α=1时，称为一般累加生成；当α是正整数时，但1＜α≤k，称为去首累加生成。
（二）累减生成（IAGO）
原始数据列中前后两个数据相减，累减是累加的逆运算。累减生成可将累加生成还原为非生成数列。
现代水文地质学
（三）映射生成
除累加与累减生成以外的其他生成方法。
四、灰色（GM）模型（GM（1，1））
首先对数据进行生成处理，使原有数据的随机性弱化。经累加生成的序列呈逐渐递增的形式，灰色理论建立模型就是要构造一个函数去逼近累加生成序列，从而达到利用原始数据序列进行预测未来的目的（邓聚龙，1987年）。
对于非负原始数据作累加生成后将呈现指数规律，因此可以考虑用微分方程的解去逼近生成数列，灰色模型的建立正是基于这一思想。
对于非负数据列：
现代水文地质学
作一次累加得：
现代水文地质学
其中：x（1）（k）=x（0）（i）。
设x（1）（t）满足一阶常微分方程
现代水文地质学
其中a，b为常系数，式（15-30）的解为：
现代水文地质学
根据式（15-31）构造迭代格式，得到预测模型为：
现代水文地质学
灰色建模的途径是依式（15-29）的序列值，用最小二乘法来计算a，b。最后求得
现代水文地质学
其中：
现代水文地质学
利用式（15-33）求得a，b值代入式（15-32）即可进行模型预报。具体在水文地质工作中，利用GM模型可以对地下水水位、降水量、水质状况等进行预报。此外，还可以用来进行灾变预测。
#灰色预测模型GM(1,1)#用法：#假设数列1 2 3 4 5.5 6 7.5 为已知数据，你要预测后面3项，gm11([1 2 3 4 5.5 6 7.5],10) # 10=7+3# 序列输入格式为：x<-c(1,2,3,4,5.5,6,7.5)gm11<-function(x,k){#x为行向量数据#做一次累加n<-length(x)x1<-numeric(n)for(i in 1:n){x1[i]<-sum(x[1:i])}#x1的均值数列z1<-numeric(n)m<-n-1for(j in 1:m){z1[j+1]<-(0.5*x1[j+1]+0.5*x1[j])}Yn=t(t(x[2:n]))B<-matrix(1,nrow=n-1,ncol=2)B[,1]<-t(t(-z1[2:n]))#solve(M)求M的逆#最小二乘法求解参数列u<-solve(t(B)%*%B)%*%t(B)%*%Yna<-u[1]b<-u[2]#预测x2<-numeric(k)x2[1]<-x[1]for(i in 1:k-1){x2[1+i]=(x[1]-b/a)*exp(-a*i)+b/a}x2=c(0,x2)#还原数据y=diff(x2)y}#调用函数x<-c(1,2,3,4,5.5,6,7.5)gm11(x,10)
