一种基于异质格网高程线性分解模型的高程值解算方法

1.本发明属于地表高程数据处理技术领域，尤其涉及一种基于异质格网高程线性分解模型的高程值解算方法。


背景技术：

2.常见的dsm是规则正方形格网结构，单位格网范围内水平方向上高低错落的裸露地面及地表覆盖物，由一个高程值统一概括会导致“削峰去谷”以及地形特征位移等一系列的误差，格网越大误差越明显，这种现象称为格网内部异质性。目前dsm水平分辨率最高为30m，而人类活动频繁的区域下垫面情况复杂，包含大量长度或宽度小于30m的地表覆盖物，相比其他区域更加容易出现同一个格网包含多种地表，并用一个高程值概括描述的情况；因此，30m分辨率的dsm无法表达密集、多变的地表情况，在地表覆盖密集区或存在相邻格网高程值相差很小，甚至完全一样的情况；dsm中的内部异质格网类似与遥感影像中的混合像元，遥感影像中由于混合像元问题，像元级遥感分类和面积测量的精度难以满足实际要求。人类活动频繁地区的下垫面情况复杂，需要在高分辨率高精度地表模型的基础上进行分析研究。目前大范围公开的免费数据无法反映该类地表的细微形态，而高精度数据获取成本高、难度大；即便有高分辨率数据，由此自动获取的dsm精度也难以保证。
3.为解决上述问题，本发明提出了异质性格网的高程线性分解模型，可以提高低分辨率dsm应用的精度,使dsm应用由格网级达到亚格网级,进入格网内部，在只有低分辨率dsm的情况下，为获取高分辨率dsm提供新思路。


技术实现要素：

4.本发明的目的在于使dsm应用由格网级达到亚格网级,进入格网内部，提供了一种基于异质格网高程线性分解模型的高程值解算方法，解决了目前低分辨率dsm应用精度低的问题。
5.本发明的目的通过以下技术方案来实现：一种基于异质格网高程线性分解模型的高程值解算方法，所述高程值解算方法包括：
6.步骤1：基于辅助数据(遥感影像、众源数据等)提取研究区域的地表覆盖分类矢量数据；
7.步骤2：提取待分解低分辨率的dsm矢量格式的格网轮廓，将其与地表覆盖分类矢量数据利用地理信息叠置分析中的联合计算方法计算出dsm每个格网中亚格网的面积比例；
8.步骤3：根据异质格网综合高程值和格网内各亚格网高程值的线性函数关系，计算剔除非连续亚格网后只包含连续亚格网的格网综合高程值；
9.步骤4：假设相邻格网(8个)中的同一种连续亚格网高程值相同，结合只包含连续亚格网的综合高程值和连续亚格网面积比例，列出待求格网和相邻格网中连续亚格网高程值和只包含连续亚格网的综合高程值的关系方程，共9个；
10.步骤5：根据多元线性回归的方法解算方程组，分解待求格网的连续亚格网综合高程值，得到各连续亚格网的高程值。
11.所述高程值解算方法还包括：以阈值分辨率dsm作为验证数据源，采用各亚格网轮廓作为掩膜对验证数据源的高程值通过均值法进行分区统计，将得到的值作为观察研究区域的真实高程值并对解算得到的亚格网高程值进行误差统计，得到异质格网高程线性分解模型的解算精度。
12.所述根据异质格网综合高程值和格网内各亚格网高程值的函数关系，计算剔除非连续亚格网后只包含连续亚格网的格网综合高程值包括：
13.将异质格网的综合高程值设置为i＝1,2,3,
…
,n，ei表示格网内第i个亚格网的高程值，ri表示第i个亚格网所占格网的面积比例，ei表示第i个亚格网高程值的误差，n表示dsm中选定的亚格网个数；
14.根据地理空间上距离越近相似度越高的规律将距离近的m个格网中相同类型的连续亚格网高程值设置为相等，得到m个相邻格网中个亚格网的高程值、面积比例与各格网的平均高程值之间关系的方程组：
[0015][0016]
其中，eci(i＝1，2，
…
，n)表示m个格网内第i类表面连续的亚格网的高程值，rc
ji
(j＝1，2，
…
，m；i＝1，2，
…
，n)表示第j个格网中第i类表面连续的亚格网所占格网的面积比例，ed
jk
(j＝1，2，
…
，m；k＝1，2，
…
，l)表示第j个格网中第k类表面非连续表面的亚格网的高程值，rd
jk
(j＝1，2，
…
，m；k＝1，2，
…
，l)表示第j个格网中第k类表面非连续表面的亚格网所占格网的面积比例，vj(j＝1，2，
…
，m)表示第j个格网的综合高程值；
[0017]
得到第j个格网中非连续地表的亚格网线性加权高程值：vdj＝∑ed
jk
rd
jk
(j＝1，2，
…
，m；k＝1，2，
…
，l)，以及第j个格网中连续地表的亚格网线性加权高程值：vcj＝v
j-vdj＝v
j-∑ed
jk
rd
jk
(j＝1，2，
…
，m；k＝1，2，
…
，l)，结合方程组得到：
[0018]
对进行变换处理得到异质格网高程线性分解模型的系数估计公式ec＝(rc
t
rc)-1
rc
t
vc；
[0019]
根据得到的第j个格网中连续地表的亚格网线性加权高程值计算剔除了第j个格网中非连续亚格网线性加权高程值后的第j个格网的综合高程值：vbj＝vcj/∑rc
jk
＝v
j-∑ed
jk
rd
jk
/∑rc
jk
，(j＝1，2，
…
，m；k＝1，2，
…
，l)，进一步变换得到：
[0020]
所述对进行变换处理得到异质格网高程线性分解模型的系数估计公式ec＝(rc
t
rc)-1
rc
t
vc包括：
[0021]
将矩阵和向量化得到矩阵向量和
[0022]
进而转换为矩阵形式ecrc＝vc，在ecrc＝vc两边同时左乘rc的转置矩阵rc
t
得到rc
t
rcec＝rc
t
vc；
[0023]
在rc
t
rcec＝rc
t
vc两边同时左乘rc
t
rc方阵的逆矩阵(rc
t
rc)-1
得到
[0024]
得到ecrc＝vc，对ecrc＝vc两边进行转置矩阵和逆矩阵变换处理得到异质格网高程线性分解模型的系数估计公式ec＝(rc
t
rc)-1
rc
t
vc。
[0025]
异质格网高程线性分解模型中每个格网的综合高程值表示为它的亚格网高程值与其所占格网面积比例的线性组合，且每个异质格网内同一个亚格网的高程值是该亚格网各处高程的平均值；各亚格网都是内部连续，相互之间不连续，将面积小能被一个格网覆盖的亚格网设置为非连续地表，将面积大且跨越多个格网的亚格网设置为连续地表。
[0026]
本发明具有以下优点：一种基于异质格网高程线性分解模型的高程值解算方法，提出了异质格网综合高程值的线性分解模型，选取研究区运用该模型进行高程分解获得亚格网高程值，将分解高程值与验证数据进行对比，使dsm应用由格网级达到亚格网级，进入格网内部，提高了低分辨率dsm的应用精度，适用于非连续地表覆盖物(如房屋、植被等)分布较少或者非连续地表覆盖物高度值易获取的地区。
附图说明
[0027]
图1为本发明方法的流程示意图；
[0028]
图2为dsm格网高程分解试验区域范围图；
[0029]
图3为高程分解实验区域地表覆盖分类图；
[0030]
图4为异质格网分解得到的亚格网高程图；
[0031]
图5为0.1m dsm分区统计均值得到的亚格网高程图；
[0032]
图6为高程分解值绝对误差与亚格网面积占比的关系图。
具体实施方式
[0033]
为使本技术实施例的目的、技术方案和优点更加清楚，下面将结合本技术实施例中附图，对本技术实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本技术一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本技术实施例的组件可以以各种不同的配置来布置和设计。因此，以下结合附图中提供的本技术的实施例的详细描述并非旨在限制要求保护的本技术的保护范围，而是仅仅表示本技术的选定实施例。基于本技术的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得
的所有其他实施例，都属于本技术保护的范围。下面结合附图对本发明做进一步的描述。
[0034]
通过实地测量、立体像对提取、激光扫描测量等技术手段获得的高程信息通常是以规则网格的形式存储在格网dsm(数字表面模型)中，而每个格网内记录的信息代表了该格网范围内的地面高程，由格网内或附近某一点的高程值代替，也可以通过周围采样点的高程值插值得到，是一个估计值。但是单个格网中有可能包含了多个地形或地表覆盖物，如悬崖、建筑物等突变地形周围的高程值呈不连续急剧变化，相差较大，由于规则格网dsm的结构局限性，只能用一个估计的高程值来概括剧烈起伏的地表，这种格网称为内部异质性格网，异质性格网的值，无论是取自某一采样点，还是其他插值方法得到的估计值，都会造成地形表达失真。若简单地用该高程值表示异质性格网内的某一地物高程，显然是不合理的。
[0035]
如果研究的dsm格网内部仅仅含有一种连续形态，那么这个格网就被认为是同质格网，否则就称该格网为异质格网。同质格网的高程值就是所含唯一地形的高程值，异质格网的高程值是多种地形高程的综合值，称为异质格网的综合高程值。异质格网内部的某一类地形或地表覆盖物也称其为亚格网。亚格网尺度小于一个格网，位于格网内部，一般情况下，用来分析应用的dsm最基本的单位是格网。dsm中的异质格网实际包含了不连续的多种地表形态，然而各地表形态通常不能直接由格网的综合高程值表示，于是就需要进行异质格网的分解，再根据各地表形态对于某种应用的重要程度或者在格网内的面积比例等因素，通过调整高程值对各地表形态进行强调或削弱。异质格网分解的主要目的是得到更加合理的高程值，使地形的表达更加符合实际或便于应用。异质格网通常出现在地表形态破碎的地方，多种不连续地表的交界处，以及下垫面情况复杂的区域(如人类活动频繁区域)。dsm的格网异质性随分辨率增高而减弱。
[0036]
如图1所示，本发明包括的内容如下：
[0037]
步骤1：用0.1m分辨率dsm均值重采样成分辨率为30m的低分辨率dsm，选取30个格网作为异质格网高程分解实验数据。
[0038]
步骤2：利用辅助数据(高分辨率正射航空影像)提取研究区的地表覆盖分类情况(图3)。
[0039]
步骤3：提取待分解低分辨率dsm矢量格式的格网轮廓，将其与地表覆盖分类矢量数据利用地理信息叠置分析中的联合计算方法计算出30个dsm格网中亚格网的面积比例。
[0040]
步骤4：对于非连续地表覆盖物——建筑物和植被，根据异质格网综合高程值和格网内各亚格网高程值的线性关系，利用公式(2)-(5)，计算剔除房屋类建筑物和植被(非连续亚格网)后的只包含裸地和连续地表覆盖物(连续亚格网)的格网综合高程值。将该综合值用于接下来的地表覆盖物高程值解算。
[0041]
步骤5：假设相邻格网(8个)中的同一种连续亚格网高程值相同，结合只包含连续亚格网的综合高程值和连续亚格网面积比例，列出待求格网和相邻格网中连续亚格网高程值和只包含连续亚格网的综合高程值的关系方程，共9个。
[0042]
步骤6；根据多元线性回归的方法解算方程组(6)-(11)，分解待求格网的连续亚格网综合高程值，得到各连续亚格网的高程值。
[0043]
步骤7：把0.1m分辨率dsm作为验证数据源，采用各亚格网轮廓作为掩膜对验证数据源的高程值采用均值法分区统计，将得到的值作为地块的真实高程值，对解算得到的亚
格网高程值进行误差统计，从而获得该线性分解模型的解算精度。
[0044]
进行dsm异质格网高程分解首先就是要提取其所包含的亚格网，亚格网是异质格网中代表某种连续地形的同质格网，是dsm中组成异质格网的基本单元。异质格网亚格网的提取，一是要遵循地形突变原则；二是要结合数字地形模型的实际应用，通过格网中各地形对实际应用的影响程度(本实例使用的是面积权重)，确定参与赋值计算的地表高程。
[0045]
由高分辨率的高程数据均值重采样得到的低分辨率dsm，其异质格网的综合高程值是亚格网表面高程值与面积比例的乘积之和。因此，均值重采样dsm的异质格网高程分解模型为线性模型。该模型把dsm异质格网的高程值表示为各亚格网高程及其所占格网面积比例的线性组合。异质格网高程线性分解模型是在假设每个异质格网内同一个亚格网的高程值是该亚格网各处高程的平均值的基础上建立的，优点是构模简单，其物理含义明确，理论上有较好的科学性，对于解决dsm格网异质性问题有一定的效果。
[0046]
dsm异质格网的高程线性分解模型中，每一个格网的综合高程值表示为它的亚格网高程值与它们各自所占格网的面积比例的线性组合。因此，异质格网综合高程值可以表示为：
[0047][0048]
其中，i＝1,2,3,
…
,n，ei表示格网内第i个亚格网的高程值，ri表示第i个亚格网所占格网的面积比例，ei表示第i个亚格网高程值的误差，n表示dsm中选定的亚格网个数；
[0049]
各亚格网都是内部连续，相互之间不连续。本发明将面积相对较小，通常能被一个格网覆盖的亚格网(如房屋等建筑物)视为非连续地表；将面积相对较大，通常跨越多个格网的亚格网(如河道、道路等)视为连续地表。根据地理空间上距离越近相似性越高的规律，可以假设距离较近的m个格网中相同类型的连续亚格网高程值相等。则m个相邻格网中各亚格网的高程值、面积比例与各格网的平均高程值之间的关系可以表示为以下方程组：
[0050][0051]
其中，eci(i＝1，2，
…
，n)表示m个格网内第i类表面连续的亚格网的高程值，rc
ji
(j＝1，2，
…
，m；i＝1，2，
…
，n)表示第j个格网中第i类表面连续的亚格网所占格网的面积比例，ed
jk
(j＝1，2，
…
，m；k＝1，2，
…
，l)表示第j个格网中第k类表面非连续表面的亚格网的高程值，rd
jk
(j＝1，2，
…
，m；l＝1，2，
…
，l)表示第j个格网中第k类表面非连续表面的亚格网所占格网的面积比例，vj(j＝1，2，
…
，m)表示第j个格网的综合高程值；
[0052]
对于第j个格网中非连续地表的亚格网的线性加权高程可以表示为：
[0053]
vdj＝∑ed
jk
rd
jk
(j＝1，2，
…
，m；k＝1，2，
…
，l)
ꢀꢀꢀ
(3)
[0054]
因此，第j个格网中连续地表的亚格网线性加权高程值表示为：
[0055]
vcj＝v
j-vdj＝v
j-∑ed
jk
rd
jk
(j＝1，2，
…
，m；k＝1，2，
…
，l)
ꢀꢀꢀ
(4)
[0056]
从而得到：
[0057][0058]
公式(5)为一个不含常数项的多元线性方程组，引入以下矩阵和向量：
[0059][0060]
将公式(5)变换为矩阵形式：
[0061]
ecrc＝vc
ꢀꢀꢀ
(7)
[0062]
可以用多元线性回归模型求解该方程组的最小二乘解。为了求解系数向量ec，需要开展一些过渡性的数学变换。在公式(7)两边同时左乘rc的转置矩阵rc
t
得到：
[0063]
rc
t
rcec＝rc
t
vc
ꢀꢀꢀ
(8)
[0064]
这样式中出现方阵rc
t
rc，在(3-8)两边同时左乘这个方阵的逆矩阵(rc
t
rc)-1
，得到：
[0065]
ec＝(rc
t
rc)-1
rc
t
vc
ꢀꢀꢀ
(9)
[0066]
这就得到了多元线性回归模型的系数估计公式。
[0067]
求得代表连续地表的亚格网的高程值后，根据公式(2)-(5)计算剔除了非连续亚格网后的第j个格网的综合高程值：
[0068]
vbj＝vcj/∑rc
jk
＝v
j-∑ed
jk
rd
jk
/∑rc
jk
，(j＝1，2，
…
，m；k＝1，2，
…
，l)
ꢀꢀꢀ
(10)
[0069]
进一步变换得到：
[0070][0071]
本发明选取了澳洲墨尔本雅拉河畔面积约为22622m2的地块作为研究区，区内有高架桥、铁路等地表覆盖物和裸露的地表，在周围区域还有少量建筑物和植被。采用的高程数据是0.1m分辨率dsm，辅助数据是分辨率为0.075m的正射航空影像，来源于墨尔本政府网站(https://data.melbourne.vic.gov.au/)。数据详情见表1：
[0072]
表1研究区域数据信息表
[0073] 0.1m分辨率dsm正射航空影像数据日期2018.52018.5密度/分辨率0.100m0.075m空间参考wgs_1984_utm_zone_55sgda94_mga_zone_55基准面d_wgs_1984d_gda_1994
[0074]
该实验将研究区范围内的0.1m分辨率dsm均值重采样得到30个约30m(1
″
)分辨率的dsm格网，再运用文中提出的异质格网高程线性分解模型分解出异质格网的亚格网高程值，最后，将各亚格网轮廓作为掩膜对0.1m分辨率dsm采用均值法分区统计，认为得到的值是各亚格网的真实高程值，并将其作为验证数据对分解得到的亚格网高程值进行误差统计，从而获得该线性分解模型的解算精度。图2中，实验区内的灰色格网是参与计算的格网，黑色格网既是参与计算的格网也是待分解格网。本实验将房屋类建筑物和植被作为非连续地表，高架桥、铁路和裸地等作为连续地表，采用基于3
×
3的模板进行高程值分解运算。
[0075]
如图3-图5所示，研究区共有包括裸地在内的10类地表，其中高架桥凌驾于铁路、裸地以及小路1上方，与小路2相接。实验区内的非连续地表覆盖物——房屋类建筑物和植被都位于周边格网。位于中间区域待分解的12个格网，包含高架桥、铁路2、小路1三种地表
覆盖物，经过格网划分，共有24个地表覆盖物高程值分解地块。12个格网中有3个同质格网，对应的是地表覆盖物是铁路2和高架桥，则这3个格网中的地表覆盖高程值就是格网的综合高程值；剩下9个格网经过分解计算，得到各自亚格网的高程值。
[0076]
通过高程分解模型计算的地表覆盖高程值和验证高程值对比如下：验证高程值在高架桥上的高程值出现明显的高程差，这主要是由地表覆盖分类误差造成的；而分解高程值在计算过程中，会根据地表类型设定来统一计算，对于分类产生的误差会均摊到周围8个格网中，因此相对来说分解高程值在同一地表覆盖物上的高程值更加连续。
[0077]
表2异质格网高程分解结果与验证数据对比信息表
[0078][0079]
表3高程分解值误差信息表
[0080]
[0081][0082]
由表2和表3发现，24个分解高程值与验证值整体接近，其绝对误差平均值约为1.479m，分解值基于验证值的分散度为2.940。绝对误差小于1m的值为18个，百分比含量83.3％；绝对值误差小于2m的值有20个，百分比含量75％。其中绝对误差值大于2的4个地块，其在格网内部面积占比都小于0.03。图6呈现了高程分解值的绝对误差与地块面积占比的关系曲线，可知随着地块面积占比增大，运用线性分解模型得到的高程值的精度越高。
[0083]
本发明针对介绍了规则格网dsm的格网异质性，将异质格网类比于遥感影像中的混合像元，提出了异质格网中的亚格网概念；提出了异质格网综合高程值的线性分解模型，选取研究区运用该模型进行高程分解获得亚格网高程值，将分解高程值与验证数据进行对比。结果证明，异质格网的高程值线性分解法可以较好地反解出剔除了非连续地表覆盖物后的格网中各类地表覆盖高程值，提供了一种较低成本下，由低分辨率dsm获取高分辨率dsm的新思路新方法。然而，亚格网的高程值解算高度依赖于待分解源数据dsm的高程精度和地表覆盖物轮廓的准确性，该模型特点会使得数据本身误差在传递过程中放大，因此，异质格网中面积比例小的亚格网高程值的解算结果容易出现较大误差。该方法适用于非连续地表覆盖物(如房屋、植被等)分布较少或者高度值易获取的地区。
[0084]
以上所述仅是本发明的优选实施方式，应当理解本发明并非局限于本发明所披露的形式，不应看作是对其他实施例的排除，而可用于各种其他组合、修改和环境，并能够在本发明所述构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离本发明的精神和范围，则都应在本发明所附权利要求的保护范围内。