CN106257948B - 一种流域降水监测无线传感器网络节点优化布局方法 - Google Patents

Abstract

本发明涉及一种流域降水监测无线传感器网络节点优化布局方法，包括：相关性计算；回归分析；去冗余站点；增设站点；增设站点布局优化计算；相似度计算；得到站点布局最佳方案。本发明对传统抽站法进行改进，基于改进的抽站法对现有降水监测站点点进行站网分析，去除冗余站点，在此基础上直接考虑路网对站点布设可行性的限制，确定增设站点的备选区，在备选区上增设站点，综合利用目标变量降水量和经纬度、高程、坡度以及坡向等环境变量之间的相关关系，以回归克里格方差为目标函数建立站点布局优化模型，基于并行设计的模拟退火算法，通过高性能求解站点布局优化模型，实现降水监测站点的优化布局。

Description

一种流域降水监测无线传感器网络节点优化布局方法

技术领域

本发明涉及一种流域水循环监测无线传感器网络节点优化布局方法，是一种合理布置降水监测站点的方法。

背景技术

降水作为流域水循环的重要驱动因子，其准确监测对防洪减灾、水资源管理、水文预报调度具有重要的实践价值。目前，降水主要依靠气象站、降水监测站点、水文站等传统地面监测和卫星遥感获取。传统站网监测中站点建设和维护成本高，难以建立高密度降水监测站网，尤其高山峡谷区，难以布设站点，实现降水的时空精细观测十分困难，且野外自然条件恶劣、气温过低以及人为破坏等多种因素均会引起仪器故障，导致无法弥补的数据损失。卫星遥感产品虽然可以获取大范围降水数据，但数据的真实性需利用地面观测数据进行检验，由于卫星和地面观测在时空尺度上存在差异，为此需要考虑降水的时空变异特征，加密布置地面监测站点，获取遥感像元尺度的地面降水观测网，利用地统计方法获取像元降水真值，以有限且具有代表性的地面观测降水数据检验遥感降水观测在像元尺度上的真实性。目前这些观测手段难以满足精细尺度下流域水循环过程观测、模拟以及遥感降水产品真实性检验对地面降水的观测需求。集传感器技术，自动控制技术，数据网络传输、存储、处理与分析技术于一体的无线传感器网络技术具有低成本、自动化、高精度、时空连续性等特性，通过各种通讯技术可将各个传感器节点动态组网，形成传感器矩阵，克服了传统单点无法连续观测区域尺度降水时空变异特性的劣势，能获取大范围、高时空分辨率的分布式实时观测数据，可以很好捕捉流域尺度下降水的时空动态分布特征及其不确定性，同时可利用通讯技术对各传感器节点的工作状态进行远程监控，基于数据自动诊断实时检查观测数据质量，最大程度避免出现无效数据导致数据缺失问题。目前以自动监测和自动传输的自动降水监测站、自动测报气象站、自动测报水文站组成的降水监测无线传感器节点在获取降水资料中发挥越来越重要的作用。降水监测无线传感器节点即降水监测站网布设密度和位置对捕捉降水的时空变异性具有决定性作用。传统的站网布设应综合考虑两个因素：(1)利用站网捕捉降水时空分布规律的能力；(2)降水监测站点建设的可行性和运行维护的便易性。由于区域发展不平衡及早期降水监测站网的不合理规划，部分区域降水监测站网难以捕捉区域降水分布规律且存在站点冗余情况，为了以最小的站网密度获取流域降水的时空分布信息，需对优化布局降水监测站点。目前降水监测站网优化布局中常采用传统抽站法对站网进行密度分析进行站点优化，在应用中存在以下不足：(1)传统抽站法以获取的降水数据均值作为该区域降水近似真值，以此为标准计算组合站点降水均值相对于近似真值的相对误差，在允许误差范围内对现有稠密降水监测站网进行精简，剔除冗余站点，该方法仅考虑降水量的大小而忽略了降水监测站点捕捉降水空间分布趋势的能力；(2)在精简现有站点基础上，如何根据高山丘陵区的降水特点，选择有效的方法并考虑合理的约束优化布局增设降水监测站点，成为数字流域信息获取研究的重点和难点；(3)模拟退火算法、遗传算法等启发式算法虽然能提高站点布局优化效率，但耗时较长，计算效率亟待进一步提高。

发明内容

为了克服现有技术的问题，本发明提出了一种流域降水监测无线传感器网络节点优化布局方法。所述的方法对传统抽站法进行了改进，基于改进的抽站法对现有降水监测站网进行密度分析，在此基础上考虑路网对站点布设可行性的限制，确定增设站点的备选区，在备选区上增设站点，综合利用目标变量降水量和经纬度、高程、坡度以及坡向等环境变量之间的相关关系，以回归克里格方差为目标函数建立站点布局优化模型，采用基于并行设计的模拟退火算法对站点布局优化模型高性能求解，实现降水监测的优化布局。该方法对传统抽站法进行改进，基于改进的抽站法去除冗余站点，在此基础上直接考虑路网对站点布设可行性的限制，基于模拟退火算法，通过高性能求解站点布局优化模型，对站点进行更加高效合理的优化布局。

本发明的目的是这样实现的：一种流域降水监测无线传感器网络节点优化布局方法，所述方法的步骤如下：

相关性计算的步骤：用于利用现有n个站点数据，以年平均降水量为目标变量，计算年平均降水量和环境变量之间的相关性，所述的环境变量包括：经度、纬度、高程、坡度、坡向；回归分析的步骤：用于选择和年平均降水量显著相关的环境变量和年平均降水量进行广义最小二乘回归分析，并计算回归残差的变异函数；

去冗余站点的步骤：用于以组合站点和现有站点获取降水量分布概率密度函数的相似度作为评价标准对现有降水监测站点进行站网分析，去除冗余站点，精简站点至m个；

基于路网约束增设站点的步骤：用于考虑路网对站点布设可行性的限制，将距道路小于等于3公里的区域作为增设站点的备选区，基于去冗余后剩余的m个站点，在备选区依次增设1个站点；

增设站点布局优化模型建立的步骤：用于分别以m+1、…、m+N_n个站点的环境变量分别计算回归克里格方差，以回归克里格方差为目标函数建立站点布局优化模型；

优化模型求解的步骤：基于并行计算设计的模拟退火算法求解站点布局优化模型，对增设的1，2，…，N_n个站点进行优化布局；

相似度计算的步骤：用于计算增设站点优化布局后优化布局站点获取的年平均降水量分布概率密度函数和整个研究区域年平均降水量分布概率密度函数的相似度；

得到站点布局最佳方案的步骤：用于根据所述增设1，2，…，N_n个站点和研究区降水分布的概率密度函数的相似度，选择相似度达到最大值时的站点数N_k，其中：N_k∈[N₁，N₂，…，N_n]，作为增设的最佳站点数，此时N_k个站点位置为研究区最佳降雨监测站点布设方案。

进一步的，所述“相关性计算的步骤”中计算相关性包括如下子步骤：

根据所述现有站点获取的降水数据计算各个站点处的年平均降水；

根据所述现有站点所在研究区的DEM数据，获取研究区坡度、坡向数据；

根据所述现有站点的经纬度信息，提取现有站点处的高程、坡度以及坡向；

根据所述现有站点处的年平均降水量数据和环境变量数据进行相关性分析，计算目标变量和环境变量之间的相关系数以及显著性水平。

进一步的，所述的“回归分析的步骤”中所述的回归分析包括如下子步骤：

根据所述目标变量和环境变量之间的相关系数以及显著性水平，选择和年平均降水量显著相关的环境变量，即显著性水平p<0.05的环境变量，和年平均降水量进行广义最小二乘回归分析，计算回归系数；

根据所述广义最小二乘回归系数，计算现有站点处降水量回归残差，计算回归残差的变异函数。

进一步的，所述“相似度计算的步骤”中计算增设站点优化布局后优化布局站点获取的年平均降水量分布概率密度函数和整个研究区域年平均降水量分布概率密度函数的相似度包括如下子步骤：

根据所述现有站点降水数据采用普通克里格插值获取整个研究区的降水量数据，计算其分布概率密度函数，即降水总体的分布概率密度函数；

根据所述增设1，2，…，N_n个站点优化布局结果，利用优化站点所在经纬度信息提取m+1、…、m+N_n个站点降水量，计算其分布概率密度函数；

根据所述优化站点和研究区降水分布的概率密度函数，计算两者的相似度。

进一步的，所述“优化模型求解的步骤”中求解包括如下子步骤：

步骤1：初始化种群，在所有可布设站点位置进行编号：1、2、…、VarMax，从中随机选择N组站点，每组站点数为1，每组站点为一个种群，提取每个种群代表的站点处和目标变量显著相关的环境变量；

步骤2：计算各个种群的目标函数值，将目标函数值最小的种群定义为全局最优解，根据残差的变异函数和回归克里格方差计算公式，计算在去冗余剩余m个站点基础上增设每个种群中的1个站点情况下共m+1个样本点的估计回归克里格方差，选择回归克里格方差最小的种群为初始全局最优解；

步骤3：随机移动站点位置，进行种群更新，计算新种群的目标函数值，和进化前的种群目标函数进行对比，若新种群目标函数值小于进化前种群，则保留新种群，此时新种群为该种群的个体最优解，若新种群目标函数值大于等于进化前种群，为了保证算法有很好的全局搜索能力，以一定概率接受新种群；随机移动站点位置，进行种群更新，该过程在MATLAB中采用并行设计，包括：

根据计算机核数P对种群更新任务进行分解，分解思路：利用MATLAB matlabpoollocal P函数开启并行计算环境，生成P个workers，采用spmd并行结构，利用workers索引labindex将步骤1中初始化的N个种群分解，给每个worker分配部分种群，利用P个workers同时对种群进行更新；每次迭代后将P个workers种群进行合并，确定全局最优解，然后继续进行种群分解，利用P个workers对种群进行更新，循环迭代，实现种群更新的并行计算，并行计算效率采用加速比和加速效率进行衡量，公式如下：

加速比：

加速效率：

式中，T_串行、T_并行分别代表串行优化和并行优化所需时间，P为workers的个数。

根据公式0.1*VarMax*randn(n)进行种群位置更新，randn(n)生成n个[-1,1]范围之间的随机数，若超出[1,VarMax]范围，则重新更新种群位置，直至新种群位于[1,VarMax]范围内；

若新种群目标函数值大于等于进化前种群，为了保证算法有很好的全局搜索能力，以一定概率接受新种群，概率P计算公式如下：

DELTA＝(newpop(i).Cost-pop(i).Cost)/pop(i).Cost

P＝exp(-DELTA/T_Iter)

式中，newpop(i).Cost为进化后种群i的目标函数值，pop(i).Cost为进化前种群i的目标函数值，T_Iter为温度，并按下式衰减：

T_Iter＝0.99*T₀

式中，T₀为初始温度，T_Iter随迭代次数的叠加逐渐减小；

概率P随着迭代次数增加会逐渐减小，局部搜索能力不断增强，使种群收敛于近似最优解；步骤4：循环迭代，直至达到最大迭代次数，最终获取的全局最优解为最佳布设方案。

本发明产生的有益效果是：本发明对传统抽站法进行改进，基于改进的抽站法对现有降水监测站点进行站网分析，去除冗余站点，直接考虑路网对站点布设可行性的限制，确定增设站点的备选区，在备选区增设站点，以回归克里格方差为目标函数建立站点布局优化模型，基于并行设计的模拟退火算法，通过高性能求解站点布局优化模型，对增设站点进行优化布局，计算增设站点优化布局后站点处和整个研究区域年平均降水量分布概率密度函数的相似度，高效率的选择相似度达到最大值时所需增设的站点数，找到降水监测站点最佳布设方案。与传统方式相比，本发明所述方法优化效率更高，优化布局结果更合理。

附图说明

下面结合附图和实施例对本发明作进一步说明。

图1是本发明的实施例一所述降水监测站点优化布局方法流程图；

图2是本发明的实施例一所述应用实例中的雅砻江流域(部分)降水监测站点分布图；

图3是本发明的实施例一所述应用实例中的雅砻江流域(部分)组合站点和现有站点降水概率密度函数相似度随组合站点数变化图；

图4是本发明的实施例一所述应用实例中的雅砻江流域现有站点去冗余剩余78个站点分布图；

图5是本发明的实施例一所述应用实例中的去冗余剩余78个站点和现有155个站点降水概率密度对比图；

图6是本发明的实施例一所述应用实例中的雅砻江流域基于路网限制的可布设降水监测站点备选区示意图；

图7是本发明的实施例一所述应用实例中的降水监测站点年平均降水量和降水量总体分布概率密度函数的相似度随新增降水监测站点数目N_n的变化图；

图8是本发明的实施例一所述应用实例中的原始站点、优化布局站点和总体的降水概率密度分布对比图；

图9是本发明的实施例一所述应用实例中的雅砻江流域降水监测站点优化布局图。

具体实施方式

实施例一：

本实施例是一种流域降水监测无线传感器网络节点优化布局方法。所述方法对传统的抽站法进行改进，基于改进的抽站法对现有站网进行密度分析，去除冗余站点，在此基础上直接考虑路网对站点布设可行性的限制，确定可增设站点备选区，在备选区域依次增设站点，综合利用目标变量降水量和经纬度、高程、坡度以及坡向等环境变量之间的相关关系，以回归克里格方差为目标函数建立站点布局优化模型，基于并行设计的模拟退火算法，通过高性能求解站点布局优化模型，实现降水监测站点的优化布局，以优化布局站点和研究区降水概率密度函数的相似度达到最大值时增设的站点数为增设的最佳站点数，此时增设站点布局即为站点的最佳布设方案。该方法对传统抽站法进行改进，基于改进的抽站法去除冗余站点，同时考虑站网对布设站点可行性的限制，基于并行设计的模拟退火算法，通过高性能求解站点布局优化模型，实现对站点更加高效合理的优化布局。

本实施例所述方法的原理为：

利用现有n个站点数据计算目标变量年平均降水量和经度、纬度、高程、坡度、坡向等环境变量之间的相关性。

选择和年平均降水量显著相关(p<0.05)的环境变量和年平均降水量进行广义最小二乘(GLS)回归分析，利用GS+软件计算回归残差的变异函数；

对传统抽站法进行改进，基于改进的抽站法对现有降水监测站点进行站网分析，去除冗余站点；

考虑路网对站点布设可行性的限制，将距道路小于等于3公里的区域作为增设站点的备选区，在去冗余后剩余m个站点基础上在备选区依次增设1个站点，以回归克里格方差为目标函数建立站点布局优化模型，基于并行设计的模拟退火算法，通过高性能求解站点布局优化模型，对增设的站点进行优化布局。

统计增设1，2，…，N_n站点后站点获取的年平均降水量分布的概率密度函数，计算优化布局后站点处和整个研究区域年平均降水量分布概率密度函数的相似度。

选择相似度达到最大值时增设的站点数为增设的最佳站点数，此时站点的布局为降水监测站点最佳布设方案。

对传统抽站法进行改进，基于改进的抽站法对现有降水监测站点进行站网分析，去除冗余站点包括：

以组合站点和现有站点获取降水量分布的概率密度函数的相似度代替传统抽站法中以降水量均值的相对误差作为评价标准对传统的抽站法进行改进，该标准综合考虑了降水量的大小以及空间分布规律。

组合站点是指随机抽取数目为1、2、…、N个站点进行组合，N为现有站点个数，为了保证随机抽站法的稳定性，每个数目的组合站点，随机抽取100个组合。

组合站点和现有站点获取降水量分布的概率密度函数的相似度是指组合站点和现有站点获取降水量分布的概率密度函数重叠部分面积，选择随机抽取的100个组合中和总体相似度最大的组合为该站点数的最终组合。

评价标准是指利用组合站点和现有站点获取的多年平均降水概率密度函数图重叠部分面积的均值作为该站点和现有站点之间的相似度，控制误差在5％以内的最小组合站点数为可以代表现有站点最佳站点数，去冗余剩余m个站点。

考虑路网对站点布设可行性的限制，将距道路小于等于3公里的区域作为增设站点的备选区，在去冗余后剩余m个站点基础上在备选区依次增设1个站点，以回归克里格方差为目标函数建立站点布局优化模型，基于并行设计的模拟退火算法，通过高性能求解站点布局优化模型，对增设的站点进行优化布局包括：

考虑路网对站点布设可行性的限制，将距道路小于等于3公里的区域作为增设站点的备选区包括：

基于研究区的路网，在Arcgis中设置缓冲区，缓冲距离设置设为3公里，得到备选区的矢量范围，利用该矢量范围对研究区DEM数据进行裁剪，获得备选区DEM数据；

回归克里格方差计算公式为：

δ²(s₀)＝C(0)-c′₀c^-1c₀+x′_α(X′c^-1X)^-1x_α

式中，x_ɑ＝x₀－X'C^-1c₀；C(0)为残差变异函数图的基台值；其中C(0)－c′₀c^-1c₀为残差的估计误差的方差；x'_ɑ(X'C^-1X)^-1x_ɑ为趋势项的估计误差的方差。

在去冗余后剩余m个站点基础上在备选区依次增设1个站点包括：

步骤1：在去冗余后剩余m个站点基础上在备选区内随机选择1个站点，以回归克里格方差为目标函数建立站点布局优化模型，基于并行设计的模拟退火算法，通过高性能求解站点布局优化模型，进行站点优化布局。

步骤2：固定优化布局后的该站点位置，和去冗余后剩余m个站点合并作为新的样本点，在m+1个样本点的基础上在备选区继续增设1个站点，以回归克里格方差为目标函数，利用基于并行设计的模拟退火算法对增设站点进行优化布局。

步骤3：重复步骤2，每次增设站点最优位置确定后，固定该站点，加入样本点，在此基础上继续增设1个站点，依次类推。

利用基于并行设计的模拟退火算法对增设的站点进行优化布局包括：

步骤1：初始化种群。

步骤2：计算各个种群的目标函数值，将目标函数值最小的种群定义为全局最优解。

步骤3：随机移动站点位置，进行种群更新，计算新种群的目标函数值，和进化前的种群目标函数进行对比，若新种群目标函数值小于进化前种群，则保留新种群，此时新种群为该种群的个体最优解，若新种群目标函数值大于等于进化前种群，为了保证算法有很好的全局搜索能力，以一定概率接受新种群。

步骤4：循环迭代，直至达到最大迭代次数，最终获取的全局最优解为最佳布设方案。

步骤1中，在所有可布设站点位置进行编号：1、2、…、VarMax，从中随机选择N组站点，每组站点数为1，每组站点为一个种群，提取每个种群代表的站点处和目标变量显著相关的环境变量。

步骤2中，根据残差的变异函数和回归克里格方差计算公式，计算在去冗余剩余m个站点基础上增设每个种群中的1个站点情况下共m+1个样本点的估计回归克里格方差，选择回归克里格方差最小的种群为初始全局最优解。

步骤3中，随机移动站点位置，进行种群更新，该过程在MATLAB中采用并行设计，包括：

根据计算机核数P对种群更新任务进行分解，分解思路：利用MATLAB matlabpoollocal P函数开启并行计算环境，生成P个workers，采用spmd并行结构，利用workers索引labindex将步骤1中初始化的N个种群分解，给每个worker分配部分种群，利用P个workers同时对种群进行更新；每次迭代后将P个workers种群进行合并，确定全局最优解，然后继续进行种群分解，利用P个workers对种群进行更新，循环迭代，实现种群更新的并行计算，并行计算效率采用加速比和加速效率进行衡量，公式如下：

加速比：

加速效率：

式中，T_串行、T_并行分别代表串行优化和并行优化所需时间，P为workers的个数。

步骤3中，根据公式0.1*VarMax*randn(n)进行种群位置更新，randn(n)生成n个[-1,1]范围之间的随机数，若超出[1,VarMax]范围，则重新更新种群位置，直至新种群位于[1,VarMax]范围内。

步骤3中，若新种群目标函数值大于等于进化前种群，为了保证算法有很好的全局搜索能力，以一定概率接受新种群，概率P计算公式如下：

DELTA＝(newpop(i).Cost-pop(i).Cost)/pop(i).Cost

P＝exp(-DELTA/T_Iter)

式中，newpop(i).Cost为进化后种群i的目标函数值，pop(i).Cost为进化前种群i的目标函数值，T_Iter为温度，并按下式衰减：

T_Iter＝0.99*T₀

式中，T₀为初始温度，T_Iter随迭代次数的叠加逐渐减小。

概率P随着迭代次数增加会逐渐减小，局部搜索能力不断增强，使种群收敛于近似最优解。

统计增设1，2，…，N_n站点后站点获取的年平均降水量分布的概率密度函数，计算优化布局后站点处和整个研究区域年平均降水量分布概率密度函数的相似度包括：

首先基于并行设计的模拟退火算法，通过高性能求解站点布局优化模型，确定增设1个站点的最优布局并固定，然后在m+1个站点基础上增设1个站点，确定该站点位置，依次确定N_n个站点位置。

基于N个现有站点降水数据采用普通克里格插值获取整个研究区的降水分布。

根据确定的N_n个站点位置以及去冗余剩余m个站点位置获取m+N_n个站点位置处获取的年平均降水量，统计增设1，2，…，N_n站点获取的年平均降水概率密度函数，计算与研究区年平均降水总体分布的概率密度函数的相似度。

选择相似度达到最大值时增设的站点数为增设的最佳站点数，此时站点的布局为降水监测站点最佳布设方案包括：

根据获取的增设站点优化布局后降水监测站点获取的年平均降水概率密度函数与研究区年年平均平均降水总体分布的概率密度函数的相似度随增设站点的变化图，确定相似度达到最大值时增设的站点数为最佳布设站点个数，此时站点的布局为降水监测站点最佳布设方案。

所述方法的具体步骤如下：

相关性计算的步骤：用于利用现有n个站点数据，以年平均降水量为目标变量，计算年平均降水量和环境变量之间的相关性，所述的环境变量包括：经度、纬度、高程、坡度、坡向。

回归分析的步骤：用于选择和年平均降水量显著相关的环境变量和年平均降水量进行广义最小二乘回归分析，并计算回归残差的变异函数。

去冗余站点的步骤：用于以组合站点和现有站点获取降水量分布的概率密度函数的相似度作为评价标准对现有降水监测站点进行站网分析，去除冗余站点，精简站点至m个。

基于路网约束增设站点的步骤：用于考虑路网对站点布设可行性的限制，将距道路小于等于3公里的区域作为增设站点的备选区，去冗余后剩余m个站点基础上，在备选区依次增设1个站点。

增设站点布局优化模型建立的步骤：用于分别以m+1、…、m+N_n个站点的环境变量分别计算回归克里格方差，以回归克里格方差为目标函数建立站点布局优化模型。

优化模型求解的步骤：基于并行设计的模拟退火算法，通过求解站点布局优化模型，对增设的1，2，…，N_n个站点进行优化布局；

相似度计算的步骤：用于计算增设站点优化布局后，优化布局站点获取的年平均降水量分布概率密度函数和整个研究区域年平均降水量分布概率密度函数的相似度。

得到站点布局最佳方案的步骤：用于根据所述增设1，2，…，N_n个站点和研究区降水分布的概率密度函数的相似度，选择相似度达到最大值时的站点数N_k，其中：N_k∈[N₁，N₂，…，N_n]，作为增设的最佳站点数，此时N_k个站点位置为研究区最佳降水监测站点布设方案。

使用实例：

以雅砻江流域(部分，下同)降水监测无线传感器网络节点即降水监测站点(图2)优化布局为例来说明。

根据雅砻江流域现有155个站点数据计算各个站点处的年平均降水。根据砻江流域DEM数据，在Arcgis中利用3D Analyst Tools工具箱中的slope和aspect工具获取研究区坡度、坡向数据。雅砻江流域DEM数据空间分辨率为30米。根据现有站点的经纬度信息，提取现有站点处的高程、坡度以及坡向。

根据现有155个站点处的目标变量年平均降水量数据和经纬度、高程、坡度、坡向等环境变量数据进行相关性分析，利用SPSS软件计算目标变量和环境变量之间的相关系数以及显著性水平，分析结果见表1。

表1目标变量与环境变量之间的相关性

rainfall	0.587**	-0.275	-0.335*	0.052	-0.034
							longitude	-0.217	-0.628**	0.001	0.305*
		latitude	0.425**	0.309*	-0.067
									elevation	0.214	-0.021
				slope	-0.151
											aspect

其中，*、**分别代表显著性水平为0.05、0.01。

选择和年平均降水量显著相关(p<0.05)的环境变量和年平均降水量进行广义最小二乘(GLS)回归分析，利用GS+软件计算回归残差的变异函数包括：

根据目标变量与环境变量相关性分析结果，选择和年平均降水量显著相关(p<0.05)的经度和高程和年平均降水量进行广义最小二乘(GLS)回归分析，计算回归系数，计算公式如下：

β＝(X'C^-1X)^-1X'C^-1z(s_n)

式中，z(s_n)为155×1的样本点处年平均降水量矩阵，C为155×155的回归残差的方差-协变量矩阵，X为155×3的样本点环境变量矩阵。

由回归系数，可得到其他位置降水的最优线性无偏估计值(BLUE)：

式中，x₀为待估计点s₀处经度、高程矩阵。c₀为s₀和155个样本点处年平均降水量的协方差矩阵。根据残差的变异函数可获取c₀和C。

计算回归残差：

ε(s)＝Z(s)-x′β

式中，Z(s)为目标变量矩阵，x为环境变量矩阵，β为回归系数。

利用GS+软件计算回归残差的变异函数，变异函数拟合选择指数模型，拟合结果如下：

式中，h为样点之间的距离。

对传统抽站法进行改进，基于改进的抽站法对现有降水监测站点点进行站网分析，去除冗余站点包括：

以组合站点和现有站点获取降水量分布的概率密度函数的相似度代替传统抽站法中以降水量均值的相对误差作为评价标准对传统抽站法进行改进，该标准综合考虑了降水量的大小及其空间分布规律。

组合站点是指随机抽取数目为1、2、…、155个站点进行组合，为了保证随机抽站法的稳定性，每个数目的组合站点，随机抽取100个组合；

组合站点和现有155个站点获取降水量分布的概率密度函数的相似度是指组合站点和现有站点获取降水量分布的概率密度函数重叠部分面积。

组合站点和现有155个站点获取降水量分布的概率密度函数的相似度随组合站点数的变化如图3所示。

评价标准是指利用组合站点和现有站点获取的多年平均降水概率密度函数图重叠部分面积的均值作为该站点和现有站点之间的相似度，控制误差在5％以内的最小组合站点数为最佳站点数，去冗余剩余78个站点。

去冗余剩余78个站点分布如图4所示，去冗余剩余78个站点和研究区(总体)的年平均降水分布概率密度函数对比如图5所示。

考虑路网对站点布设可行性的限制，将距道路小于等于3公里的区域作为增设站点的备选区(图6)，在去冗余后剩余78个站点基础上在备选区依次增设1，2，…，N_n站点，以回归克里格方差为目标函数建立站点布局优化模型，基于模拟退火算法，通过站点布局优化模型高性能求解对增设的站点进行优化布局包括：

回归克里格方差计算公式为：

δ²(s₀)＝C(0)-c′₀c^-1c₀+x′_α(X′c^-1X)^-1x_α

式中，x_ɑ＝x₀－X'C^-1c₀；C(0)为残差变异函数图的基台值；其中C(0)－c'₀c^-1c₀为残差的估计误差的方差；x'_ɑ(X'C^-1X)^-1x_ɑ为趋势项的估计误差的方差。

在去冗余后剩余78个站点基础上在备选区依次增设1，2，…，N_n站点包括：

步骤1：在去冗余后剩余78个站点基础上在备选区增设1个站点，以回归克里格方差为目标函数建立站点布局优化模型，基于模拟退火算法，通过站点布局优化模型高性能求解进行站点优化布局；

步骤2：固定优化布局后的该站点位置，和去冗余后剩余78个站点合并作为新的样本点，在79个样本点的基础上继续增设1个站点，以回归克里格方差为目标函数建立站点布局优化模型，基于模拟退火算法，通过站点布局优化模型高性能求解对增设站点进行优化布局；

步骤3：重复步骤2，每次增设站点最优位置确定后，固定该站点，加入样本点，在此基础上继续增设1个站点，依次类推。

利用采用并行设计的模拟退火算法对增设的站点进行优化布局包括：

步骤1：初始化20个种群；

步骤2：计算各个种群的目标函数值，将目标函数值最小的种群定义为全局最优解；

步骤3：随机移动站点位置，进行种群进化，计算新种群的目标函数值，和进化前的种群目标函数进行对比，若新种群目标函数值小于进化前种群，则保留新种群，此时新种群为该种群的个体最优解，若新种群目标函数值大于等于进化前种群，为了保证算法有很好的全局搜索能力，以一定概率接受新种群；

步骤4：循环迭代，直至达到最大迭代次数，最终获取的全局最优解为最佳布设方案。

步骤1中，在所有可布设站点位置进行编号：1、2、…、VarMax，VarMax为可布设站点最大编号，取值42775，从中随机选择20组个站点，每组站点数为1，每组站点为一个种群，提取每个种群代表的站点处和目标变量显著相关的环境变量；

步骤2中，根据残差的变异函数和回归克里格方差计算公式，计算在去冗余剩余78个站点基础上增设每个种群中的1个站点情况下共79个样本点的估计回归克里格方差，选择回归克里格方差最小的种群为初始全局最优解；

步骤3中，对随机移动站点位置，进行种群优化过程，在MATLAB中进行并行设计，包括：

根据计算机核数P＝8对优化任务进行分解，分解思路：利用MATLAB matlabpoollocal 8函数开启并行计算环境，生成8个workers，采用spmd并行结构，利用workers索引labindex将步骤1中初始化的20个种群进行分配，给前6个worker各分配2个种群，第8、9个worker各分配4个种群，利用8个workers同时对种群进行更新；每次迭代后将8个workers种群进行合并，确定全局最优解，然后继续进行种群分解，利用8个workers对种群进行更新，循环迭代，实现种群更新的并行计算，并行计算效率采用加速比和加速效率进行衡量，公式如下：

加速比：

加速效率：

式中，T_串行、T_并行分别代表串行优化和并行优化所需时间，P为workers的个数。

步骤3中，根据公式0.1*VarMax*randn(1)进行种群位置更新，randn(n)生成1个[-1,1]范围之间的随机数，若超出[1,VarMax]范围，则重新更新种群位置，直至新种群位于[1,VarMax]范围内；

步骤3中，若新种群目标函数值大于等于进化前种群，为了保证算法有很好的全局搜索能力，以一定概率接受新种群，概率P计算公式如下：

DELTA＝(newpop(i).Cost-pop(i).Cost)/pop(i).Cost

P＝exp(-DELTA/T_Iter)

式中，newpop(i).Cost为进化后种群i的目标函数值，pop(i).Cost为进化前种群i的目标函数值，T_Iter为温度，并按下式衰减：

T_Iter＝0.99*T₀

式中，T₀为初始温度，取0.01，T_Iter随迭代次数的叠加逐渐减小。

概率P随着迭代次数增加会逐渐减小，局部搜索能力不断增强，使种群收敛于近似最优解。

加速比以增设第37个站点为例，计算机CPU为Intel(R)Core(TM)，八核，加速比为3.91，加速效率为48.81％，加速效果较好，基于模拟退火算法，通过站点布局优化模型高性能求解可以成倍节省优化布局所需时间。

统计增设1，2，…，N_n站点后站点获取的年平均降水量分布的概率密度函数，计算优化布局后站点处和雅砻江流域年平均降水量分布概率密度函数的相似度包括：

首先基于模拟退火算法，通过站点布局优化模型高性能求解确定增设第1个站点的最优布局并固定，然后在79个站点基础上增设第2个站点，确定该站点位置，依次确定N_n个站点位置；

确定的N_n个站点位置以及去冗余剩余78个站点位置获取78+N_n个站点位置处获取的年平均降水量，统计获取年平均降水的概率密度函数，计算与雅砻江流域年平均降水总体分布概率密度函数的相似度；

78+N_n个站点位置处获取年平均降水量和雅砻江流域的年平均降水总体分布概率密度函数的相似度随N_n的变化如图7所示。

选择相似度达到最大值时增设的站点数为增设的最佳站点数，此时站点的布局为降水监测站点最佳布设方案包括：

根据获取的增设站点后降水监测站点获取的概率密度函数与研究区年平均降水总体分布的概率密度函数的相似度作相似度随增设站点的变化图，确定相似度达到最大值时增设的28个站点为最佳布设站点个数，原始站点、优化布局站点和总体的降水概率密度函数与对比如图8所示，相比原始站点，优化布局站点可以更好的捕捉整个研究区降水量分布规律，此时站点的布局(图9)为降水监测站点最佳布设方案。

实施例二：

本实施例是实施例一的改进，是实施例一关于“相关性计算的步骤”中计算相关性的细化。本实施例所述“相关性计算的步骤”中计算相关性包括如下子步骤：

根据所述现有站点获取的降水数据计算各个站点处的年平均降水。

根据所述现有站点所在研究区的DEM数据，获取研究区坡度、坡向数据；

根据所述现有站点的经纬度信息，提取现有站点处的高程、坡度以及坡向；

根据所述现有站点处的年平均降水量数据和环境变量数据进行相关性分析，计算目标变量和环境变量之间的相关系数以及显著性水平。

实施时，利用现有站点数据计算目标变量年平均降水量和经纬度、高程、坡度、坡向等环境变量之间的相关性包括：

根据所述现有站点获取的降水数据计算各个站点处的年平均降水。

根据现有站点所在研究区的DEM数据，在Arcgis中利用3D Analyst Tools工具箱中的slope和aspect工具获取研究区坡度、坡向数据。

根据现有站点的经纬度信息，提取现有站点处的高程、坡度以及坡向。

根据现有站点处的目标变量年平均降水量数据和经纬度、高程、坡度、坡向等环境变量数据进行相关性分析，利用SPSS软件计算目标变量和环境变量之间的相关系数以及显著性水平。

实施例三：

本实施例是上述实施例的改进，是上述实施例关于“回归分析的步骤”中所述的回归分析的细化。本实施例所述的“回归分析的步骤”中所述的回归分析包括如下子步骤：

根据目标变量和环境变量之间的相关系数以及显著性水平，选择和年平均降水量显著相关的环境变量，即相关系数p<0.05的环境变量，和年平均降水量进行广义最小二乘回归分析，计算回归系数。

根据广义最小二乘回归系数，计算现有站点处降水量回归残差，利用GS+软件计算回归残差的变异函数。

实施时，选择和年平均降水量显著相关(p<0.05)的环境变量和年平均降水量进行广义最小二乘(GLS)回归分析，利用GS+软件计算回归残差的变异函数包括：

选择和年平均降水量显著相关(p<0.05)的环境变量和年平均降水量进行广义最小二乘(GLS)回归分析，计算回归系数，计算公式如下：

β＝(X'C^-1X)^-1X'C^-1z(s_n)

式中，z(s_n)为n×1的样本点处目标变量矩阵，C为n×n的残差的方差-协变量矩阵，X为n×m+1的样本点环境变量矩阵。

根据回归系数，可得到其他位置降水的最优线性无偏估计值(BLUE):

式中，x₀为待估计点s₀处环境变量矩阵；c₀为s₀和n个样本点处目标变量的协方差矩阵；根据残差的变异函数可获取c₀和C。

计算回归残差：

ε(s)＝Z(s)-x′β式中，Z(s)为目标变量矩阵，x为环境变量矩阵，β为回归系数。

实施例四：

本实施例是上述实施例的改进，是上述实施例关于“相似度计算的步骤”中年平均降水量分布概率密度函数的相似度计算的细化。本实施例所述“相似度计算的步骤”中年平均降水量分布概率密度函数的相似度计算包括如下子步骤：

根据所述现有站点降水数据采用普通克里格插值获取整个研究区的降水量数据，计算其分布概率密度函数。

根据所述增设1，2，…，N_n个站点优化布局结果，根据优化站点所在经纬度信息提取m+1、…、m+N_n个站点降水量，计算其分布概率密度函数。

根据所述优化站点和研究区降水分布的概率密度函数，计算两者的相似度。

实施时，统计增设1，2，…，N_n站点后站点获取的年平均降水量分布的概率密度函数，计算优化布局后站点处和整个研究区域年平均降水量分布概率密度函数的相似度包括：首先基于模拟退火算法，通过高性能求解站点布局优化模型，确定增设1个站点的最优布局并固定，然后在m+1个站点基础上增设1个站点，确定该站点位置，依次确定N_n个站点位置；

根据N个现有站点降水数据采用普通克里格插值获取整个研究区的降水分布；

确定的N_n个站点位置以及去冗余剩余m个站点位置获取m+N_n个站点位置处获取的年平均降水量，统计增设1，2，…，N_n站点获取的年平均降水概率密度函数，计算与研究区年平均降水总体分布的概率密度函数的相似度。

选择相似度达到最大值时增设的站点数为增设的最佳站点数，此时站点的布局为降水监测站点最佳布设方案包括：

根据获取的增设站点后降水监测站点获取的概率密度函数与研究区年平均降水总体分布的概率密度函数的相似度作相似度随增设站点的变化图，确定相似度达到最大值时增设的站点数为最佳布设站点个数，此时站点的布局为降水监测站点最佳布设方案。

最后应说明的是，以上仅用以说明本发明的技术方案而非限制，尽管参照较佳布置方案对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案(比如所使用的软硬件平台、步骤的先后顺序等)进行修改或者等同替换，而不脱离本发明技术方案的精神和范围。

Claims (4)

1.一种流域降水监测无线传感器网络节点优化布局方法，其特征在于，所述方法的步骤如下：

相关性计算的步骤：用于利用现有n个站点数据，以年平均降水量为目标变量，计算年平均降水量和环境变量之间的相关性，所述的环境变量包括：经度、纬度、高程、坡度、坡向；回归分析的步骤：用于选择和年平均降水量显著相关的环境变量和年平均降水量进行广义最小二乘回归分析，并计算回归残差的变异函数；

去冗余站点的步骤：用于以组合站点和现有站点获取降水量分布的概率密度函数的相似度作为评价标准对现有降水监测站点进行站网分析，去除冗余站点，精简站点至m个；

基于路网约束增设站点的步骤：用于考虑路网对站点布设可行性的限制，将距道路小于等于3公里的区域作为增设站点的备选区，基于去冗余后剩余m个站点，在备选区依次增设1个站点；

增设站点布局优化模型建立的步骤：用于分别以m+1、…、m+N_n个站点的环境变量分别计算回归克里格方差，以回归克里格方差为目标函数建立站点布局优化模型；

优化模型求解的步骤：基于并行设计的模拟退火算法，通过求解站点布局优化模型，对增设的1，2，…，N_n个站点进行优化布局；

本步骤包括如下子步骤：

子步骤1：初始化种群，在所有可布设站点位置进行编号：1、2、…、VarMax，从中随机选择N组站点，每组站点数为1，每组站点为一个种群，提取每个种群代表的站点处和目标变量显著相关的环境变量；

子步骤2：计算各个种群的目标函数值，将目标函数值最小的种群定义为全局最优解，根据残差的变异函数和回归克里格方差计算公式，计算在去冗余剩余m个站点基础上增设每个种群中的1个站点情况下共m+1个样本点的估计回归克里格方差，选择回归克里格方差最小的种群为初始全局最优解；

子步骤3：随机移动站点位置，进行种群更新，计算新种群的目标函数值，和进化前的种群目标函数进行对比，若新种群目标函数值小于进化前种群，则保留新种群，此时新种群为该种群的个体最优解，若新种群目标函数值大于等于进化前种群，为了保证模拟退火算法有很好的全局搜索能力，以概率P₀接受新种群；

随机移动站点位置，进行种群更新，种群更新在MATLAB中采用并行设计，包括：

根据计算机核数P对种群更新任务进行分解，分解思路：利用MATLAB matlabpoollocal P函数开启并行计算环境，生成P个workers，采用spmd并行结构，利用workers索引labindex将步骤1中初始化的N个种群分解，给每个worker分配部分种群，利用P个workers同时对种群进行更新；每次迭代后将P个workers种群进行合并，确定全局最优解，然后继续进行种群分解，利用P个workers对种群进行更新，循环迭代，实现种群更新的并行计算，并行计算效率采用加速比和加速效率进行衡量，公式如下：

加速比：

加速效率：

式中，T_串行、T_并行分别代表串行优化和并行优化所需时间，P为workers的个数；

根据公式0.1*VarMax*randn(n)进行种群位置更新，randn(n)生成n个[-1,1]范围之间的随机数，若超出[1,VarMax]范围，则重新更新种群位置，直至新种群位于[1,VarMax]范围内；

若新种群目标函数值大于等于进化前种群，为了保证模拟退火算法有很好的全局搜索能力，以概率P₀接受新种群，概率P₀计算公式如下：

DELTA＝(newpop(i).Cost-pop(i).Cost)/pop(i).Cost

P₀＝exp(-DELTA/T_Iter)

式中，newpop(i).Cost为进化后种群i的目标函数值，pop(i).Cost为进化前种群i的目标函数值，T_Iter为温度，并按下式衰减：

T_Iter＝0.99*T₀

式中，T₀为初始温度，T_Iter随迭代次数的叠加逐渐减小；

概率P₀随着迭代次数增加会逐渐减小，局部搜索能力不断增强，使种群收敛于近似最优解；

子步骤4：循环迭代，直至达到最大迭代次数，最终获取的全局最优解为最佳布设方案；

相似度计算的步骤：用于计算增设站点优化布局后优化布局站点获取的年平均降水量分布概率密度函数和整个研究区域年平均降水量分布概率密度函数的相似度；

得到站点布局最佳方案的步骤：用于根据所述增设1，2，…，N_n个站点和研究区降水分布的概率密度函数的相似度，选择相似度达到最大值时的站点数N_k，其中：N_k∈[N₁，N₂，…，N_n]，作为增设的最佳站点数，此时N_k个站点位置为研究区最佳降水监测站点布设方案。

2.根据权利要求1所述的流域降水监测无线传感器网络节点优化布局方法，其特征在于，所述“相关性计算的步骤”中计算相关性包括如下子步骤：

根据所述现有站点获取的降水数据计算各个站点处的年平均降水；

根据所述现有站点所在研究区的DEM数据，获取研究区坡度、坡向数据；

根据所述现有站点的经纬度信息，提取现有站点处的高程、坡度以及坡向；

根据所述现有站点处的年平均降水量数据和环境变量数据进行相关性分析，计算目标变量和环境变量之间的相关系数以及显著性水平。

3.根据权利要求1所述的流域降水监测无线传感器网络节点优化布局方法，其特征在于，所述的“回归分析的步骤”中所述的回归分析包括如下子步骤：

根据所述目标变量和环境变量之间的相关系数以及显著性水平，选择和年平均降水量显著相关的环境变量，即显著性水平p<0.05的环境变量，和年平均降水量进行广义最小二乘回归分析，计算回归系数；

根据所述广义最小二乘回归系数，计算现有站点处降水量回归残差，计算回归残差的变异函数。

4.根据权利要求1所述的流域降水监测无线传感器网络节点优化布局方法，其特征在于，所述“相似度计算的步骤”中计算增设站点优化布局后优化布局站点获取的年平均降水量分布概率密度函数和整个研究区域年平均降水量分布概率密度函数的相似度包括如下子步骤：

根据所述现有站点降水数据采用普通克里格插值获取整个研究区的降水量数据，计算其分布概率密度函数，即降水总体的分布概率密度函数；

根据所述增设1，2，…，N_n个站点优化布局结果，根据优化站点所在经纬度信息提取m+1、…、m+N_n个站点降水量，计算其分布概率密度函数；

根据所述优化站点和研究区降水分布的概率密度函数，计算两者的相似度。