1. 研究背景与意义
青海湖作为青藏高原最大的湖泊,是我国西部生态安全屏障和水资源调配的核心组成部分[1]。流域内水文过程长期变化的定量化模拟对于研究生态环境演化具有十分重要的科学意义。特别地,流域径流量的准确评估在灌溉调度、水资源管理系统的设计和实施、生态环境保护、洪水预警和防治等领域中发挥着关键作用[2]-[5]。然而,径流受多种相互制约的因素控制,如降水、融雪、气温、土地利用、河流形态、湖面面积和集水区特征等[6] [7],因此径流的准确模拟和预测一直是水文学领域的挑战性问题。此外,径流变化规律受气候、流域自然地理条件以及人类活动等多元因素影响,呈现出复杂多变的特性,包括多时间尺度、突变性和非线性等特点。
近年来,水文过程模拟方法通常包括基于过程的物理模型,如SWAT、HBV、VIC、新安江模型等[8],以及基于机器学习/深度学习的数据驱动模型,如SVM、LSTM、GRU、Transformer等[9]。水文过程驱动模型通常需要多种输入数据,如气象、土壤、地形、植被、土地利用等。水文过程模型能够对预测结果的成因做出物理解释,然而模型的简化程度、数据误差以及对流域水文过程概念性模型的偏差带来了水文模型的误差和不确定性。数据驱动模型不考虑这些复杂因素与水文过程之间直接的物理关系和影响机制,旨在建立最佳的输入和输出之间的数学关系,避免了模型参数和模型的不确定性带来的建模误差,在水文过程模拟中逐渐展现出较强的应用潜力。然而,这类方法通常被称为黑箱模型,难以清晰展示模型输入对输出的约束以及两者的物理关系。因此,耦合物理模型和数据驱动模型,整合并利用两者各自的优势,为水文过程模拟和预测带来了新的思路。特别是随着深度学习方法快速发展,建立耦合径流预测模型已成为一种广泛的研究趋势[10] [11],特别是对高寒缺资料地区具有更大的应用潜力[12] [13]。
基于上述背景,本文分析了青海湖流域气象水文要素的演变规律,构建了基于水文模型FLEX和深度学习算法GRU的混合模型,模拟预测了2000~2011年青海湖流域逐日径流演变规律,并利用SHAP可解释性算法揭示了影响径流变化的关键驱动因子。FLEX (FluxExchange)水文模型是集总式概念性水文模型[12],由主衰退曲线推导出蓄存–排泄关系,克服了通常对低流量的忽略,能够高性能地模拟峰值流量。这一模型适合于干旱、寒冷区域,展示了模型框架在刻画流域多样性多尺度性方面的潜力[14]。门控循环单元(GRU)隶属递归神经网络(Recurrent Neural Network, RNN),是长短期记忆(Long Short-Term Memory, LSTM)网络的一种有效变种[15]。GRU在表达上更优化和简洁,擅长捕捉长时间序列数据特征,并已成功应用于水文预报领域[16]。
同时,为了增强FLEX-GRU混合模型的性能,本文采用了三种策略:1) 应用差分进化自适应算法(DiffeRential Evolution Adaptive Metropolis algorithm, DREAM(zs)) [17]反演FLEX水文模型参数从而优化模型;2) 在深度神经网络训练模型之前采用变分模态分解(Variational Model Decomposition, VMD) [18]提取径流序列数据的隐含信息和特征;3) 利用麻雀搜索算法(Sparrow Search Algorithm, SSA) [19]对深度学习GRU的超参数进行优化。此外应用Shapley Additive Explanations (SHAP)方法[20]分析输入特征对输出的贡献从而揭示气象要素对径流的影响及重要度。
本文以青海湖流域的第一大支流布哈河为研究区域,基于2000~2011年布哈河口水文站的逐日气象和径流观测数据,首先采用DREAM(zs)方法优化FLEX水文模型对多年径流进行模拟,然后通过变分模态分解VMD对FLEX水文模型得到的径流数据预处理,并将其与构建的SSA-GRU多维神经网络模型进行结合进行模型训练,构建组合模型FLEX-VMD-SSA-GRU,最后利用SHAP结果分析影响日径流的关键气象因子。模拟结果显示,将水文模型和深度学习方法相结合,并引入更多的气象因子作为输入,显著提高了对径流峰值的预测精度,为高海拔寒旱区流域水文过程的模拟和预测提供了一种新方法,并阐释了影响水文过程的关键气象要素。
2. 研究区概况及数据收集
2.1. 研究区概况
研究区为我国青藏高原青海湖流域内的第一大支流布哈河,如图1所示。布哈河流域地处青海湖流域西北部,是青海湖流域内最大的子流域,研究区域主要涵盖天峻县及流域周边的刚察县和共和县。布哈河源头为疏勒南山,自西北向东南蜿蜒流动,最终汇入青海湖。流域内海拔变化显著,跨越3185 m到5229 m,平均海拔为4172 m,4000 m至4500 m之间的海拔区间面积占据整个流域的47.69%。布哈河流域的集水面积近15,000 km2,流域宽度为50 km,河流长度约为286 km,多年平均流量为7.85 × 108 m3,占据青海湖入湖总量的50%~60%。
布哈河流域多年平均降水量在300~500 mm之间,潜在蒸发量在1000~2000 mm之间,地区分布不均,年际变化较小。西格尔区上游分布有少量冰川,分布在布哈河支流中的冰川总面积约13.29 km2,每年消融水量在0.07 × 108 m3左右,占布哈河每年径流量的近0.9%。布哈河年内径流分布极不均匀,在汛期(6月至9月)占全年的82.7%,而在枯水期(12月至次年3月)仅占全年的3.24%。
2.2. 数据来源及预处理
径流数据选取布哈河口水文站2000~2011年的日径流观测数据。气象资料来源于中国气象数据网和青海省气象局。采用距离权重反比法对布哈河流域周边的托勒、野牛沟、张掖、祁连、山丹、德令哈、刚察、门源、都兰、共和、西宁、贵德12个气象站点的观测数据插值,得到布哈河流域的气象数据。模型构建中,2000~2001年为FLEX模型预热期,2002~2008年为FLEX模型的率定期(深度学习模型的训练期),2009~2011年为FLEX模型的验证期(深度学习模型的测试期)。
图1. 研究区概况图:左图为青海湖流域位置图,右图为布哈河流域地形图(根据高程数据用ArcGiss所画)
3. 主要研究方法
3.1. FLEX-DREAM(ZS)模型
FLEX (FluxExchange)模型是一种集总式概念性水文模型,通过模拟流域内主要水文过程实现径流估算。本文采用FLEX建模框架构建日尺度降雨–径流模型,输入为日降水量与潜在蒸发量,输出为逐日径流量。模型结构由三部分组成:① 截留层,以截留水箱模拟冠层截留与蒸发;② 非饱和层,表征土壤蓄水能力及水分垂直运移;③ 径流层,由快速响应水箱(产生地表径流与壤中流)和慢速响应水箱(产生地下径流)耦合而成。产流机制以饱和过剩产流为主,径流量被参数化为流域相对蓄水量的非线性函数,弱化对降雨强度的瞬时依赖。图2左侧展示了FLEX模型的基本框架和计算流程。图中P表示降雨,E表示蒸发,Pe表示有效降雨,Ea表示实际蒸发,Su表示排泄至非饱和土壤层的径流,Sf表示排泄至快速响应水库的径流,Ss表示表示排泄至慢速响应水库的径流,Q表示总径流,Qf表示快速径流;Qs表示慢速径流。
FLEX模型的主要参数(表1中所示)由DREAM(ZS)算法进行反演。DREAM(ZS) (DiffeRential Evolution Adaptive Metropolis algorithm)是一种改进的多链马尔可夫链蒙特卡罗(Markov Chain Monte Carlo, MCMC)算法,结合了差分进化、子空间探索、过去状态采样以及多重尝试的Metropolis采样策略。该算法通过多条交互式马尔可夫链并行探索参数空间法,从而动态调整采样策略,更高效地收敛到目标分布,展现出强大的收敛能力和鲁棒性,常被用于水文模型中的贝叶斯推断和参数估计。
3.2. VMD-SSA-GRU深度神经网络
研究表明,模态分解在处理非线性、非平稳信号时表现出色,能有效提取信号中的本征模态函数,从而去除噪声、分离不同来源的振动成分。变分模态分解(VMD) [18]是一种将信号分解过程转化为变分优化问题的创新方法,将信号分解转化为变分优化问题,是一个完全非递归的变分模式分解模型,一次性提取K个本征模态,同时估计其中心频率。VMD在其非递归框架内置维纳滤波,它表现出强大的抗干扰能力,可抑制高斯白噪声并避免模态混叠;通过约束各模式在基带解调后的平滑性,实现信号–噪声自适应分离,具有鲁棒性强、无混叠、可并行等优点。VMD与GRU网络结合的流程在图2中间示意。
麻雀搜索算法(SSA)是一种受自然启发的元启发式优化算法。该算法是以麻雀的觅食和防捕食行为为模型的群体智能仿生算法[19],模拟了麻雀种群在觅食过程中不同个体间的分工与社会互动。此外,算法中预警机制的引入有效防止了算法陷入局部最优。因为算法的结构简单、收敛快,已被广泛用于时间序列预测、神经网络参数优化等任务。图2右侧展示了SSA的主要步骤。
门控循环单元(GRU)通常被视为更复杂的LSTM网络的简化版本。它专门设计用于处理序列数据并解决基础RNN的常见问题。GRU通过合并LSTM中的遗忘门和输入门,保留原有的重置门,从而简化了模型结构。这一简化使得GRU在计算时减少了矩阵乘法的计算,因此在大规模样本时,能够在保持预测效果的同时,显著降低计算时间。GRU的基本结构展示在图2中间。
3.3. FLEX-VMD-SSA-GRU混合模型
基于上述方法将FLEX模型的模拟结果与其相关数据一起作为深度学习的输入进行模型训练,构建了FLEX-VMD-SSA-GRU混合模型,结构如图2所示。
图2. FLEX-VMD-SSA-GRU混合模型结构示意图
3.4. SHAP可解释性方法
SHAP是以合作博弈论为基石的机器学习结果解释模型,通过计算Shapley值精确量化其对单次预测结果的边际贡献。与传统“全局重要性”仅能排序不同,SHAP同时给出方向与强度。SHAP值的最大优势是能反映每个变量特性的影响,并且还显示出正负影响。SHAP为正,则贡献为正影响推高预测,否则贡献为负影响抑制预测,从而揭示影响模型结果的原因和程度。SHAP的绝对值越大,说明特征i的贡献越大。其计算公式为:
(1)
式中:φi为特征在样本中的SHAP值,S是特征的子集;F是所有特征子集的集合;
表示集合S中输入特征的值;
是用特征i训练模型f。
3.5. 模型评价指标
为了评估所构建的混合在径流模拟中的性能,选用了纳什效率系数NSE、平均绝对误差MAE、平均相对误差MAPE、均方根误差RMSE,以及决定系数R2等五种评价指标。这些指标从不同角度反映了模型的预测能力、稳定性和准确性。其中,NSE能够综合评估模型的预测效能,而MAPE和R2则用以衡量模型的稳定程度;MAE和RMSE则侧重于反映模型的预测准确程度。它们的表达式如下:
(2)
(3)
(4)
(5)
(6)
上式中,
表示径流模拟值,
表示径流实测值,
表示实测径流的平均值,n为径流序列样本个数。
4. 结果与分析
4.1. 基于FLEX-DREAM(ZS)模型的径流模拟
首先结合FLEX水文模型和参数优化方法DREAM(ZS)对2000~2011年布哈河径流量进行模拟。输入数据为日降水和日蒸发数据。为了削弱初始状态的不确定性对模型预测结果的影响,将2000年1月1日~2001年12月31日作为模型的预热期,2002年1月1日~2011年12月31日的数据对模型进行率定和验证。其中,2002年1月1日~2008年12月31日为率定期,2009年1月1日~2011年12月31日为验证期。
采用DREAM(ZS)反演FLEX水文模型的七个水文参数,包括最大截留量Imax,土壤蓄水量Smax,最大下渗率Qmax,蒸发量参数αE,径流参数αF,快速水库响应时间常数KF和慢速水库响应时间常数KS,率定结果如表1所示。
基于反演得到的水文模型参数,模拟和预测布哈河的径流变化。图3中展示了FLEX模型模拟的日径流和实际观测值的对比。可以看到,无论是在率定期还是验证期,模拟值与实测值的变化趋势大致相同,同时模型总体上能够捕捉到径流量的峰值,特别是对低谷值有很好的刻画。然而FLEX模型在高峰值模拟方面仍然面临挑战,特别是在2002、2005、2006和2008年。因此尝试将FLEX模型和深度学习方法相结合,旨在提高径流模拟的精度。
表1. FLEX模型参数的率定结果
参数 |
单位 |
范围 |
率定值 |
最大截留量(Imax) |
mm |
0~1 |
0.292 |
土壤蓄水量(Smax) |
mm |
1~5000 |
4593.396 |
最大下渗率(Qsmax) |
mm/d |
0~100 |
75.758 |
蒸发系数(αE) |
- |
0~200 |
169.696 |
径流系数(αF) |
- |
−100~100 |
47.996 |
快速水库响应时间(KF) |
days |
0~100 |
1.282 |
慢速水库响应时间(KS) |
days |
0~200 |
5.783 |
图3. FLEX模型模拟的日径流和实际观测值的对比
4.2. 基于FLEX-VMD-SSA-GRU的径流模拟
4.2.1. VMD分解结果
在利用GRU训练径流模型之前,首先采用VMD分解算法对FLEX模型模拟得到的径流序列进行预处理,提取径流序列的主要特征。VMD算法核心参数为模态分量个数K,由图4所示,针对布哈河径流数据的K值为13。
图4. 基于VMD的布哈河流域径流数据的分解图
4.2.2. VMD-SSA-GRU模型参数及输入设置
本文采用麻雀搜索算法SSA优化深度神经网络GRU模型的超参数。SSA的种群规模为100,最大迭代次数为50,20%为发现者,其余为追随者,警告值为0.8。设置GRU层神经元个数N1、N2和每批次样本数量的取值范围分别为[10, 200],[10, 200]和[100, 500]。优化得到的GRU的隐藏层个数为2,隐含层节点为10,初始学习率为0.001。
为了筛选出合理有效的气象因子作为混合模型的输入,本研究采用通过皮尔逊、肯德尔、斯皮尔曼三种相关性分析方法评估了各气象因子与径流的相关性,确定了相关性较高的日降水P、日蒸发E、基于日温度的最大、最小和平均温度(Tmax, Tmin, Tmean)、相对湿度RH和地表净辐射Rn七个气象要素。各气象变量与模拟日径流的相关性如表2所示。这七个气象要素将与FLEX模型模拟的日径流通过设置不同的组合方式一同作为深度神经网络模型的输入。与FLEX模型保持一致,2002年1月1日~2008年12月31日为模型训练期(占全部数据集的70%),2009年1月1日~2011年12月31日为模型测试期(占全部数据集的30%)。
表2. 各气象因子与日径流的相关性
相关性 分析方法 |
日蒸发/E |
日降水/P |
日均温度/Tmean |
日最高温度/Tmax |
日最低温度/Tmin |
相对湿度/RH |
地表净辐射/Rn |
日照时长/N |
风速/U |
皮尔逊 |
0.45 |
0.37 |
0.57 |
0.57 |
0.65 |
0.56 |
0.42 |
0.03 |
0.12 |
肯德尔 |
0.3 |
0.17 |
0.45. |
0.42 |
0.49 |
0.42 |
0.24 |
0.03 |
0.07 |
斯皮尔曼 |
0.45 |
0.27 |
0.66 |
0.62 |
0.70 |
0.63 |
0.37 |
0.04 |
0.11 |
以这七个气象变量与FLEX模拟结果进行不同组合,共设置14组不同输入的情景。其中情景FLEX为仅以FLEX的输出结果QF作为深度学习的输入;情景1~情景7是在QF的基础上,增加单一气象变量;情景8~情景13是在QF的基础上,依次累加一个气象变量。表3列出了用于实现FLEX模型和神经网络结合的混合模型的输入组合。
表3. 混合模型的输入设置(F-VSG中F代指FLEX,VSG代指VMD-SSA-GRU)
情景 |
模型 |
输入 |
0 |
FLEX |
QF |
1 |
F-VSG1 |
QF, P |
2 |
F-VSG2 |
QF, E |
3 |
F-VSG3 |
QF, RH |
4 |
F-VSG4 |
QF, Rn |
5 |
F-VSG5 |
QF, Tmean |
6 |
F-VSG6 |
QF, Tmax |
7 |
F-VSG7 |
QF, Tmin |
8 |
F-VSG8 |
QF, P, E |
9 |
F-VSG9 |
QF, P, E, Tmean |
10 |
F-VSG10 |
QF, P, E, T (Tmean, Tmax) |
11 |
F-VSG11 |
QF, P, E, T (Tmean, Tmax, Tmin) |
12 |
F-VSG12 |
QF, P, E, T (Tmean, Tmax, Tmin), RH |
13 |
F-VSG13 |
QF, P, E, T (Tmean, Tmax, Tmin), RH, Rn |
4.2.3. FLEX-VMD-SSA-GRU模型的径流模拟结果
14种情景下的径流模拟结果显示在图5(a)和图5(b)中,可以看到所有模型能够刻画径流的总体趋势。为了更清晰地对比不同输入的模拟结果,将2010年12月的结果进行了局部放大,如图中所示。F-VSG13相比于其他组合更接近实测值,特别在低径流时期,围绕实测值波动,且波动幅度更小。
为了进一步定量化评估混合模型在不同输入情景下的模拟结果,将上面14个模型模拟的径流与实测数据对比,并统计了评价指标,详见表4。总体来看,随着输入中气象变量的增加,FLEX模拟的日径流精度逐步提高,其中输入因子最多的情景F-VSG13的模拟精度也最高。与单独FLEX水文模型结果相比,F-VSG13得到的测试集对应的NSE和R2分别提高了10%和5%,MAE、MAPE和RMSE分别降低21%、19%和15%。由此可见,水文模型对深度学习的优化策略具有引导作用,同时利用了深度学习对数据特征较强的学习能力,相较于水文模型,混合模型在模拟和预测径流过程中展示出模拟精度的优势。
表4. 混合模型模拟结果的评价指标(F-VSG中F代指FLEX,VSG代指VMD-SSA-GRU)
|
训练集 |
测试集 |
评价指标 |
NSE |
MAE |
MAPE |
RMSE |
R2 |
NSE |
MAE |
MAPE |
RMSE |
R2 |
FLEX |
0.755 |
0.069 |
0.586 |
0.141 |
0.871 |
0.719 |
0.080 |
0.477 |
0.157 |
0.851 |
F-VSG1 |
0.744 |
0.070 |
0.604 |
0.145 |
0.864 |
0.721 |
0.080 |
0.477 |
0.156 |
0.851 |
F-VSG2 |
0.742 |
0.071 |
0.731 |
0.145 |
0.864 |
0.723 |
0.080 |
0.572 |
0.155 |
0.852 |
F-VSG3 |
0.744 |
0.070 |
0.640 |
0.145 |
0.864 |
0.723 |
0.079 |
0.510 |
0.156 |
0.852 |
F-VSG4 |
0.742 |
0.070 |
0.608 |
0.145 |
0.864 |
0.722 |
0.079 |
0.458 |
0.156 |
0.852 |
F-VSG5 |
0.772 |
0.068 |
0.621 |
0.136 |
0.880 |
0.760 |
0.074 |
0.427 |
0.145 |
0.874 |
F-VSG6 |
0.744 |
0.070 |
0.603 |
0.144 |
0.865 |
0.759 |
0.074 |
0.419 |
0.145 |
0.873 |
F-VSG7 |
0.744 |
0.070 |
0.602 |
0.145 |
0.864 |
0.759 |
0.074 |
0.425 |
0.145 |
0.873 |
F-VSG8 |
0.744 |
0.071 |
0.758 |
0.145 |
0.864 |
0.725 |
0.080 |
0.605 |
0.155 |
0.852 |
F-VSG9 |
0.772 |
0.069 |
0.675 |
0.136 |
0.880 |
0.758 |
0.075 |
0.454 |
0.145 |
0.873 |
F-VSG10 |
0.781 |
0.068 |
0.637 |
0.134 |
0.884 |
0.765 |
0.074 |
0.476 |
0.143 |
0.876 |
F-VSG11 |
0.789 |
0.066 |
0.635 |
0.131 |
0.891 |
0.771 |
0.073 |
0.451 |
0.142 |
0.879 |
F-VSG12 |
0.799 |
0.065 |
0.586 |
0.128 |
0.895 |
0.780 |
0.072 |
0.434 |
0.138 |
0.885 |
F-VSG13 |
0.806 |
0.064 |
0.550 |
0.126 |
0.898 |
0.793 |
0.068 |
0.367 |
0.134 |
0.892 |
![]()
![]()
图5. 不同输入条件下得到的日径流模拟值与实测值之间的对比:左图为情景F-VGS1~F-VGS7与情景FLEX的对比,右图为情景F-VGS8~F-VGS13与情景FLEX的对比
为了更加全面评估FLEX-VMD-SSA-GRU模型对大小不同程度径流量的模拟能力,将实测径流数据按照<5%、5%~95%、>95%的标准划分为三个区间:低径流(0~0.01 mm/d)、中径流(0.01~0.8 mm/d)和高径流(>0.8 mm/d),并通过评估指标对比了FLEX模型和F-VSG13模型在不同区间对径流模拟的结果,如图6(a)、图6(b)和图6(c)所示。总体上,对于高中低三个区间的径流动态,F-VSG13的模拟精度都高于FLEX模型。更进一步,相较于FLEX模型,F-VSG13对高径流的模拟精度提升幅度较大,MAE、MAPE和RMSE分别下降50.7%、50.9%和54.7%;中径流模拟精度的提高幅度次之,MAE、MAPE和RMSE分别下降13.6%、21.7%和14.9%;低径流模拟精度提升的幅度最小,MAE、MAPE和RMSE分别下降17%、16.7%和13.6%。这表明混合模型对于刻画高径流更有优势。
4.3. 基于SHAP的可解释性分析
为了解释模型的输入和输出之间的决策关联,提高混合模型的透明度,采用SHAP可解释性方法计算了各输入变量对径流长期趋势的贡献,结果如图7所示。图7(a)为各气象变量对应的SHAP值,按照从大到小排列,图7(b)为各气象变量的SHAP平均值,蓝色表示低值,红色表示高值。从图中可见,气温(包括Tmean、Tmax、Tmin)对径流长期趋势的影响最为显著,且Tmean和Tmin与径流量之间存在较强的正相关关系。其原因可归于布哈河流域处于青藏高原这一高海拔地区,上游分布有部分冰川,因此在解冻期内冰川融雪增加了布哈河径流量。同时,温度引起的季节性冻土融化也会导致径流增加。影响径流的因素其次有蒸发E、降水P和相对湿度RH,其中蒸发主要表现为负相关关系,降水主要表现为正相关关系。地表辐射Rn影响程度最低。
图6. FLEX模型和F-VSG13模型模拟的日径流对应的误差指标直方图:(a) MAE,(b) MAPE,(c) RMSE
图7. 不同气象变量对径流长期趋势的贡献:(a) 各气象变量的SHAP值,(b) 各气象变量的平均贡献值
5. 结论
本研究基于概念性水文模型FLEX和深度学习方法GRU构建了一种径流模拟预测的混合模型FLEX-VMD-SSA-GRU,并将模型应用于我国青海湖流域的第一大支流布哈河。该模型基于优化的FLEX模拟得到的径流结果,并利用变分模态分解VMD对模拟数据进行预处理,然后和GRU神经网络(通过麻雀搜索算法SSA优化网络超参数)相结合得到径流模拟的混合模型。模拟结果显示,混合模型能够有效模拟径流的总体局势和峰值,并且更多的气象输入因子有利于获得更高的模拟精度。此外,采用SHAP可解释性方法评估了日降水量、日潜在蒸发量、日平均、日最低和日最高气温、相对湿度和地表净辐射七种气象变量对径流长期趋势的贡献。研究表明气温对于径流有着最为重要的影响,因此输入中加入气温能够有效提高径流模拟精度。
基金项目
国家自然科学基金项目“基于多变量多尺度数据反演变饱和带水力参数的理论与方法”(项目编号:42272280);中央高校基本科研业务费“深时数字地球”前沿科学中心“深时数字地球”中央高校科技领军人才团队项目(项目编号:2652023001);中国科学院战略性先导科技专项(项目编号:XDA20100103)。
致 谢
感谢硕士研究生:樊雨萱、邢琳婧在图1和图2制作中的贡献。
NOTES
作者简介:李娜(1981-),女,陕西汉中人,博士,副教授,主要从事非饱和带土壤水文学研究。Email: lina@cugb.edu.cn
*通讯作者Email: liangsh@cugb.edu.cn