灰色聚类分析
灰色聚类是灰色系统理论中的重要分支，通过灰色关联度或白化权函数将观测对象划分为若干可定义类别。该方法特别适用于信息不完全、样本量小、数据贫乏的分类问题，在社会经济系统评价、环境质量评估、工程决策等领域具有广泛应用。
灰色系统理论基础
灰色系统理论由邓聚龙教授于1982年创立，研究对象是“部分信息已知、部分信息未知“的不确定性系统，即灰色系统。
基本概念
灰色系统理论的核心思想在于对“少数据、贫信息“的不确定性系统进行研究。根据系统信息的完备程度，可将系统划分为：
- 白色系统：信息完全明确的系统
- 黑色系统：信息完全未知的系统
- 灰色系统：部分信息已知、部分信息未知的系统
灰数与灰度
灰数是灰色系统理论中的基本元素，记为 \( \otimes \)，表示取值范围已知但确切值未知的数。
设灰数 \( \otimes \in [a, b] \)，其灰度定义为：
\[ g^\circ(\otimes) = \frac{b - a}{b} \]
当 \( g^\circ(\otimes) = 0 \) 时，灰数退化为白数（确定值）；当 \( g^\circ(\otimes) \to \infty \) 时，趋向黑数。
灰色聚类的分类
灰色聚类主要包含两大类方法：
- 灰色关联聚类：基于灰色关联度进行对象或指标的归类
- 灰色白化权函数聚类：基于白化权函数计算聚类系数，进行对象的分类评价
灰色关联聚类
灰色关联聚类以灰色关联分析为基础，通过计算序列间的关联度构建关联矩阵，再按照一定阈值进行聚类划分。
灰色关联度
设参考序列为 \( X_0 = (x_0(1), x_0(2), \ldots, x_0(n)) \)，比较序列为 \( X_i = (x_i(1), x_i(2), \ldots, x_i(n)) \)，则 \( X_i \) 对 \( X_0 \) 在第 \( k \) 个观测点的关联系数为：
\[ \xi_i(k) = \frac{\min_i \min_k |x_0(k) - x_i(k)| + \rho \max_i \max_k |x_0(k) - x_i(k)|}{|x_0(k) - x_i(k)| + \rho \max_i \max_k |x_0(k) - x_i(k)|} \]
其中 \( \rho \in [0, 1] \) 为分辨系数，通常取 \( \rho = 0.5 \)。
关联度为关联系数的均值：
\[ r_i = \frac{1}{n} \sum_{k=1}^{n} \xi_i(k) \]
关联聚类步骤
步骤一：数据无量纲化处理。常用初值化 \( x_i’(k) = x_i(k) / x_i(1) \) 或均值化 \( x_i’(k) = x_i(k) / \bar{x}_i \)。
步骤二：对 \( m \) 个序列两两计算关联度，得到关联矩阵 \( R = (r_{ij})_{m \times m} \)。
步骤三：选取临界值 \( r^* \)，若 \( r_{ij} \geq r^* \)，则认为序列 \( i \) 与序列 \( j \) 属于同一类。
关联聚类的性质
关联矩阵 \( R \) 满足对称性 \( r_{ij} = r_{ji} \) 和规范性 \( 0 < r_{ij} \leq 1 \)，\( r_{ii} = 1 \)。
灰色白化权函数聚类
白化权函数聚类是灰色聚类中应用最广泛的方法，其核心思想是对各聚类指标构造相应的白化权函数，通过计算聚类系数确定对象所属类别。
基本设定
设有 \( m \) 个聚类对象，\( n \) 个聚类指标，\( s \) 个灰类（等级）。记第 \( i \) 个对象关于第 \( j \) 个指标的观测值为 \( x_{ij} \)（\( i = 1, 2, \ldots, m \)；\( j = 1, 2, \ldots, n \)）。
白化权函数的构造
对于第 \( j \) 个指标的第 \( k \) 个灰类，白化权函数 \( f_j^k(x) \) 通常采用以下三种典型形式：
（1）上限测度白化权函数（第一类）
\[ f_j^1(x) = \begin{cases} x / \lambda_j^1 & x \in [0, \lambda_j^1] \\ 1 & x \in [\lambda_j^1, +\infty) \\ \end{cases} \]
（2）适中测度白化权函数（中间类）
\[ f_j^k(x) = \begin{cases} x / \lambda_j^k & x \in [0, \lambda_j^k] \\ (\lambda_j^{k+1} - x) / (\lambda_j^{k+1} - \lambda_j^k) & x \in [\lambda_j^k, \lambda_j^{k+1}] \\ 0 & x \notin [0, \lambda_j^{k+1}] \end{cases} \]
（3）下限测度白化权函数（末类）
\[ f_j^s(x) = \begin{cases} 1 & x \in [0, \lambda_j^{s-1}] \\ (\lambda_j^s - x) / (\lambda_j^s - \lambda_j^{s-1}) & x \in [\lambda_j^{s-1}, \lambda_j^s] \\ 0 & x > \lambda_j^s \end{cases} \]
其中 \( \lambda_j^k \) 为第 \( j \) 个指标第 \( k \) 个灰类的阈值（转折点）。
聚类权重
设第 \( j \) 个指标对第 \( k \) 个灰类的权重为 \( \eta_j^k \)，通常采用如下方式确定：
\[ \eta_j^k = \frac{\lambda_j^k}{\sum_{j=1}^{n} \lambda_j^k} \]
聚类系数
第 \( i \) 个对象关于第 \( k \) 个灰类的聚类系数为：
\[ \sigma_i^k = \sum_{j=1}^{n} f_j^k(x_{ij}) \cdot \eta_j^k \]
聚类判定
对第 \( i \) 个对象，若满足：
\[ \sigma_i^{k^*} = \max_{1 \leq k \leq s} {\sigma_i^k} \]
则判定对象 \( i \) 属于第 \( k^* \) 个灰类。
实际案例分析
以下通过一个完整的环境质量评价案例，演示灰色白化权函数聚类的全部计算过程。
问题描述
某地区对4个监测站点的空气质量进行评价，选取3个污染指标：PM2.5浓度（\(\mu g/m^3\)）、SO\(_2\)浓度（\(\mu g/m^3\)）、NO\(_2\)浓度（\(\mu g/m^3\)），将空气质量分为3个等级：优（灰类1）、良（灰类2）、差（灰类3）。
原始数据
| 站点 | PM2.5 (\(x_1\)) | SO\(_2\) (\(x_2\)) | NO\(_2\) (\(x_3\)) |
|---|---|---|---|
| A | 30 | 40 | 35 |
| B | 55 | 70 | 60 |
| C | 80 | 50 | 90 |
| D | 45 | 85 | 50 |
确定灰类阈值
根据国家标准和专家经验，确定各指标对应各灰类的阈值：
| 指标 | 优 (\(\lambda^1\)) | 良 (\(\lambda^2\)) | 差 (\(\lambda^3\)) |
|---|---|---|---|
| PM2.5 | 35 | 75 | 115 |
| SO\(_2\) | 50 | 150 | 250 |
| NO\(_2\) | 40 | 80 | 120 |
步骤一：构造白化权函数
以PM2.5指标为例（\( j = 1 \)）：
灰类1（优）的白化权函数：
\[ f_1^1(x) = \begin{cases} 1 & x \in [0, 35] \\ (75 - x) / (75 - 35) & x \in (35, 75] \\ 0 & x > 75 \end{cases} \]
灰类2（良）的白化权函数：
\[ f_1^2(x) = \begin{cases} x / 75 & x \in [0, 75] \\ (115 - x) / (115 - 75) & x \in (75, 115] \\ 0 & x > 115 \end{cases} \]
灰类3（差）的白化权函数：
\[ f_1^3(x) = \begin{cases} x / 115 & x \in [0, 115] \\ 1 & x > 115 \end{cases} \]
同理可构造SO\(_2\)和NO\(_2\)的白化权函数。
步骤二：计算白化权函数值
以站点A为例演示详细计算，其余站点同理：
站点A（PM2.5=30, SO\(_2\)=40, NO\(_2\)=35）：
- \( f_1^1(30) = 1 \)，\( f_1^2(30) = 30/75 = 0.400 \)，\( f_1^3(30) = 30/115 = 0.261 \)
- \( f_2^1(40) = 1 \)，\( f_2^2(40) = 40/150 = 0.267 \)，\( f_2^3(40) = 40/250 = 0.160 \)
- \( f_3^1(35) = 1 \)，\( f_3^2(35) = 35/80 = 0.438 \)，\( f_3^3(35) = 35/120 = 0.292 \)
站点B（PM2.5=55, SO\(_2\)=70, NO\(_2\)=60）：
- \( f_1^1(55) = (75-55)/40 = 0.500 \)，\( f_1^2(55) = 55/75 = 0.733 \)，\( f_1^3(55) = 55/115 = 0.478 \)
- \( f_2^1(70) = (150-70)/100 = 0.800 \)，\( f_2^2(70) = 70/150 = 0.467 \)，\( f_2^3(70) = 70/250 = 0.280 \)
- \( f_3^1(60) = (80-60)/40 = 0.500 \)，\( f_3^2(60) = 60/80 = 0.750 \)，\( f_3^3(60) = 60/120 = 0.500 \)
站点C（PM2.5=80, SO\(_2\)=50, NO\(_2\)=90）：
- \( f_1^1(80) = 0 \)，\( f_1^2(80) = (115-80)/40 = 0.875 \)，\( f_1^3(80) = 80/115 = 0.696 \)
- \( f_2^1(50) = 1 \)，\( f_2^2(50) = 50/150 = 0.333 \)，\( f_2^3(50) = 50/250 = 0.200 \)
- \( f_3^1(90) = 0 \)，\( f_3^2(90) = (120-90)/40 = 0.750 \)，\( f_3^3(90) = 90/120 = 0.750 \)
站点D（PM2.5=45, SO\(_2\)=85, NO\(_2\)=50）：
- \( f_1^1(45) = (75-45)/40 = 0.750 \)，\( f_1^2(45) = 45/75 = 0.600 \)，\( f_1^3(45) = 45/115 = 0.391 \)
- \( f_2^1(85) = (150-85)/100 = 0.650 \)，\( f_2^2(85) = 85/150 = 0.567 \)，\( f_2^3(85) = 85/250 = 0.340 \)
- \( f_3^1(50) = (80-50)/40 = 0.750 \)，\( f_3^2(50) = 50/80 = 0.625 \)，\( f_3^3(50) = 50/120 = 0.417 \)
步骤三：计算聚类权重
灰类1（优）：\( \eta_1^1 = 35/125 = 0.280 \)，\( \eta_2^1 = 50/125 = 0.400 \)，\( \eta_3^1 = 40/125 = 0.320 \)
灰类2（良）：\( \eta_1^2 = 75/305 = 0.246 \)，\( \eta_2^2 = 150/305 = 0.492 \)，\( \eta_3^2 = 80/305 = 0.262 \)
灰类3（差）：\( \eta_1^3 = 115/485 = 0.237 \)，\( \eta_2^3 = 250/485 = 0.515 \)，\( \eta_3^3 = 120/485 = 0.247 \)
步骤四：计算聚类系数
站点A：\( \sigma_A^1 = 1 \times 0.280 + 1 \times 0.400 + 1 \times 0.320 = 1.000 \) \( \sigma_A^2 = 0.400 \times 0.246 + 0.267 \times 0.492 + 0.438 \times 0.262 = 0.345 \) \( \sigma_A^3 = 0.261 \times 0.237 + 0.160 \times 0.515 + 0.292 \times 0.247 = 0.216 \)
站点B：\( \sigma_B^1 = 0.500 \times 0.280 + 0.800 \times 0.400 + 0.500 \times 0.320 = 0.620 \) \( \sigma_B^2 = 0.733 \times 0.246 + 0.467 \times 0.492 + 0.750 \times 0.262 = 0.607 \) \( \sigma_B^3 = 0.478 \times 0.237 + 0.280 \times 0.515 + 0.500 \times 0.247 = 0.381 \)
站点C：\( \sigma_C^1 = 0 \times 0.280 + 1 \times 0.400 + 0 \times 0.320 = 0.400 \) \( \sigma_C^2 = 0.875 \times 0.246 + 0.333 \times 0.492 + 0.750 \times 0.262 = 0.576 \) \( \sigma_C^3 = 0.696 \times 0.237 + 0.200 \times 0.515 + 0.750 \times 0.247 = 0.453 \)
站点D：\( \sigma_D^1 = 0.750 \times 0.280 + 0.650 \times 0.400 + 0.750 \times 0.320 = 0.710 \) \( \sigma_D^2 = 0.600 \times 0.246 + 0.567 \times 0.492 + 0.625 \times 0.262 = 0.590 \) \( \sigma_D^3 = 0.391 \times 0.237 + 0.340 \times 0.515 + 0.417 \times 0.247 = 0.371 \)
步骤五：聚类判定
| 站点 | \(\sigma^1\) | \(\sigma^2\) | \(\sigma^3\) | 判定结果 |
|---|---|---|---|---|
| A | 1.000 | 0.345 | 0.216 | 优 |
| B | 0.620 | 0.607 | 0.381 | 优 |
| C | 0.400 | 0.576 | 0.453 | 良 |
| D | 0.710 | 0.590 | 0.371 | 优 |
结论：站点A、B、D的空气质量属于“优“等级，站点C的空气质量属于“良“等级。其中站点B的聚类系数 \(\sigma^1\) 与 \(\sigma^2\) 十分接近，说明其空气质量处于优良交界状态，需要重点关注。
Python代码实现
以下代码实现了灰色白化权函数聚类的完整流程，可直接用于实际数据分析。
import numpy as np
def whitening_type1(x, lam1, lam2):
"""上限测度白化权函数（第一灰类）"""
if x <= lam1:
return 1.0
elif x <= lam2:
return (lam2 - x) / (lam2 - lam1)
return 0.0
def whitening_mid(x, lam_k, lam_k1):
"""适中测度白化权函数（中间灰类）"""
if x <= lam_k:
return x / lam_k
elif x <= lam_k1:
return (lam_k1 - x) / (lam_k1 - lam_k)
return 0.0
def whitening_last(x, lam_s):
"""下限测度白化权函数（末灰类）"""
return x / lam_s if x <= lam_s else 1.0
class GreyClusteringModel:
"""灰色白化权函数聚类模型"""
def __init__(self, data, thresholds):
"""
参数:
data: shape (m, n)，m个对象，n个指标
thresholds: shape (n, s)，n个指标，s个灰类的阈值
"""
self.data = np.array(data, dtype=float)
self.thresholds = np.array(thresholds, dtype=float)
self.m, self.n = self.data.shape
self.s = self.thresholds.shape[1]
self.weights = None
self.cluster_coefficients = None
self.results = None
def compute_weights(self):
"""计算各灰类的聚类权重"""
self.weights = np.zeros((self.n, self.s))
for k in range(self.s):
col_sum = np.sum(self.thresholds[:, k])
self.weights[:, k] = self.thresholds[:, k] / col_sum
return self.weights
def _whitening_value(self, x, j, k):
"""计算第j个指标第k个灰类的白化权函数值"""
if k == 0:
return whitening_type1(x, self.thresholds[j, 0],
self.thresholds[j, 1])
elif k == self.s - 1:
return whitening_last(x, self.thresholds[j, k])
else:
return whitening_mid(x, self.thresholds[j, k],
self.thresholds[j, k + 1])
def compute_cluster_coefficients(self):
"""计算聚类系数矩阵"""
if self.weights is None:
self.compute_weights()
self.cluster_coefficients = np.zeros((self.m, self.s))
for i in range(self.m):
for k in range(self.s):
sigma = sum(
self._whitening_value(self.data[i, j], j, k)
* self.weights[j, k]
for j in range(self.n)
)
self.cluster_coefficients[i, k] = sigma
return self.cluster_coefficients
def run(self):
"""运行完整的聚类分析流程"""
self.compute_weights()
self.compute_cluster_coefficients()
self.results = np.argmax(self.cluster_coefficients, axis=1) + 1
return self.results
def print_results(self, class_names=None, object_names=None):
"""打印聚类结果"""
if self.results is None:
self.run()
if class_names is None:
class_names = [f"灰类{k+1}" for k in range(self.s)]
if object_names is None:
object_names = [f"对象{i+1}" for i in range(self.m)]
print("对象\t" + "\t".join(class_names) + "\t判定")
for i in range(self.m):
coeffs = "\t".join(f"{self.cluster_coefficients[i, k]:.4f}"
for k in range(self.s))
print(f"{object_names[i]}\t{coeffs}"
f"\t{class_names[self.results[i] - 1]}")
# ============ 案例运行 ============
if __name__ == "__main__":
data = np.array([
[30, 40, 35], # 站点A
[55, 70, 60], # 站点B
[80, 50, 90], # 站点C
[45, 85, 50], # 站点D
])
thresholds = np.array([
[35, 75, 115], # PM2.5
[50, 150, 250], # SO2
[40, 80, 120], # NO2
])
model = GreyClusteringModel(data, thresholds)
model.run()
model.print_results(
class_names=["优", "良", "差"],
object_names=["站点A", "站点B", "站点C", "站点D"]
)
使用时只需准备好数据矩阵（形状 \(m \times n\)）和阈值矩阵（形状 \(n \times s\)），实例化 GreyClusteringModel 并调用 run() 方法即可获得聚类结果。
应用注意事项与局限性
灰色聚类方法虽然在处理小样本不确定性问题上具有独特优势，但在实际应用中仍需注意诸多问题。
应用注意事项
1. 阈值的确定
白化权函数中的阈值 \( \lambda_j^k \) 对聚类结果影响重大。确定方式包括：
- 参照国家标准或行业规范
- 借助专家经验判断
- 利用历史数据统计分析
不同阈值选取可能导致完全不同的聚类结论，因此应进行敏感性分析。
2. 权重的选取
除了基于阈值比例确定权重外，还可采用：
- 层次分析法（AHP）确定主观权重
- 熵权法确定客观权重
- 组合赋权法综合确定权重
权重的合理性直接决定了聚类结果的可靠性。
3. 灰类数目的确定
灰类数目 \( s \) 通常根据实际问题背景确定（如优良中差四级），但也应注意：
- 灰类过少则分辨力不足
- 灰类过多则各类边界模糊，增加判定难度
- 一般建议3至5个灰类为宜
4. 数据预处理
- 指标方向一致化：确保所有指标具有相同的优劣方向
- 异常值处理：极端值可能严重影响白化权函数的取值
- 缺失值处理：灰色系统虽然容许一定程度的信息缺失，但大量缺失仍需插补
5. 白化权函数形式的选择
除标准三角形外，也可采用梯形、高斯型或端点修正型白化权函数，根据具体问题灵活选择。
局限性分析
- 主观性较强：阈值设定、灰类划分、权重确定等环节均依赖主观判断
- 对阈值敏感：阈值的微小变化可能导致边界对象的分类发生改变
- 难以处理高维数据：指标数量较多时，权重分配和白化权函数构造变得复杂
- 缺乏统计检验：与统计聚类方法不同，灰色聚类缺乏严格的显著性检验手段
- 硬划分问题：基于最大聚类系数进行硬判定，对聚类系数接近的情况存在误判风险
与其他方法的比较
| 特征 | 灰色聚类 | 模糊聚类 | K-means |
|---|---|---|---|
| 样本需求 | 小样本适用 | 中等样本 | 大样本 |
| 先验信息 | 需要阈值 | 需要隶属函数 | 需要K值 |
| 不确定性处理 | 灰色不确定性 | 模糊不确定性 | 较弱 |
| 可解释性 | 强 | 较强 | 中等 |
| 适用场景 | 评价分级 | 软分类 | 探索性分析 |
改进方向
- 组合赋权：将主观权重与客观权重结合，降低主观因素影响
- 区间灰数扩展：将确定观测值推广到区间灰数，更好地描述数据不确定性
- 与其他方法融合：如灰色-模糊聚类、灰色-TOPSIS组合评价等
- 动态灰色聚类：考虑时间序列变化的动态评价问题
参考文献
- 邓聚龙. 灰色系统理论教程. 华中理工大学出版社, 1990.
- 刘思峰等. 灰色系统理论及其应用(第七版). 科学出版社, 2014.
- Liu S F, Lin Y. Grey Systems: Theory and Applications. Springer, 2010.