一种基于BP神经网络模型的年天文辐射量计算方法与流程

本发明涉及地理信息系统的技术领域，尤其涉及一种基于bp神经网络模型的年天文辐射量计算方法。

背景介绍

太阳辐射是地球上一切能量的最初来源，其变化对全球地表能量流动、物质循环、生物生长以及人类活动等都具有决定性的影响。在不考虑大气影响的情形下，由日地天文关系和地形遮蔽因素决定的地表太阳辐射称为天文辐射，它是计算地表实际入射太阳辐射(直接辐射，总辐射等)所必需的基础数据，也是进行太阳能资源评估及其他相关领域研究的重要基础数据。

目前天文辐射的主要计算方法如下：1)使用数字高程模型(digitalelevationmodel，dem)并根据纬度和日期(积日)信息计算天文辐射。但计算的数学模型结构复杂，计算量大，效率低。2)以气象站点的观测值(如气温)作为输入值以相应站点的天文辐射观测值为输出值，训练了对应的bp(backpropagation)神经网络模型，用于估算其它气象站点的天文辐射。但所得到的bp神经网络模型仅可使用于气象站点(其数量极为有限)，而在非气象站点的广大地区无法使用，极大地限制了天文辐射的基础数据作用。

技术实现要素：

本发明的目的在于克服现有技术的不足，提供一种基于bp神经网络模型计算年天文辐射量的方法，能够快速准确地计算出指定地点天文辐射量，解决了以往计算方法中模型复杂、计算量大的问题。

一种基于bp神经网络模型的年天文辐射量的计算方法，包括以下步骤：

步骤1，数据下载以及预处理：下载待计算地区的dem影像，将dem影像转换为数值矩阵，填充数值矩阵中的空缺值；

步骤2，将所述待计算地区的dem影像构建成输入矩阵；

步骤3，对所述输入矩阵归一化处理；

步骤4，将所述归一化后的矩阵输入bp神经网络模型，计算得到输出矩阵；

步骤5，对输出矩阵进行反归一化，转换成所述计算地区的年天文辐射数值矩阵。

进一步地，所述的bp神经网络模型为5层网络，由一个输入层、三个隐含层、一个输出层构成，所述隐含层的神经元数量均为15，所述输出层的神经元数量为1；所述的bp神经网络模型的训练，包括以下步骤：

第一步：选取训练样本；

第二步：构建bp神经网络模型的输入矩阵和期望输出矩阵；

第三步：将上述第二步中的输入矩阵和期望输出矩阵分别归一化；

第四步：构建bp神经网络模型的框架，设置网络基本参数；

第五步：初始化各层权值矩阵和阈值向量；

第六步：将归一化后的输入矩阵，输入bp神经网络模型的输入层，计算出其输出矩阵，求出输出矩阵与归一化后的期望输出矩阵之间的训练误差e；

第七步：判断训练误差e是否达到设定的要求，若e达到要求，则结束训练，并转向第八步；若没有达到要求，则反向修正bp神经网络模型的各层权值矩阵和阈值向量，再转向第六步；

第八步：保存各层权值矩阵和阈值向量，bp神经网络模型构建完成；

进一步地，所下载的dem影像，为tiff格式，所包含像素的数量为p1×p2，像素值表示该像素所对应位置的高程值。根据地理信息系统专业的表达习惯，将dem影像的像素称为栅格。为了方便描述，将每个栅格进行编号(编号从1到p1×p2)，编号按照先从上到下，后从左到右的顺序进行。

进一步地，所述“将所述待计算地区的dem影像构建成输入矩阵”，包括以下步骤：

步骤21：提取待计算地区dem影像每个栅格的高程、坡度、坡向因子以及所述栅格上下左右四个方向邻域栅格的高程、坡度和坡向因子；

步骤22：将每个栅格的经度、纬度、高程、坡度、坡向，以及所述栅格上下左右四个方向邻域栅格的高程、坡度、坡向作为该栅格的17个因子，并将每个栅格对应的17个因子作为一个列向量；

步骤23：将p1×p2个栅格对应的列向量依次按照栅格的编号组合为一个矩阵，该矩阵即为待计算地区dem影像对应的输入矩阵inputm×n

其中，m为输入矩阵input的行数，n为input的列数；i＝1，2，3，…，m，该矩阵中m＝17；j为dem栅格的编号，j＝1，2，3，…，n，该矩阵中n＝p1×p2；x1，j和x2，j为第j个栅格的经度和纬度，取该幅dem影像中心栅格的经度和纬度；x3，j、x4，j、x5，j、x6，j、x7，j分别为第j个栅格和该栅格四个方向邻域栅格的高程；x8，j、x9，j、x10，j、x11，j、x12，j分别为第j个栅格和该栅格四个方向邻域栅格的坡度；x13，j、x14，j、x15，j、x16，j、x17，j分别为第j个栅格和该栅格四个方向邻域栅格的坡向；

进一步地，所述“对所述输入矩阵归一化处理”，使每一行的元素值限定在[-1,1]区间范围内。

进一步地，所述“对输出矩阵进行反归一化，转换成所述计算地区的年天文辐射数值矩阵”，包括以下步骤：

步骤51：对bp神经网络模型的输出矩阵进行反归一化处理；

步骤52：将反归一化后的单行输出矩阵重新组合为p1行p2列的矩阵，组合时每p1个元素为一列，即第1到第p1个元素为第一列，第p1+1到第2*p1个元素为第二列，以此类推，组合成一个p1行p2列的矩阵。

综上所述，本发明具有以下有益效果：

1.与数学模型的计算相比，使用本bp神经网络模型计算年天文辐射具有模型简单易懂、计算量小、计算时间短，计算方法简单等优点。

2.使用本发明所述的bp神经网络模型能够精确计算指定地点，包括非气象站点的年天文辐射，为后续的研究提供更丰富的数据。

3.本发明可为精细化农业提供重要的基础数据即天文辐射，为农业种植区域的规划、农业资源的评估提供依据；为准确、科学地计算评估太阳能资源提供基础数据，为大规模开发利用可再生太阳能资源提供依据；为气象气候学的发展提供基础数据，为区域气候的空间模拟与可视化表达提供了新的研究思路。

4.本发明的bp神经网络训练方法可应用于任意空间分辨率的数字高程模型，训练相应空间分辨率的bp神经网络模型，可推广性高。

附图说明

图1为本发明基于bp神经网络模型的年天文辐射量计算方法的步骤流程图；

图2为本发明将待计算地区的dem影像构建成输入矩阵的步骤流程图；

图3为本发明bp神经网络模型的结构图；

图4为本发明bp神经网络模型训练步骤流程图；

图5为本发明实施例中包含了北京市顺义区的单幅dem影像图；

图6为本发明实施例中包含了北京市顺义区的年天文辐射量分布图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1，是基于bp神经网络模型的年天文辐射量计算方法的步骤流程图，包括如下步骤：

步骤1，数据下载以及预处理：下载待计算地区的dem影像，将dem影像转换为数值矩阵，填充数值矩阵中的空缺值；

在本发明实施例中，数字高程模型(digitalelevationmodel，dem)影像数据来源于中国科学院计算机网络信息中心地理空间数据云平台(http://www.gscloud.cn)。下载待计算地区包含的30m空间分辨率的dem影像，所下载的每一幅dem影像的经度与纬度间隔均为1°，影像格式为tiff格式，每一幅dem影像的像素为3601×3601，根据地理信息系统专业的表达习惯，将一幅dem影像中的每一个像素称为栅格，则一幅30m分辨率的dem影像共有12967201(3601×3601)个栅格。对单幅dem影像的栅格进行编号，按照先从上到下，后从左到右的顺序依次进行编号，例如：第1行第1列的栅格编号为1，第2行第1列的栅格编号为2，第3行第1列的栅格编号为3，…，第3601行3601列的栅格编号为12967201，以此形成一个矩阵。该矩阵的数值为dem影像每一个栅格编号所对应的像素值，即该dem影像像素所在位置的高程值。对数值矩阵存在空缺值的栅格处使用邻近法进行填充，即填入最邻近栅格的高程值。

另外，计算天文辐射时需要dem影像每个栅格的经纬度，由于每幅dem影像的经度和纬度间隔仅为1°，为方便计算，各栅格的经纬度坐标均取dem影像中心栅格的经纬度坐标。

在一个优选的实施例中，可以借助计算机编程语言辅助构建数值矩阵，包括：python、c++、matlab脚本等。

步骤2，将所述待计算地区的dem影像构建成输入矩阵；

如图2是本发明将待计算地区的dem影像构建成输入矩阵的步骤流程图，包括以下步骤：

步骤21：提取待计算地区dem影像每个栅格的高程、坡度、坡向因子以及所述栅格上下左右四个方向邻域栅格的高程、坡度和坡向因子；

步骤22：将每个栅格的经度、纬度、高程、坡度、坡向，以及所述栅格上下左右四个方向邻域栅格的高程、坡度、坡向作为该栅格的17个因子，并将每个栅格对应的17个因子作为一个列向量；

步骤23：将dem影像每个栅格对应的列向量依次按照栅格的编号组合为一个矩阵，该矩阵即为待计算地区dem影像对应的输入矩阵inputm×n

其中，m为输入矩阵input的行数，n为input的列数；i＝1，2，3，…，m，该矩阵中m＝17；j为dem栅格的编号，j＝1，2，3，…，n，该矩阵中n＝p1×p2；x1，j和x2，j为第j个栅格的经度和纬度，取该幅dem影像中心栅格的经度和纬度；x3，j、x4，j、x5，j、x6，j、x7，j分别为第j个栅格和该栅格四个方向邻域栅格的高程；x8，j、x9，j、x10，j、x11，j、x12，j分别为第j个栅格和该栅格四个方向邻域栅格的坡度；x13，j、x14，j、x15，j、x16，j、x17，j分别为第j个栅格和该栅格四个方向邻域栅格的坡向；

步骤3，对所述输入矩阵归一化处理；

在本实施例中，对输入矩阵进行归一化处理，使每一行的元素值限定在[-1，1]区间范围内，可使用如下归一化公式：

式中，i＝1，2，3，…，17；j＝1，2，3，…，n；xi，j是矩阵input第i行第j列元素的值；是矩阵input第i行所有元素中的最大值，是矩阵input第i行所有元素中的最小值，pi，j是归一化后的输入矩阵p第i行第j列元素的值，取值范围是[-1，1]；

在一个优选的实施例中，可使用matlab软件调用归一化函数mapminmax完成，函数调用方式如下：

[p，pn]＝mapminmax(input，pmin，pmax)

式中input为输入矩阵，p为输出矩阵，pn为调用mapminmax函数时得到的参数，主要包含了归一化前矩阵每一行全部元素的最小值和最大值，以及输出矩阵p的最大值pmax和最小值pmin(pmin默认为-1，pmax默认为1)。

步骤4，将所述归一化后的矩阵输入bp神经网络模型，计算得到输出矩阵；

如图3是本发明bp神经网络模型结构图，所述的bp神经网络模型为5层网络，由一个输入层、三个隐含层、一个输出层构成，所述隐含层的神经元数量均为15，所述输出层的神经元数量为1。

在本实施例中，p为归一化后的输入矩阵；wl和bl为bp神经网络模型第l层隐含层的权值矩阵和阈值向量；r为输入层神经元数量即p的行数，该模型中r＝17；s1，s2，s3为隐含层神经元数量，该模型中s1＝s2＝s3＝15；s4为输出层神经元数量，该模型中s4＝1；q为p的列数即dem中栅格的个数，该模型中q＝12967201；h1＝1，2，3，…，s1；h2＝1，2，3，…，s2；h3＝1，2，3，…，s3；h4＝1；

f1，f2，f3为隐含层的激活函数，是非线性函数，公式为：

式中l＝1，2，3；al为bp神经网络模型第l层隐含层的输出矩阵，al-1为bp神经网络模型第l层隐含层的输入矩阵，二者的取值范围是[-1，1]，当l＝1时，a0＝p，p为bp神经网络模型输入层的输入矩阵；

f4为输出层的激活函数，是线性函数，公式为：

y＝f4(wl×x+bl)＝wl×x+bl

式中，y为bp神经网络模型输出层的输出矩阵，x为bp神经网络模型输出层的输入矩阵；该式中l＝4，则o＝w4*a3+b4＝a4；

o和a4为bp神经网络模型输出层的输出矩阵(1行12967021列)，a3为bp神经网络模型输出层的输入矩阵(1行12967021列)。

如图4是本发明bp神经网络模型训练步骤流程图，包括以下步骤：

第一步：选取训练样本；

在本实施例中，选取的训练样本，包括250幅dem影像，编号为1-250，并在每个dem影像随机抽取40个栅格，此后称抽取的栅格为“采样栅格”；

第二步：构建bp神经网络模型的输入矩阵和期望输出矩阵；

在本实施例中，构建的bp神经网络模型输入矩阵的具体步骤为：

构建编号为1的dem影像对应的矩阵l1，

式中，l1为17行40列的矩阵；l1的下标表示dem影像的编号；d＝1，2，3，...，40；矩阵l1中第d列的元素x1d，x2d，...，x17,d依次为第d个采样栅格的经度值、纬度值、高程值、四个方向邻域栅格的高程值、坡度值、四个方向邻域栅格的坡度值、坡向值、四个方向邻域栅格的坡向值；

以此类推，构建全部250幅dem影像中采样栅格的矩阵l1，l2，...，l250，将其按如下方式进行组合，得到bp神经网络模型的输入矩阵p，

p＝[l1，l2，...，l250]

式中，p为17行，10000列的矩阵；

在本实施例中，构建的bp神经网络模型期望输出矩阵的具体步骤为：

该期望输出矩阵的每一个元素值为该元素所对应栅格的年天文辐射值，是通过数学模型推导计算得到的；

构建编号为1的dem影像中所有采样栅格的年天文辐射矩阵t1，

t1＝[y1，y2，...，y40]

式中，t1为1行40列的矩阵，下标表示dem影像的编号；y1，y2，...，y40分别为上述所有采样栅格的年天文辐射值；

以此类推，构建全部250幅dem影像所有采样栅格的天文辐射矩阵t1，t2，...，t250，将其按如下方式进行组合，得到bp神经网络模型的输出矩阵t，

t＝[t1，t2，...，t250]

式中，t为1行10000列的矩阵；

需保存期望输出矩阵归一化前所有元素中的最大值和最小值，记为ymax和ymin；

第三步：将上述第二步中的输入矩阵和期望输出矩阵分别归一化；

在一个优选的实施例中，使用matlab软件，调用归一化函数mapminmax，以期望输出矩阵的归一化为例，调用方式为：

[tt，tn]＝mapminmax(t，tmax，tmin)

式中，tt为期望输出矩阵归一化后的矩阵，t为期望输出矩阵归一化前的矩阵；tn为归一化时所获得的参数，它主要包含了期望输出矩阵归一化前所有元素的最大值和最小值，以及归一化后的矩阵tt所有元素的最大值tmax和最小值tmin(tmin默认为-1，tmax默认为1)。

第四步：构建bp神经网络模型的框架，设置网络基本参数；

在本实施例中，设置网络的基本参数具体为：最大训练次数为10000；学习率设置为0.01；最大失败次数设置为20；

第五步：初始化各层权值矩阵和阈值向量；

在一个优选的实施例中，可以借助计算机编程语言辅助实现初始化bp神经网络各层的权值矩阵和阈值向量，比如：matlab软件的rands函数。

第六步：将归一化后的输入矩阵，输入bp神经网络模型的输入层，计算出其输出矩阵，求出输出矩阵与归一化后的期望输出矩阵之间的训练误差e；

在本实施例中，所述训练误差的公式为：

式中，n为训练样本即栅格的数量；τk为归一化后的期望输出矩阵第k个元素的值；ok为bp神经网络模型的输出矩阵第k个元素的值。

第七步：判断训练误差e是否达到设定的要求，本实施例中，将训练最小误差e设置为100，若e<e，则结束训练，并转向第八步；反之，则反向修正bp神经网络模型的各层权值矩阵和阈值向量，再转向第六步；

第八步：保存各层权值矩阵和阈值向量，bp神经网络模型构建完成；

步骤5，对输出矩阵进行反归一化，转换成所述计算地区的年天文辐射数值矩阵。

在本实施例中，所述“对输出矩阵进行反归一化，转换成所述计算地区的年天文辐射数值矩阵”，包括以下步骤：

步骤51：对bp神经网络模型的输出矩阵进行反归一化处理；

在本实施例中，所述的对输出矩阵进行反归一化的公式为：

式中：j＝1，2，3，...，12967201；outputj为反归一化后的值；oj为bp神经网络模型的输出值；ymax和ymin是期望输出矩阵归一化前所有元素的最大值和最小值。

在一个优选的实施例中，使用matlab软件进行反归一化，调用mapminmax函数完成，函数的调用方式如下：

output＝mapminmax(‘reverse’，o，tn)

式中，output为反归一化后的矩阵，o为归一化前的矩阵；tn为期望输出矩阵归一化时所获得的参数。

步骤52：将反归一化后的单行输出矩阵重新组合为p1行p2列的矩阵，组合时每p1个元素为一列，即第1到第p1个元素为第一列，第p1+1到第2*p1个元素为第二列，以此类推，组合成一个p1行p2列的矩阵。

在本实施例中，将反归一化后的单行矩阵output重新组合为3601行3601列的矩阵，组合时每3601个元素为一列，即第1到第3601个元素为第一列，第3602到第7202个元素为第二列，以此类推，组合成一个3601行3601列的矩阵。

实施例1：

假设需计算北京市顺义区某区域的年天文辐射。该区域的经纬度范围约为40°n～40°18＇n，116°28＇e～116°58＇e，如图5为本发明实施例中包含了北京市顺义区的单幅dem影像图。

具体计算步骤如下：

1)根据该区域的经纬度范围，下载包含该区域的dem影像。由于dem影像的经纬度范围为40°n～41°n，116°e～117°e，跨度仅为1°，为方便计算，各栅格的经纬度坐标均取dem影像中心栅格的经纬度坐标(40°30′n，116°30′e)，转换为弧度制为(0.7069n，2.0333e)。将dem影像转换为数值矩阵；随后需对dem影像中存在高程空缺值的栅格进行填充(此处使用邻近法进行填充，即将空缺值处填入最邻近栅格的高程值)。

2)如步骤二所述，提取输入矩阵input17×12967201；

3)如步骤三所述，对input17×12967201进行归一化，得到p17×12967201矩阵；

4)如步骤四所述，将p17×12967201输入图2所示bp神经网络模型，计算出输出矩阵o1×12967201。

5)如步骤五所述，对输出矩阵o1×12967201进行反归一化得到output1×12967201。将单行output1×12967201矩阵重新组合为3601×3601的矩阵，得到该样区的年天文辐射矩阵output3601×3601。如图6为本发明实施例中包含了北京市顺义区的年天文辐射量分布图。

上述具体实施方式不能认为是对本发明的进一步限定，本领域技术人员根据本发明内容作出的非实质性改变均应落入本发明保护范围内。