一种InSAR地形产品水体高程编辑处理方法

一种insar地形产品水体高程编辑处理方法
技术领域
1.本发明涉及基于insar技术生产dem技术领域，特别涉及一种insar地形产品水体高程编辑处理方法。


背景技术：

[0002]“数字地球”是一个将真实世界以及其相关现象统一用数字化手段重现并进行多分辨率、多尺度、多时空和多种类三维描述，并利用它作为工具来支持和改善人类活动和生活质量。“数字地球”涉及到多种类型的空间数据，其中最基础和最重要的是数字地形数据。常见的数字地形数据有dsm和dem两种。dsm(digital surface model)是指描述包含建筑物、桥梁、树木等所有可见地物高程变化情况的二维数组或栅格图像，它可以反映整个自然或人造地表表面，并为各种三维建模和可视化提供基础数据。dem(digital elevation model)则只包含裸露或去除了其他覆盖物后的自然或人造地面高程变化情况的二维数组或栅格图像，它可以反映真实的地形起伏，并为各种科学研究和工程应用提供基础数据。
[0003]
dsm和dem可以通过多种方式获取，如野外测量、摄影测量、激光扫描和insar等。insar是一种利用微波遥感技术进行全天候、全天时地表观测的方法，它具有高精度、大范围、低成本等优势，是生产全球dem的最佳手段之一。但insar获取的是dsm而不是dem，因此需要从dsm中提取出dem，这是“数字地球”中的一个重要问题。
[0004]
由于镜面反射导致的低后向散射以及时间去相干，水体在由insar技术得到的原始dsm中表现为随机噪声或高程异常，因此需要对水体高程进行编辑处理，以使其在保证高程真实的前提下，达到平滑的效果。


技术实现要素：

[0005]
本发明的目的是：根据不同水体展现出的不同特性以及实际数据的特点，将水体分成海洋、湖泊和河流三类，并采用不同的算法对不同水体进行高程编辑处理，得到水体高程真实且平滑的dem。
[0006]
为了达到上述目的，本发明提供了一种insar地形产品水体高程编辑处理方法，包括如下步骤：
[0007]
s1，对水体掩膜进行数据预处理；
[0008]
s2，利用全球高分辨率海岸线矢量数据，对水体掩膜中的海洋进行自动化分类；
[0009]
s3，根据湖泊和河流的形态面积的不同，对水体掩膜中剩余的湖泊和河流进行自动化分类；
[0010]
s4，对海洋高程进行自动化编辑处理，计算海洋像素的大地水准面起伏高度；
[0011]
s5，对湖泊高程进行自动化编辑处理；对河流高程进行自动化编辑处理；
[0012]
s6，基于水体高程初步编辑结果，自动识别连接在一起的异类水体，并在接边位置进行高程的平滑过渡处理。
[0013]
进一步地，s1中将水体提取结果处理成二值掩膜图；对水体掩膜中的空洞进行填
补处理；对水体掩膜进行包括平滑、去噪和增强形态学处理。
[0014]
进一步地，s2中根据dsm图幅的地理范围确定矩形矢量，通过矩形矢量和全球海岸线矢量数据进行矢量求交运算，得到图幅地理范围内的海岸线面状矢量；将海岸线面状矢量转换成与dsm图幅分辨率一致的栅格文件，得到与dsm图幅的地理范围及分辨率一致的海洋掩膜栅格文件；根据海洋掩膜文件将水体中的海洋区分出来。
[0015]
进一步地，s3中将已分类成海洋的水体在水体掩膜中排除，确定剩余水体连通域的范围；再根据湖泊和河流表现出的一般形态的不同，设计指标和阈值，对剩余水体中的湖泊和河流进行自动化分类，并在水体分类结果图中分别进行标记，最终得到水体分类结果。
[0016]
进一步地，根据分类水域面积与水体最大延伸的比值指标t，指标大于阈值a的连通域分类成湖泊，指标小于阈值a的连通域分类成河流。
[0017]
进一步地，s4中利用egm2008地球重力场模型，依据下式计算大地水准面起伏高度：
[0018][0019]
其中，n
gm
为模型大地水准面高，r，θ，λ分别为地心距离、余纬和经度，gm为引力常数与地球质量的乘积，a为参考椭球长半轴，和为完全规格化系数，为完全规格化缔合legendre函数，n
max
为地球重力场模型展开的最高阶数；
[0020]
其中，完全规格化缔合legendre函数采用跨阶次递推法进行计算。
[0021]
进一步地，s5中根据连通域膨胀原理获得湖泊岸边像素的高程值，绘制出分布直方图，选择直方图中的峰值下降的四分之一处的高程值，作为湖泊的统一高程。
[0022]
进一步地，s5中根据河流两岸像素的高程值确定河流高程，且河流高程由高到低平滑变化；对河流两岸像素的高程值进行去噪处理；根据去噪后的河流两岸像素的高程值进行自动化插值；对河流两岸像素的高程值进行多次局部平滑滤波处理。
[0023]
进一步地，运用克里金插值法对河流两岸像素的高程值进行插值，克里金插值法的公式定义为：
[0024][0025]
其中，是点(x0,y0)的估值，即z0＝(x0,y0)；λi是权重系数，用空间上所有已知点的数据加权求和来估计未知点的值，满足估值与真值之差最小以及无偏估计条件的一套最优系数。
[0026]
进一步地，s6中根据水体分类结果，识别出连接在一起的异类水体，并找出接边位置；基于水体高程初步编辑结果，在异类水体的接边位置处进行高程的平滑过渡处理，消除高程跳变现象。
[0027]
本发明的上述方案有如下的有益效果：
[0028]
本发明提供的insar地形产品水体高程编辑处理方法，充分应用了不同类别的水体展现出的不同特性以及实际数据的特点，对存在噪声的insardsm中的水体进行高程编辑
处理，在保证水体高程真实的前提下，实现各类水体的平滑处理以及异类水体接边的平滑过渡，整个流程结构清晰，具有流程清晰、实现简单、自动化程度高等优点；
[0029]
本发明的其它有益效果将在随后的具体实施方式部分予以详细说明。
附图说明
[0030]
图1是本发明的流程框图；
[0031]
图2是本发明实施例中的水体掩膜示意图；
[0032]
图3是本发明实施例中的水体分类结果示意图；
[0033]
图4是本发明实施例中由insar技术得到的原始dsm示意图；
[0034]
图5是本发明实施例中水体高程编辑结果示意图。
具体实施方式
[0035]
以下通过特定的具体实例说明本公开的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本公开的其他优点与功效。显然，所描述的实施例仅仅是本公开一部分实施例，而不是全部的实施例。本公开还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本公开的精神下进行各种修饰或改变。需说明的是，在不冲突的情况下，以下实施例及实施例中的特征可以相互组合。基于本公开中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本公开保护的范围。
[0036]
需要说明的是，下文描述在所附权利要求书的范围内的实施例的各种方面。应显而易见，本文中所描述的方面可体现于广泛多种形式中，且本文中所描述的任何特定结构及/或功能仅为说明性的。基于本公开，所属领域的技术人员应了解，本文中所描述的一个方面可与任何其它方面独立地实施，且可以各种方式组合这些方面中的两者或两者以上。举例来说，可使用本文中所阐述的任何数目个方面来实施设备及/或实践方法。另外，可使用除了本文中所阐述的方面中的一或多者之外的其它结构及/或功能性实施此设备及/或实践此方法。
[0037]
还需要说明的是，以下实施例中所提供的图示仅以示意方式说明本公开的基本构想，图式中仅显示与本公开中有关的组件而非按照实际实施时的组件数目、形状及尺寸绘制，其实际实施时各组件的型态、数量及比例可为一种随意的改变，且其组件布局型态也可能更为复杂。另外，在以下描述中，提供具体细节是为了便于透彻理解实例。然而，所属领域的技术人员将理解，可在没有这些特定细节的情况下实践所述方面。
[0038]
如图1所示，本发明的实施例提供了一种insar地形产品水体高程编辑处理方法，包括如下步骤：
[0039]
s1，对水体掩膜进行数据预处理。
[0040]
在本实施例中，先将水体提取结果处理成二值掩膜图；对水体掩膜中的空洞(像素值为空)进行填补处理；对水体掩膜进行包括平滑、去噪和增强等形态学处理，以提高水体分类精度，完成对水体掩膜的数据预处理。
[0041]
s2，利用全球高分辨率海岸线矢量数据，对水体掩膜中的海洋进行自动化分类。
[0042]
在本实施例中，先根据水体掩膜文件，确定各个水体连通域的范围；在将全球高分
辨率海岸线矢量数据处理成与待处理的原始dsm图幅的地理范围以及分辨率一致的栅格文件，根据栅格文件自动化识别水体掩膜中的海洋，并将其在水体分类结果图中进行标记。
[0043]
需要说明的是，海洋和湖泊在形态面积上较为接近，很难区分，因此需要引入外部数据加以辅助分类。目前公开发布的关于水域的产品中，虽然有湖泊和河流的水域掩膜，但是这些水系掩膜往往并未分类成多种水体，数据更新周期也较长，或者只存在某个地区的水域掩膜数据，因此难以适应全球dem产品生产的需求。但是全球海岸线的数据产品是存在的，虽然海岸线也经常随着时间有所变化，但是相对于整个海洋的面积，这些变化可以忽略不计，而且在实际海洋分类时，只利用了海岸线数据进行辅助分类，并不影响实际水体提取的精度，其导致海洋错误分类的影响也很小。所以在实际海洋分类时，采用了全球高分辨率海岸线矢量数据，首先对水体掩膜中的海洋进行了自动化分类。
[0044]
具体到本实施例中，将全球高分辨率海岸线矢量数据裁剪成与待处理dsm图幅的地理范围以及分辨率一致的栅格文件：首先根据dsm图幅的地理范围确定一个矩形矢量，然后用这个矩形矢量和全球海岸线矢量数据进行矢量求交运算，得到图幅地理范围内的海岸线面状矢量，然后将面状矢量转换成与dsm图幅分辨率一致的栅格文件，这样便可得到与dsm图幅的地理范围及分辨率一致的海洋掩膜栅格文件，再根据海洋掩膜文件将水体中的海洋分类出来即可。
[0045]
s3，根据湖泊和河流的形态面积的不同，对水体掩膜中剩余的湖泊和河流进行自动化分类。
[0046]
在本实施例中，先将已分类成海洋的水体在水体掩膜中排除，确定剩余水体连通域的范围；再根据湖泊和河流表现出的一般形态的不同，即湖泊的形态一般接近圆形，河流的形态一般是细长型，设计相应的指标和阈值，据此对剩余水体中的湖泊和河流进行自动化分类，并在水体分类结果图中分别进行标记，最终得到水体分类结果。
[0047]
需要说明的是，湖泊是静止水体，河流是流动水体，理论上可以从高程变化幅度的角度出发，对湖泊和河流进行分类。但是由于镜面反射导致的低后向散射以及时间去相干，水体在由insar技术得到的原始dsm中表现为随机噪声或高程异常，无法反映真实的水体高程，所以很难从高程上分辨出湖泊和河流。因此在实际分类中，主要是根据湖泊和河流的形态面积的不同而进行分类，据此本实施例设置了相应的分类指标。因为湖泊的形态一般接近圆形，河流的形态一般是细长型的，所以可以用水域面积与水体最大延伸的比值这一指标对湖泊和河流进行分类，再根据大量实际数据总结出来合适的阈值，从而将除海洋外的剩余水体掩膜自动化分类成湖泊和河流。
[0048]
具体地，根据分类水域面积与水体最大延伸的比值指标t，指标大于阈值a(即t＞a)的连通域，分类成湖泊，反之分类成河流。
[0049]
s4，利用egm2008地球重力场模型，对海洋高程进行自动化编辑处理，计算海洋像素的大地水准面起伏高度。
[0050]
在本实施例中，考虑到真实海水面的高度无法直接测得，先对海洋像素赋值大地水准面的起伏高度，以表示静止状态下的平均海水面的高度；再利用egm2008地球重力场模型自动化计算海洋像素的大地水准面起伏高度；计算大地水准面起伏高度的关键在于计算超高阶缔合legendre函数，本实施例在计算超高阶缔合legendre函数时，采用了相较于传统计算方法精度更高且稳定性更强的跨阶次递推方法。
[0051]
需要说明的是，大地水准面起伏高度也称为大地水准面差距，是指大地水准面上的某一点，沿法线投影到参考椭球面上的相应位置之间的距离。大地准面是指地球重力场中，与处于自由静止状态的平均海水面相重合或最为接近的重力等位面，利用重力数据计算的大地水准面，称为重力大地水准面。重力大地水准面起伏高度常利用地球重力场模型进行计算。实际计算时采用了egm2008地球重力场模型进行计算，egm2008(earth gravitational model 2008)是地球重力场模型的简称，它是由美国国家地理空间情报局(nga)egm开发团队公开发布的官方版地球重力场模型。该模型的球谐展开阶数达到了2159，并且还提供扩展到2190阶的扩展参数，它结合了精度高、覆盖面广的地面重力数据，全球5
′×5′
网格大地水准面的估算精度为：误差最小为3.045cm，最大为102.194cm，标准差为
±
13cm，因此egm2008在全球具有很高的精度。
[0052]
在计算大地水准面起伏高度时，超高阶缔合legendre函数的计算最为关键，本实施例采用了跨阶次递推算法对超高阶缔合legendre函数进行计算，相比于传统的标准向前列递推法计算法，不仅计算精度更高，而且在超高阶次(3060阶甚至更高阶)缔合legendre函数的计算时，稳定性也更强。因此，可以为海洋高程计算得到高精度的大地水准面起伏高度。
[0053]
具体到本实施例中，利用egm2008地球重力场模型，依据下式计算大地水准面起伏高度：
[0054][0055]
其中，n
gm
为模型大地水准面高，r，θ，λ分别为地心距离、余纬和经度，gm为引力常数与地球质量的乘积，a为参考椭球长半轴，和为完全规格化系数，为完全规格化缔合legendre函数，n
max
为地球重力场模型展开的最高阶数。
[0056]
上式中的完全规格化缔合legendre函数是采用跨阶次递推法进行计算，原理如下：
[0057]
中n是阶数，m是次数，它的初始值如下所示：
[0058][0059]
当m＝0,1时，采用基本递推公式计算，即下式：
[0060][0061]
上式中
[0062]
[0063][0064]
当m≥2时，采用跨阶次递推公式，即下式：
[0065][0066]
上式中：
[0067][0068][0069][0070]
当m≥2时，的值可由线性组合得到，且递推系数α
nm
、β
nm
、γ
nm
的值均小于1，所以跨阶次递推是稳定可靠的。
[0071]
s5，根据统计分析原理，对湖泊高程进行自动化编辑处理；根据实际数据特点，对河流高程进行自动化编辑处理。
[0072]
在本实施例中，对于水体分类得到的湖泊，视其为静止水体，为其赋值统一高程；对湖泊岸边像素的高程进行自动化统计分析，根据统计分析结果为湖泊像素赋值合适的统一高程。
[0073]
需要说明的是，合适的统一高程是通过对湖泊岸边像素的高程值进行统计分析计算得到的。在对不同湖泊进行实证分析后，发现由于山脉、悬崖和水体噪声的存在，湖岸部分像素的高程值过高或者过低，呈现出异常状态，这些像素的高程不适合作为湖泊的统一高程，并且这些像素在整个湖岸像素的中的占比并不高，而那些高程值相对正常的湖岸像素的占比较高，这就导致湖岸像素高程分布直方图呈现出单峰分布的特点。因此，通过对湖泊岸边像素高程值的分布直方图进行自动化分析，即可为湖泊确定合适的统一高程。
[0074]
具体到本实施例中，根据连通域膨胀原理，获得湖泊岸边像素的高程值，绘制出其分布直方图，选择直方图中的峰值下降(朝向较低的高程值)的四分之一处的高程值，作为湖泊的统一高程。
[0075]
对于水体分类得到的河流，视其为流动水体，根据河流两岸像素的高程值确定河流高程，且要求高程由高到低平滑变化；考虑到河流两岸常有植被存在，导致由insar技术获得的河流两岸像素高程中存在着许多噪声和异常值，所以先对其进行去噪处理；根据去噪后的河岸像素高程对河流两岸像素进行自动化高程插值；考虑到插值后的河流两岸像素高程中仍存在噪声，最后对河流两岸像素高程进行多次局部平滑滤波处理。
[0076]
需要说明的是，河流视为流动水体，其高程应是从高到低平滑变化的。为了实现河流的平滑流动，最直接的办法就是根据河流两岸像素的高程对其进行插值。但是考虑到河流两岸往往存在着许多树木、草木等植被，根据insar成像原理，植被区域的相干性较差，所
以在由insar技术得到的原始dem中，河流两岸常常表现出很多噪声或者高程异常值，这对河流插值效果影响很大。因此为了改善河流插值效果，在河流插值前，需要对河流两岸像素的高程进行去噪处理。然后根据去噪后的河流两岸像素的高程值，运用克里金插值法对河流两岸像素进行高程插值。最后，考虑到插值得到的河流两岸像素高程中仍存在部分噪声，所以对河流区域进行多次局部平滑滤波，以达到河流平滑流动的效果。
[0077]
其中，克里金模型是1951年由南非矿业工程师d.g.krige和统计学家h.s.sichel等人提出的用来预测未知点上响应的一种无偏估计模型。根据均值的不同，克里金插值法可以分为简单克里金法(simple kringing，sk)、普通克里金法(ordinary kriging，ok)。普通克里金法使用较广泛，本实施例以普通克里金法进行计算。
[0078]
克里金插值法的公式可以定义为：
[0079][0080]
其中，是点(x0,y0)的估值，即z0＝(x0,y0)；λi是权重系数，用空间上所有已知点的数据加权求和来估计未知点的值，但并非是距离的倒数，而是满足估值与真值之差最小以及无偏估计条件的一套最优系数：
[0081][0082][0083]
其中，式(12)表示估值与真值之差最小条件；式(13)表示无偏估计条件。假设空间属性z是均一的，即对于空间上任意一点，都具有同样的期望与方差：
[0084]
e[z(x,y)]＝e[z]＝c
ꢀꢀꢀ
(14)
[0085]
var[z(x,y)]＝σ2ꢀꢀꢀ
(15)
[0086]
因此，空间上任意一点的值z(x,y)都可以由平均值和随机误差组成：
[0087]
z(x,y)＝e[z(x,y)]+r(x,y)]＝c+r(x,y)
ꢀꢀꢀ
(16)
[0088]
其中，r(x,y)为点(x,y)处的偏差，其方差为一常数：
[0089]
var[r(x,y)]＝σ2ꢀꢀꢀ
(17)
[0090]
将式(11)代入式(13)有：
[0091][0092]
基于式(14)，式(13)可以改写为：
[0093][0094]
化简得：
[0095][0096]
式(20)就是λi的约束条件之一。令代价函数则有：
[0097][0098]
为了简化描述，令c
ij
＝cov(zi,zj)＝cov(ri,rj)，其中ri＝z
i-c，即点(xi,yi)处的属性值与平均值的差值。则有：
[0099][0100]
令r
ij
＝σ
2-c
ij
，将其与式(20)代入式(22)中，有：
[0101][0102]
为了求得使代价函数j达到最小值时的一组λi，且需要满足式(22)，因此构造拉格朗日函数：
[0103][0104]
其中，φ是拉格朗日乘数，需要求解φ及λi(i＝1,2,...,n)，即：
[0105][0106]
考虑c
ij
＝cov(zi,zj)＝c
ji
，r
ij
＝r
ji
，求解式(25)可得：
[0107][0108]
将式(26)写成线性方程组有：
[0109][0110]
将式(27)写成矩阵形式有：
[0111][0112]
对矩阵求逆即可得到最优解。
[0113]
在对河流两岸像素高程值进行克里金插值后，再对其插值后高程进行多次局部平滑滤波，直到河流高程平滑。
[0114]
s6，在水体高程初步编辑结果的基础上，自动识别连接在一起的异类水体，并在它们的接边位置进行高程的平滑过渡处理。
[0115]
在本实施例中，考虑到前述步骤中采用不同的算法对不同类别的水体进行高程编辑后，异类水体之间的接边处会出现高程跳变的现象，所以根据水体分类结果，识别出连接在一起的异类水体，并找出它们的接边位置；在水体高程初步编辑结果的基础上，在异类水体的接边位置处进行高程的平滑过渡处理，消除高程跳变现象。
[0116]
在异类水体接边时，它们的轮廓像素是相连的，所以可以通过遍历水体的轮廓像素，判断是否有其他类别的水体与其相连，如果有，则记录接边像素的位置。然后在接边位置处，采用反距离加权法来平滑水体高程，实现异类水体接边的自动化平滑过渡处理。
[0117]
综上所述，采用本实施例提供的insar地形产品水体高程编辑处理方法，充分应用了不同类别的水体展现出的不同特性以及实际数据的特点，对存在噪声的insardsm中的水体进行高程编辑处理，在保证水体高程真实的前提下，实现各类水体的平滑处理以及异类水体接边的平滑过渡，整个流程结构清晰，具有流程清晰、实现简单、自动化程度高等优点。
[0118]
以下通过具体案例进一步说明本方法的效果，以澳门特别行政区、广东省珠海市作为实验区(21.97
°
n～22.52
°
n，113.26
°
e～113.63
°
e)，采用条带模式获得的实验区升轨影像一景sar强度影像、一景原始dsm数据进行案例说明。数据详细信息如表1所示。
[0119]
表1实验数据详细信息
[0120][0121]
水体掩膜如图2所示，水体分类结果如图3所示，可以看出水体被分成了海洋、湖泊和河流三类水体。
[0122]
由insar技术得到原始dsm如图4所示，水体编辑后结果如图5所示，可以看出，原始dsm中水体存在了许多噪声，或者高程值是空值(异常值)，经过水体高程编辑处理后，海洋、湖泊和河流各类水体都赋予了与其岸边像素相近的高程值，且效果都十分平滑。
[0123]
本方案旨在保证编辑后水体高程真实性的前提下，消除原始dsm中存在的噪声和异常值(包括空值)，使得编辑后的水体高程呈现出平滑变化的效果，这会导致水体高程值变化幅度减小，因此可以用水体高程编辑前后的标准差变化来定量化评定水体高程编辑效果。
[0124]
表2展示了实验区海洋、湖泊、河流以及整个水体的高程编辑前后的标准差对比。可以看出，在水体编辑后，各类水体以及全部水体的高程标准差在水体编辑后均有所降低，其中海洋和河流的标准差效果提升显著，说明水体高程编辑起到了十分良好的效果。
[0125]
表2水体编辑前后高程标准差对比
[0126][0127]
基于同一发明构思，本实施例还提供了一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现前述的insar地形产品水体高程编辑处理方法。
[0128]
该计算机可读介质包括但不限于任何类型的盘(包括软盘、硬盘、光盘、cd-rom、和磁光盘)、rom、ram、eprom(erasable programmable read-only memory，可擦写可编程只读存储器)、eeprom、闪存、磁性卡片或光线卡片。也就是说，可读介质包括由设备(例如计算机)以能够读的形式存储或传输信息的任何介质。
[0129]
本实施例提供的计算机可读存储介质，与前述的方法具有相同的发明构思及相同
的有益效果，在此不再赘述。
[0130]
以上所述是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明所述原理的前提下，还可以作出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

技术特征：
1.一种insar地形产品水体高程编辑处理方法，其特征在于，包括如下步骤：s1，对水体掩膜进行数据预处理；s2，利用全球高分辨率海岸线矢量数据，对水体掩膜中的海洋进行自动化分类；s3，根据湖泊和河流的形态面积的不同，对水体掩膜中剩余的湖泊和河流进行自动化分类；s4，对海洋高程进行自动化编辑处理，计算海洋像素的大地水准面起伏高度；s5，对湖泊高程进行自动化编辑处理；对河流高程进行自动化编辑处理；s6，基于水体高程初步编辑结果，自动识别连接在一起的异类水体，并在接边位置进行高程的平滑过渡处理。2.根据权利要求1所述的一种insar地形产品水体高程编辑处理方法，其特征在于，s1中将水体提取结果处理成二值掩膜图；对水体掩膜中的空洞进行填补处理；对水体掩膜进行包括平滑、去噪和增强形态学处理。3.根据权利要求2所述的一种insar地形产品水体高程编辑处理方法，其特征在于，s2中根据dsm图幅的地理范围确定矩形矢量，通过矩形矢量和全球海岸线矢量数据进行矢量求交运算，得到图幅地理范围内的海岸线面状矢量；将海岸线面状矢量转换成与dsm图幅分辨率一致的栅格文件，得到与dsm图幅的地理范围及分辨率一致的海洋掩膜栅格文件；根据海洋掩膜文件将水体中的海洋区分出来。4.根据权利要求1所述的insar地形产品水体高程编辑处理方法，其特征在于，s3中将已分类成海洋的水体在水体掩膜中排除，确定剩余水体连通域的范围；再根据湖泊和河流表现出的一般形态的不同，设计指标和阈值，对剩余水体中的湖泊和河流进行自动化分类，并在水体分类结果图中分别进行标记，最终得到水体分类结果。5.根据权利要求4所述的一种insar地形产品水体高程编辑处理方法，其特征在于，根据分类水域面积与水体最大延伸的比值指标t，指标大于阈值a的连通域分类成湖泊，指标小于阈值a的连通域分类成河流。6.根据权利要求1所述的一种insar地形产品水体高程编辑处理方法，其特征在于，s4中利用egm2008地球重力场模型，依据下式计算大地水准面起伏高度：其中，n
gm
为模型大地水准面高，r，θ，λ分别为地心距离、余纬和经度，gm为引力常数与地球质量的乘积，a为参考椭球长半轴，和为完全规格化系数，为完全规格化缔合legendre函数，n
max
为地球重力场模型展开的最高阶数；其中，完全规格化缔合legendre函数采用跨阶次递推法进行计算。7.根据权利要求6所述的一种insar地形产品水体高程编辑处理方法，其特征在于，s5中根据连通域膨胀原理获得湖泊岸边像素的高程值，绘制出分布直方图，选择直方图中的峰值下降的四分之一处的高程值，作为湖泊的统一高程。8.根据权利要求7所述的一种insar地形产品水体高程编辑处理方法，其特征在于，s5
中根据河流两岸像素的高程值确定河流高程，且河流高程由高到低平滑变化；对河流两岸像素的高程值进行去噪处理；根据去噪后的河流两岸像素的高程值进行自动化插值；对河流两岸像素的高程值进行多次局部平滑滤波处理。9.根据权利要求8所述的一种insar地形产品水体高程编辑处理方法，其特征在于，运用克里金插值法对河流两岸像素的高程值进行插值，克里金插值法的公式定义为：其中，是点(x0,y0)的估值，即z0＝(x0,y0)；λ
i
是权重系数，用空间上所有已知点的数据加权求和来估计未知点的值，满足估值与真值之差最小以及无偏估计条件的一套最优系数。10.根据权利要求1所述的一种insar地形产品水体高程编辑处理方法，其特征在于，s6中根据水体分类结果，识别出连接在一起的异类水体，并找出接边位置；基于水体高程初步编辑结果，在异类水体的接边位置处进行高程的平滑过渡处理，消除高程跳变现象。

技术总结
本发明提供了一种InSAR地形产品水体高程编辑处理方法，包括：对水体掩膜进行数据预处理；对水体掩膜中的海洋进行自动化分类；对水体掩膜中剩余的湖泊和河流进行自动化分类；对海洋高程进行自动化编辑处理，计算海洋像素的大地水准面起伏高度；对湖泊高程进行自动化编辑处理；对河流高程进行自动化编辑处理；自动识别连接在一起的异类水体，并在接边位置进行高程的平滑过渡处理。本发明充分应用了不同类别的水体展现出的不同特性以及实际数据的特点，对存在噪声的InSARDSM中的水体进行高程编辑处理，在保证水体高程真实的前提下，实现各类水体的平滑处理以及异类水体接边的平滑过渡具有流程清晰、实现简单、自动化程度高等优点。点。点。


技术研发人员：胡俊 苑艺博 桂容 李志伟
受保护的技术使用者：中南大学
技术研发日：2023.05.24
技术公布日：2023/8/24