CN114527465A - 一种基于永久散射体的地面控制点自动选取方法 - Google Patents

Abstract

本发明涉及一种基于永久散射体的地面控制点自动选取方法，包括：获取多时序SAR影像；从该多时序SAR影像提取满足永久散射体选取策略的像素，作为地面控制点；永久散射体选取策略包括：计算每个像素的平均相干性，选取平均相干性大于预设的低相干阈值的像素点；选取振幅大于最小振幅均值的像素点；选取振幅离差指数小于预设的振幅离差指数阈值的像素点；选取平均相干性小于预设的高相干阈值的像素点；选取坡度小于预设的地形坡度阈值的像素点；根据点的密度、距离和地形条件进一步进行抽稀处理。与现有技术相比，本发明综合考虑了相位稳定性、散射强度、地形坡度和GCP分布等影响因素，在不同区域获得的DEM精度提高了约20％～30％，在平原地区表现最佳。

Description

一种基于永久散射体的地面控制点自动选取方法

技术领域

本发明涉及地面控制点选取技术领域，尤其是涉及一种基于永久散射体的地面控制点自动选取方法。

背景技术

干涉合成孔径雷达(InSAR)是一种从单视复数(SLC)图像中获取相位信息来提取地面高程和监测形变信息的高精度技术。InSAR生成数字高程模型(DEM) 的本质是相位测量和相位到高度的转换。由于SAR的工作范围是微波波段，因此 SAR可以在任何气象条件下使用，具有全天时全天候工作的能力。虽然InSAR具有大面积、空间连续覆盖的优势，但InSAR生成DEM的过程极易受各种因素的影响，例如干扰对的选择、基线和大气环境。在轨道精炼和重去平的过程中引入高质量的地面控制点(GCP)可以校正基线参数并去除平地效应，有效地提高DEM精度，然而，手动选择等常用的GCP选择方法可能会降低DEM的精度。有必要研究新方法来选择合适的GCP进行轨道精炼和重去平，以获得高精度的DEM。

手动选择、现存测量数据、人工角反射器和永久散射体(PS)等多种方法可用于GCP的选择进行轨道精炼和重去平。在沉降监测和DEM生成中人工手动选择 GCP进行轨道精炼和重去平是一种常规的方法。这样的GCP通常是由操作者从 SAR图像中提取的，由于在SAR图像中识别同名点是比较困难的，因此，人工选择GCP的质量与操作者对SAR图像的学习程度、研究区域的熟悉程度以及SAR 的处理经验高度相关。此外，在处理大量SAR图像时手动选择非常耗时。

现存的测量数据也可以用于提取高质量的GCP以用于改进轨道精炼的基线参数并去除干涉相位误差，生成高精度的InSAR DEM。虽然这是一种可靠的方法，但并非每一景SAR图像都有测量数据存在，并且在对每个研究地点进行实地测量时，既费时费力又成本高昂。人工角反射器可以为SAR几何校准、辐射校准和 InSAR处理提供最高质量的GCP，但是人工角反射器需要针对不同的卫星进行专业的设计和制作，并在现场放置，获取它们准确的坐标和姿态。

发明内容

本发明的目的就是为了克服上述现有技术存在在使用干涉合成孔径雷达(InSAR)生成数字高程模型(DEM)的过程中，一般通过手动选择或实地测量来确定地面控制点(GCP)进行轨道精炼和重去平，但这样无法确保为后续的干涉处理提供足够且合适的GCP的缺陷而提供一种基于永久散射体的地面控制点自动选取方法。

本发明的目的可以通过以下技术方案来实现：

一种基于永久散射体的地面控制点自动选取方法，包括以下步骤：

获取多时序SAR影像；

从该多时序SAR影像提取满足永久散射体选取策略的像素，作为地面控制点；

所述永久散射体选取策略包括：

第一选取子策略：计算每个像素的平均相干性，选取平均相干性大于预设的低相干阈值的像素点；

第二选取子策略：选取振幅大于最小振幅均值的像素点；

第三选取子策略：选取振幅离差指数小于预设的振幅离差指数阈值的像素点；

第四选取子策略：选取平均相干性小于预设的高相干阈值的像素点；

第五选取子策略：选取坡度小于预设的地形坡度阈值的像素点。

进一步地，所述低相干阈值的取值在0.37-0.43范围以内。

进一步地，所述振幅离差指数阈值的取值在0.15至0.25范围以内。

进一步地，所述高相干阈值的取值在0.90至0.98范围以内。

进一步地，所述地形坡度阈值为13°至17°范围以内。

进一步地，所述平均相干性的计算表达式为：

式中，γ为相干性，M(i,j)为主图像的复数，S(i,j)为辅图像的复数，所述多时序SAR影像包括主图像和辅图像，i和j分别为方形窗口中的行和列，星号是复数的共轭算子，a和b代表窗口大小，K为多时序SAR影像中SAR图像数量， k为第k张图像，γ_pix为平均相干性，γ_pix-thd为低相干阈值。

进一步地，所述最小振幅均值的计算表达式为：

式中，(m,n)表示覆盖同一区域的SAR图像的总行列数，A_k(e,f)为第k景SAR 图像中e行f列像素的振幅值，A_img-thd为振幅均值最小值，A_img为振幅值。

进一步地，所述振幅离差指数的计算表达式为：

式中，A_pix-disp为振幅离差指数，A_pix-σ是振幅的标准偏差，A_pix是所有图像中相同位置的每个像素的平均振幅，A_pix-disp-thd为振幅离差指数阈值。

进一步地，所述方法还包括：

根据获取的地面控制点，计算各地面控制点间的距离和数量，从而进一步提取均匀分布的地面控制点，然后计算提取的每个地面控制点的坡度，并与相邻的在预设的第一区域内的地面控制点进行坡度比较，保留坡度较低的达到预设的第一数量地面控制点，若该第一区域内地面控制点坡度相同，则随机选择，获取最终的地面控制点。

进一步地，选取的所述地面控制点用于生成数字高程模型的过程中进行轨道精炼和重去平。

与现有技术相比，本发明具有以下优点：

在使用干涉合成孔径雷达(InSAR)生成数字高程模型(DEM)的过程中，一般通过手动选择或实地测量来确定地面控制点(GCP)进行轨道精炼和重去平，但这样无法确保为后续的干涉处理提供足够且合适的GCP。

因此，本发明提出了一种基于永久散射体(PS)的GCP自动选择方法，用于 InSARDEM的生成。auto-PS-GCP方法通过自动定义相干性、振幅和振幅离差指数的阈值来选择PS作为GCP。auto-PS-GCP-Thin方法在auto-PS-GCP的基础上根据点的密度、距离和地形条件进一步进行抽稀处理。

本发明在实施例中使用三段式评价方法来评价结果，包括：1)评价GCP的相位稳定性和强度，2)GCP与参考DEM中同名点高程差的RMSE，以及3)生成的InSAR DEM精度。本发明选取平原、丘陵和山区三个具有代表性的地区验证所提出的方法。结果表明，与手动选择相比，auto-PS-GCP-Thin在不同区域获得的 DEM精度提高了约20％～30％，在平原地区表现最佳，DEM精度为4.7051m。此外，与auto-PS-GCP相比，auto-PS-GCP-Thin得到的DEM精度提高了约20％-40％。可以得出结论，所提出的方法在具有不同地形条件的区域是有效且可靠的。

附图说明

图1为本发明实施例中提供的一种星载InSAR的几何关系示意图；

图2为本发明实施例中提供的一种基于永久散射体的地面控制点自动选取方法，并用于InSAR DEM生成的轨道精炼和重去平的流程图；

图3为本发明实施例中提供的一种GCP抽稀流程图；

图4为本发明实施例中提供的三个研究区域示意图，图4中，(a)中国的DEM 地图，(b)平原地区的北京市中心，(c)丘陵地区的平凉市，以及(d)山地地区的张家口市；(b-d)的底图是Landsat-8OLI图像的合成图；

图5为本发明实施例中提供的三个研究区域的空间基线(a-c)和时间基线(d-f)示意图；平原地区主图像获取日期为2020年1月29日，丘陵地区为2020年10 月20日，山地地区为2020年9月1日；

图6为本发明实施例中提供的干涉图(a-c)和相干图(d-f)；平原地区的干涉对是2020年1月29日的主图像和2020年2月22日的辅图像，丘陵地区的干涉对是2020年10月20日的主图像和2020年11月13日的辅图像，以及山地地区的干涉对是2020年9月1日的主图像和2020年9月25日的辅图像；

图7为本发明实施例中提供的手动、auto-PS-GCP和auto-PS-GCP-Thin方法的GCP分布情况图；选择的GCP覆盖在相应的相干系数图上，图7中对于三个研究区域，(a-c)表示手动选择，(d-f)表示auto-PS-GCP选择，以及(g-i)表示 auto-PS-GCP-Thin选择；

图8为本发明实施例中提供的手动(a-c)、auto-PS-GCP(d-f)和auto-PS-GCP-Thin(g-i)得到的GCP的相位直方图；

图9为本发明实施例中提供的选择得到的三种GCP生成的DEM，图9中，(a-c) 为手动选点生成的DEM，(d-f)为auto-PS-GCP选点生成的DEM，以及(g-i)为 auto-PS-GCP-Thin选点生成的DEM；

图10为本发明实施例中提供的随机点叠加在外部SRTM DEM的示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。

因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。

在本发明的描述中，需要说明的是，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括一个或者更多个该特征。在本申请的描述中，“多个”的含义是两个或两个以上，除非另有明确具体的限定。

实施例1

GCP也可以从PS中提取。PS是指具有强散射特性、高相干性和稳定相位的各种地面目标。在SAR图像中，典型的PS通常包括桥梁、建筑和裸露岩石的顶角。此类PS已应用于GB-SAR以生成高精度DEM。由于PS的时空稳定性，它们可用于SBAS-InSAR进行沉降监测，PS技术还与干扰处理相结合，从TanDEM-X 提取DEM。由于PS具有强散射特性且能在长时间序列中保持相位稳定，因此该方法可以产生足够的GCP进行轨道精炼和重去平，以减少干涉相位误差，从而提高DEM精度。通常基于相干性、振幅均值和振幅离差指数识别和提取PS。一些学者定义了提取PS的相干性阈值，而其他学者则使用合适的振幅离差指数作为阈值。事实上，例如相干性或振幅均值之类的单一准则只考虑了PS的一个方面。例如，相位色散指数阈值仅考虑了相位稳定性，因此提取的PS可能对失相关和系统噪声比较敏感，相干性考虑了高信噪比(SNR)，振幅离差指数主要考虑相位稳定性，这样可能会得到一些在轨道精炼和重去平方面表现不佳的PS。

PS的选择需要考虑相位稳定性和散射强度，因此应预先定义适当的阈值以提取合适的点并丢弃不可靠的点。并且干涉图的生成与地形影响有关，有效的GCP 通常分布在平坦的区域，可以将合适的地形坡度作为度量标准以帮助选择高质量的 GCP。所选GCP的分布和密度也会显着影响InSAR中的轨道精炼和重去平，通过点的剔除来校准GCP可能会提高应用的效果。因此，GCP的选择应充分考虑地形条件和GCP分布的影响。此外，在InSAR DEM生成中有几个与GCP选择相关的阶段，每个阶段都应进行充分评价以确保其验证。本研究旨在解决GCP选择的关键问题：1)本实施例能否将地形条件与SAR图像相结合，以选择合适的GCP进行轨道精炼和重去平；2)如何校准所选GCP的分布；3)如何在InSAR DEM 生成的每个阶段评价GCP。

为了解决上述问题，本实施例提供了一种新方法：基于永久散射体的地面控制点自动选取方法(auto-PS-GCP)，考虑了相位、散射特性、相干性、地形坡度和点的分布来提取更加可靠且有用的GCP，应用于InSAR流程中生成高精度DEM。该方法应用于平原地区(北京市中心)、丘陵地区(平凉)和山区(张家口)三个试验区的InSAR DEM生成。然后使用三段式评价方法来评价GCP的相位和强度、 GCP本身的精度和生成DEM的精度。

本实施例提供的一种基于永久散射体的地面控制点自动选取方法，包括以下步骤：

获取多时序SAR影像；

从该多时序SAR影像提取满足永久散射体选取策略的像素，作为地面控制点；

永久散射体选取策略包括：

第一选取子策略：计算每个像素的平均相干性，选取平均相干性大于预设的低相干阈值的像素点；

第二选取子策略：选取振幅大于最小振幅均值的像素点；

第三选取子策略：选取振幅离差指数小于预设的振幅离差指数阈值的像素点；

第四选取子策略：选取平均相干性小于预设的高相干阈值的像素点；

第五选取子策略：选取坡度小于预设的地形坡度阈值的像素点。

上述第一选取子策略、第二选取子策略、第三选取子策略、第四选取子策略和第五选取子策略的执行顺序不做限定，各后续子策略执行时均对前一子策略的选取结果进行进一步像素点提取。

优选的，第一选取子策略、第二选取子策略、第三选取子策略、第四选取子策略和第五选取子策略依次执行。

作为一种优选的实施方式，低相干阈值的取值在0.37-0.43范围以内，以去除低相干性点。

作为一种优选的实施方式，振幅离差指数阈值的取值在0.15至0.25范围以内，以提取稳定的GCP。

作为一种优选的实施方式，高相干阈值的取值在0.90至0.98范围以内，以进一步选择GCP。

作为一种优选的实施方式，地形坡度阈值为13°至17°范围以内，以进一步识别平坦地区的GCP。将坡度结合到GCP选择中去除受叠掩和阴影效应影响的点。

作为一种优选的实施方式，方法还包括：

根据获取的地面控制点，计算各地面控制点间的距离和数量，从而进一步提取均匀分布的地面控制点，然后计算提取的每个地面控制点的坡度，并与相邻的在预设的第一区域内的地面控制点进行坡度比较，保留坡度较低的达到预设的第一数量地面控制点，若该第一区域内地面控制点坡度相同，则随机选择，获取最终的地面控制点。

作为一种优选的实施方式，选取的地面控制点用于生成数字高程模型的过程中进行轨道精炼和重去平。

具体的，本实施例的实施过程包括以下步骤：

1、方法执行过程

1.1、InSAR DEM生成方法

利用InSAR技术提取DEM的基本原理是利用具有干涉成像能力的两副SAR 天线(或一副天线重复观测)来获取同一地区具有一定视角差的两幅SLC图像，并根据其干涉相位信息反演地表的高程信息，从而生成相应地区的DEM。主辅SAR 影像间存在相位差，可结合观测卫星的轨道参数生成高精度、高分辨率的地面高程信息。如图1所示，由于两次观测时卫星与目标的距离不同，同一目标(P)记录的相位(

和

)不同。干涉相位

可由下式给出：

其中φ1和φ2是来自目标的回波信号的相位，λ是雷达波长，R1和R2是两个天线到观测目标P的斜距，ΔR是斜距差。

干扰相位(φint)由地形相位(φtop)、地形形变相位(φdef)、平地相位(φflat)、大气相位(φatm)和噪声相位(φnoi)组成：

去除平地效应、大气和噪声相位后，可获得目标的真实地形相位

基于图1所示的几何关系，目标(P)的高程(h)可由下式给出：

其中B是SAR天线A₁和A₂之间的空间基线，H和θ是A₁的高度和视角，α是基线倾角。

在InSAR DEM生成过程中，稳定、高质量的GCP可以有效地提高轨道精炼和重去平的性能，从而提高生成的DEM精度。本实施例提出了用于GCP选取的 auto-PS-GCP方法，工作流程如图2所示：a)InSAR DEM生成；b)使用 auto-PS-GCP和auto-PS-GCP-thin方法选择GCP；c)三段式评价。最后将auto-PS-GCP和auto-PS-GCP-Thin方法应用于三个试验区(分别位于平原、丘陵和山区)的多时相SAR图像以验证提出的方法。本实施例根据三阶段评价方法评价了所选GCP的相位稳定性和散射强度、它们的高程精度以及生成的InSAR DEM 精度。

1.2、auto-PS-GCP选择方法

auto-PS-GCP选择方法综合考虑时间稳定性和空间连续性来选择合适数量的 PS。如图2所示，通过自动确定相干性、振幅均值和振幅离差指数阈值以及抽稀来选择GCP。为确保所有图像的振幅能够真实反映地面目标的散射情况，本实施例在选择GCP之前进行了辐射校准。

1.2.1、相干系数阈值

相干性(γ)是提取PS的第一个准则。γ是通过方形窗口(例如3×3)中相邻像素的复数来计算的。计算方法可由下式给出：

其中M(i,j)为主图像的复数，S(i,j)为辅图像的复数，这两幅图像构成了干涉像对，i和j分别为方形窗口中的行和列，星号是复数的共轭算子，a和b代表窗口大小(这里a＝b＝3)。像素的SNR可由其相干性来估计：

较高的相干系数表示较高的SNR，能够得到较好的干涉图质量。考虑到时间序列中的所有K景图像中相同像素的相干性(γ)，每个像素都有一个平均相干性 (γ_pix)：

其中K表示SAR图像数量，k为第k张图像。选择大于预定义的相干系数阈值(γ_pix-thd)的点作为PS。

低的阈值只能去除严重失相干的点，而高的阈值可能会错过有用的PS，因此，本实施例使用低相干阈值(γ_pix-thd1)去除水体和植被覆盖区域等区域(图2中的 Step1)，然后使用高相干阈值(γ_pix-thd2)进一步选择PS(图2中的Step4)。

1.2.2、振幅均值阈值

每幅SAR图像的振幅均值可以表示其散射强度，对于长时间序列的SAR图像，在这些振幅均值中都有一个最小值(A_img-thd)。本实施例定义最小振幅均值(图2 中的Step 2)作为提取PS的阈值，保留高于阈值的像素作为PS：

其中(m,n)表示覆盖同一区域的SAR图像的总行列数，A_k(e,f)为第k景SAR 图像中e行f列像素的振幅值。振幅均值最小值(A_img-thd)可以去除Step1中振幅较低但相干性虚高的点。

1.2.3、振幅离差指数阈值

振幅离差指数是选择PS的第三个标准(图2中的Step 3)。该指数可用于反映相位稳定性。振幅离差指数(A_pix-disp)可以由下式给出：

其中A_pix-σ是振幅的标准偏差，A_pix是所有图像中相同位置的每个像素的平均振幅。本实施例定义一个合适的阈值(A_pix-disp-thd)来提取振幅离差指数小于阈值的点作为PS。

这种方法不受相邻像素的影响，因为它是由同一位置不同图像的像素计算的，并且本实施例考虑了每个像素的相位稳定性，从PS候选点中删除了具有严重失相干的点。公式(8)中的这种方法适用于具有高SNR的像素，这些像素可根据公式 (4)提取。

1.2.4、地形坡度阈值

本实施例将地形坡度作为PS选择的标准之一(图2中的Step 5)，以减少陡峭地形对干涉精度的影响。在Step1-4生成的具有良好相位稳定性和强散射强度的 PS的基础上进一步提取位于平坦区域中的点。因此，定义合适的阈值(SLOPE_thd) 来提取坡度小于SLOPE_thd的PS。坡度可以通过以下方式计算：

其中p和q是高程梯度，分别表示东西方向和南北方向的高程变化率，180/π是从弧度到角度的转换因子。

1.2.5、GCP抽稀

基于相干性、振幅均值、振幅离差指数和坡度选择得到初始GCP(即PS)。当应用于InSAR流程生成DEM时，这些GCP用于校正基线参数并去除平地效应。此外，GCP的数量和分布也会影响基线参数和平地效应的去除。

对于使用上述方法得出的高密度GCP的区域，本实施例考虑点的密度，并将坡度作为优先级，通过剔除其中一些点来对点进行优化(图3)。首先考虑点距离和点数量，从候选点中选择合适的GCP，确保点尽可能均匀分布，然后计算每个 GCP的坡度并与其相邻像素进行比较，保留坡度较低的点。当这些GCP具有相似的坡度时，本实施例从候选点中随机选择。

应用平均最近邻(ANN)度量来分析抽稀前后GCP的分布，该指标根据每个 GCP到其最近邻GCP的平均距离计算得到。观察到的平均距离(NNObserved)由每个GCP与其最近邻GCP之间的距离计算得出。z-score是统计显著性的度量，表明GCP是否随机分布。在99％的置信水平下，如果z-score大于2.58，则GCP是分散的，如果z-score小于-2.58，则GCP是聚集的。

1.3、三段式评价

本实施例提出了一种三段式评价方法来评价选择GCP和生成的DEM(图2c)。第一阶段评价选择的GCP的相位和强度。如果所选GCP的干涉相位聚集为0rad 附近，并且GCP的强度比所在图像的强度均值高5倍以上，则认为所选GCP是有效的。第二阶段比较所选GCP和参考DEM中同名点的高程，本实施例认为在平坦区域的RMSE低于10m和在丘陵、山地区域低于20m的GCP在InSAR处理中是有效的。第三阶段将外部高精度DEM与由auto-PS-GCP和auto-PS-GCP-Thin方法生成的InSAR DEM进行比较。

3、研究区和数据集

3.1、研究区和数据集

本实施例选择平原、丘陵和山地的三个研究区域(图4)来测试auto-PS-GCP 和auto-PS-GCP-Thin方法。平原地区为位于华北平原的北京市中心(图4b)，丘陵地区为位于黄土丘陵的甘肃平凉市(图4c)，山地地区为位于燕山的河北张家口市 (图4d)。北京市区约500平方公里，平凉、张家口地区约600平方公里。三个地区的植被稀疏，减少了植被覆盖对干涉测量的影响。

本实施例使用了来自双极化C波段的Sentinel-1A SLC数据。表1为三个试验区的24景IW模式下的升轨数据，分辨率为5m×20m，极化方式为VV (scihub.copernicus.eu)。轨道参数精度对InSAR处理的结果有很大影响，轨道误差会以残余条纹的形式存在于干涉图中，因此本实施例使用POD Precise Orbit Ephemerides数据(slqc.asf.alaska.edu)来校正轨道参数。

本实施例使用日本宇宙航空研究开发机构(search.asf.alaska.edu)的高级陆地观测卫星1号项目(ALOS)的PALSAR传感器获得的具有高精细空间分辨率DEM 作为参考DEM。该参考DEM在Fine Beam Single偏振模式下的地形分辨率为12.5 m，垂直分辨率为12m。本实施例选择SRTM DEM作为外部DEM进行精度分析，SRTM在为期11天的任务中配备了两个天线，并免费提供近乎全球规模的高分辨率DEM。

表1 Sentinel-1A SAR数据集

2.2、数据预处理

InSAR DEM生成的工作流程如图2a所示。本实施例同时考虑时间和空间基线，从每个研究区域的Sentinel-1A数据(18个月内)的24景SAR图像中选择了主图像(图5)，其他23张SAR图像是辅图像，然后为每个研究区域生成23景干涉图。对于平原地区和山地地区，平均空间基线分别为48.5435m和47.5130m，两个研究区的时间基线比较合适。但是对于丘陵地区，空间基线平均为48.2428m，是比较适合生成干涉图的，然而，时间基线最大至479天，这可能会降低干涉图的相干性。

将主从图像配准后，进行干涉处理生成干涉图，并采用Goldstein滤波法对干涉图进行干涉滤波。然后通过最小费用流方法对滤波后的干涉图进行相位解缠，以得到与地形信息相对应的真实相位。然后，导入选择的GCP进行轨道精炼和重去平来改进干涉相位。由于GCP的选择是在地理坐标系中进行的，因此本实施例将 GCP转换为斜距坐标系以生成InSAR DEM。最后，通过相位高程转换和地理编码获得三个研究区的DEM。

图6为23景干涉对中获得的较为典型的干涉图和相干图，两者共同反映了起伏的地形信息。较高的相关系数和较清晰的干涉条纹表明干涉结果的可信度较高。对于平原(图6a和d)和山区(图6c和f)，干涉条纹清晰，大部分区域具有较高的相干性。丘陵区域整体的相干性低于其他两个区域，但大部分区域仍显示出清晰的干涉条纹(图6b)和高相干性(图6e)。因此，可以使用这些干涉图和相干图有效地生成InSAR DEM。

4、试验与结果

4.1、基于auto-PS-GCP方法选择的GCP

在对解缠的干涉相位进行地理编码前，将选择得到的GCP导入InSAR流程中进行轨道精炼和重去平以改善相位。根据提出的auto-PS-GCP方法，本实施例综合 GCP的相干系数、振幅均值、振幅离差指数、地形坡度和分布来为每个研究区域选择GCP，如表2所示。相干系数(γ_pix-thd1)的低阈值被定义为0.4，振幅均值阈值(A_img-thd)为振幅均值最小值A_img，三个区域的振幅离差指数阈值(A_pix-disp-thd) 分别为0.25、0.15和0.20，然后定义相干系数的高阈值(γ_pix-thd2)以从上述得到的 PS候选点中选择更有用的GCP。其中，平原研究区是一个有大量强散射物体的城区，因此北京市中心的相干系数的高阈值(γ_pix-thd2)设置为0.98，以得到合适的 GCP数量。最后，选择坡度低于15°的GCP以确保这些点位于平坦区域。使用 auto-PS-GCP导出的GCP显示出相对较强的聚集，使用抽稀方法可以弱化聚集情况。

表2 auto-PS-GCP方法定义的阈值和GCPs的分布情况评价

图7为手动、auto-PS-GCP和auto-PS-GCP-Thin方法的GCP分布情况。本实施例使用手动和auto-PS-GCP-Thin方法的选择得到相同数量的GCP，便于比较它们的分布情况和后续生成的DEM精度。北京市中心有许多不变的散射体，导致平原地区有大量的GCP(图7d)。在对auto-PS-GCP的GCP进行抽稀处理后(图7d-f)，auto-PS-GCP-Thin的点显示出更合理的分布(图7g-i)；具体来说，平原、丘陵和山区的GCP数量减少了70％左右。表2显示手动和auto-PS-GCP-Thin的z-score 大于2.58，表明它们的GCP是分散分布的；而auto-PS-GCP的z-score小于-2.58，表明它们的GCP是聚集分布的。抽稀处理将自动选择的GCP从聚集转化为分散，从而改善GCP分布。

3.2、GCP的第一阶段评价

本实施例对所有SAR图像对进行了干涉处理，进行强度和相位的评价(第一阶段评价)由于页数限制，本实施例仅针对每个研究区域呈现了一个图像对的结果 (表3)。由于每个图像对之间的时间间隔很短，地形形变的影响，天气条件等外部因素的影响很小。表3比较了整幅SAR图像和三种方法选择的GCP之间的平均强度。由auto-PS-GCP和auto-PS-GCP-Thin选择的GCP的平均强度高于整个SAR 图像，手动选择的GCP的平均强度与整个图像的平均强度相差不大。以平原地区为例，表3显示手动选择的GCPs的平均强度是整幅图像的两倍，而auto-PS-GCP 和auto-PS-GCP-Thin的GCPs强度要大得多，是整个图像的20倍以上。丘陵和山区的平均强度也有类似的情况，表明这两种方法选择的GCPs具有很强的散射性。

表3 GCPs平均强度的评价

图8显示了GCP在干涉结果中的相位直方图。手动选择的GCP的相位在三个研究区域中没有显示出一致的分布，而auto-PS-GCP和auto-PS-GCP-Thin选择的GCP相位都围绕在0rad附近。这表明auto-PS-GCP和auto-PS-GCP-Thin选择的GCPs具有良好的相位稳定性，证明了两种方法的有效性和可行性。对于 auto-PS-GCP-Thin方法，选定的GCP不仅具有高散射和相位稳定性，而且在平坦区域中表现为近似分散的分布，表明这些GCP可用在随后的InSAR DEM生成中。

3.3、GCP的第二阶段评价

将三种方法选择的GCP导入InSAR流程中进行轨道精炼和重去平生成DEM。图9显示了使用InSAR和所选择的GCP生成的DEM。本实施例将生成的DEM 中GCP的高程与参考DEM中的同名点进行比较，计算它们的RMSE(表4)。 auto-PS-GCP-Thin选择得到的GCP的RMSE最小，表现出最好的性能；然而， auto-PS-GCP选择得到的GCP的RMSE大于手动，表明GCP的聚集可能会导致 DEM精度降低。尽管手动和auto-PS-GCP-Thin选择得到的GCP数量相同，但在使用具有强散射特性和分散分布的GCP生成的InSAR DEM效果更好。这表明提出的auto-PS-GCP-Thin方法具有良好的性能。

表4选择得到的GCP与参考DEM中的同名点的RMSE

3.4、生成的DEM的第三阶段评价

为了评价DEM的精度，本实施例使用分辨率为30m的SRTM DEM作为外部 DEM，并在每个研究区域生成100个随机点(图10)。从这两种类型的DEM中提取这些随机点的高程，并进行比较以分析精度(表5)。与auto-PS-GCP-Thin生成的GCP相关的DEM在所有三个研究区域中都具有最高的精度。由于GCP聚集的影响，与auto-PS-GCP生成的GCP相关的DEM在所有三个研究区域中的准确度最低。抽稀前后的精度比较表明，抽稀后的GCP产生了精度更高的DEM。三种地貌类型中，平原地区精度最好，丘陵地区精度较差；然而，正如第二阶段评价所证实的那样，auto-PS-GCP-Thin在所有三个研究领域都表现最佳。

表5生成的DEMs和外部DEMs之间的高程差值精度分析

4、讨论

具有稳定相位和强散射特性的像素对于InSAR中的轨道精炼和重去平至关重要。本实施例提出了一种基于PS的自动选择方法(auto-PS-GCP-Thin)来提取GCP 以进行轨道精炼和重去平。考虑相干性、振幅、振幅离差指数和地形坡度进行GCP 选择。特别是，考虑到点密度和坡度优先级，将聚集的GCP进行抽稀处理。选择平原、丘陵和山区三个地形区域作为研究区域，并将这些选择得到的GCP导入 InSAR处理以生成DEM。使用三段式评价来验证提出的方法和得到的结果。第一阶段表明auto-PS-GCP-Thin得到的GCPs具有良好的相位稳定性和较强的强度，第二阶段表明auto-PS-GCP-Thin得到的GCPs对于不同地形区域的RMSE约为6～12 m，第三阶段表明，auto-PS-GCP-Thin生成的DEM的RMSE在不同区域约为4～15m。

阈值的定义对于GCP的选择至关重要，本实施例考虑了相干性、振幅、振幅离差指数、地形坡度和GCP分布等五个方面。相干性的低阈值设置为0.4，以去除低相干性点，振幅均值阈值设置为振幅均值最小值来去除表现出虚假高相干性的点。振幅离差指数的阈值通常设置为低于0.25，以确保得到能够在长时间序列中保持较好稳定性的点，在本研究中，本实施例将振幅离差指数阈值定义为0.15～0.25以提取稳定的GCP。然后本实施例将高相干性阈值定义为0.90～0.98，以进一步选择GCP，这些选定的阈值适用于不同领域。此外，GCP的选择还受到SAR成像叠掩和阴影效应的影响，这与地形坡度有很大的关系，因此，本实施例参考早期的地形分类研究，将坡度阈值定义为15°，以进一步识别平坦地区的GCP。将坡度结合到GCP 选择中去除受叠掩和阴影效应影响的点。

综上，GCP的选择由相位稳定性、散射强度和GCP分布的共同决定。 auto-PS-GCP选择的GCP具有稳定的相位和高散射强度，而手动选择的GCP显示分散(即更合适)的分布。根据第二阶段评价，例如在平原地区，auto-PS-GCP选择得到的GCP的RMSE为9.2138m，手动选择得到的GCP的RMSE为8.5438m， auto-PS-GCP-Thin选择得到的GCP的RMSE仅为6.5562m。结果表明，GCP分布是更重要的影响后续InSAR DEM生成的因素。auto-PS-GCP-Thin的GCP平均强度低于auto-PS-GCP(参见表3)，这可能是由于删除了一些聚集在一个非常小的区域的高散射强度的点，这两种方法得到的GCP的平均散射强度仍是大于手动选择的GCP的平均散射强度。

众所周知，地形条件对InSAR DEM精度有很大影响，这一点已被本研究证实。结果表明，这三种方法在平原地区的精度最高，其中auto-PS-GCP-Thin表现最好，高程精度为4.7051m。然而，这种方法在丘陵地区得到的DEM的RMSE为14.9749 m，在所有三个地区中报告的精度最低(参见表5)。此外，丘陵地区的相干性在所有三个地区中是最弱的(参见图6)，相应地，选择得到的GCP也是最少的。与其他方法相比，auto-PS-GCP-Thin方法可以产生的DEM具有更高的精度。例如，使用Sentinel-1A数据生成布什尔市的DEM的案例得到的RMSE为18.53m。

为了更好地比较这三种方法，本实施例采用了三段式评价方法来分析GCP的相位和强度、RMSE以及DEM精度。对于手动选择，虽然它的GCP是可用的，但人为选择是可变的，因此每次的RMSE都不同。auto-PS-GCP-Thin考虑了相位稳定性、散射强度、地形坡度和GCP分布，因此结果是一致的，可以应用于多种情况。在本研究中，本实施例将auto-PS-GCP-Thin应用于平原、丘陵和山地三个不同地形的区域，这些区域包括了密集的建筑和植被覆盖的区域，验证了所提出的 GCP选择和三段式评价方法。

5、结论

用于轨道精炼和重去平的高质量GCP对于InSAR DEM生成至关重要。GCP 可以从PS候选点中提取，这些候选点通常是在长时间序列图像中保持强相位稳定性和散射强度的地面目标。本研究提出了两种GCP选择方法(auto-PS-GCP和 auto-PS-GCP-Thin)，在平原、丘陵和山地三个区域进行了测试。首先利用低相干系数和振幅均值获得PS候选点，然后采用振幅离差指数和高相干系数从上述PS 候选点中提取更高质量的PS，并进一步考虑地形坡度的影响，去除高坡度点。最后，提出了一种点抽稀方法，通过剔除一些聚集在很小区域内的点来校正GCP的分布。结果表明，与手动相比，auto-PS-GCP-Thin得到的DEM精度在不同区域提高了约20％～30％，其中在平原区域DEM精度为4.7051m，效果最好。此外，与 auto-PS-GCP相比，auto-PS-GCP-Thin将DEM精度提高了约20％-40％。

本实施例的GCP选择方法综合考虑了相位稳定性、散射强度、地形坡度和GCP 分布等影响因素，三个案例研究表明，auto-PS-GCP-Thin是一种可靠的可应用与 InSAR DEM生成的GCP选取方法。因此，所提出的方法在具有不同地形条件的各个区域都是有效的。

以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术人员无需创造性劳动就可以根据本发明的构思做出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。

Claims (10)

1.一种基于永久散射体的地面控制点自动选取方法，其特征在于，包括以下步骤：

获取多时序SAR影像；

从该多时序SAR影像提取满足永久散射体选取策略的像素，作为地面控制点；

所述永久散射体选取策略包括：

第一选取子策略：计算每个像素的平均相干性，选取平均相干性大于预设的低相干阈值的像素点；

第二选取子策略：选取振幅大于最小振幅均值的像素点；

第三选取子策略：选取振幅离差指数小于预设的振幅离差指数阈值的像素点；

第四选取子策略：选取平均相干性小于预设的高相干阈值的像素点；

第五选取子策略：选取坡度小于预设的地形坡度阈值的像素点。

2.根据权利要求1所述的一种基于永久散射体的地面控制点自动选取方法，其特征在于，所述低相干阈值的取值在0.37-0.43范围以内。

3.根据权利要求1所述的一种基于永久散射体的地面控制点自动选取方法，其特征在于，所述振幅离差指数阈值的取值在0.15至0.25范围以内。

4.根据权利要求1所述的一种基于永久散射体的地面控制点自动选取方法，其特征在于，所述高相干阈值的取值在0.90至0.98范围以内。

5.根据权利要求1所述的一种基于永久散射体的地面控制点自动选取方法，其特征在于，所述地形坡度阈值为13°至17°范围以内。

6.根据权利要求1所述的一种基于永久散射体的地面控制点自动选取方法，其特征在于，所述平均相干性的计算表达式为：

式中，γ为相干性，M(i,j)为主图像的复数，S(i,j)为辅图像的复数，所述多时序SAR影像包括主图像和辅图像，i和j分别为方形窗口中的行和列，星号是复数的共轭算子，a和b代表窗口大小，K为多时序SAR影像中SAR图像数量，k为第k张图像，γ_pix为平均相干性，γ_pix-thd为低相干阈值。

7.根据权利要求1所述的一种基于永久散射体的地面控制点自动选取方法，其特征在于，所述最小振幅均值的计算表达式为：

式中，(m,n)表示覆盖同一区域的SAR图像的总行列数，A_k(e,f)为第k景SAR图像中e行f列像素的振幅值，A_img-thd为振幅均值最小值，A_img为振幅值。

8.根据权利要求1所述的一种基于永久散射体的地面控制点自动选取方法，其特征在于，所述振幅离差指数的计算表达式为：

式中，A_pix-disp为振幅离差指数，A_pix-σ是振幅的标准偏差，A_pix是所有图像中相同位置的每个像素的平均振幅，A_pix-disp-thd为振幅离差指数阈值。

9.根据权利要求1所述的一种基于永久散射体的地面控制点自动选取方法，其特征在于，所述方法还包括：

根据获取的地面控制点，计算各地面控制点间的距离和数量，从而进一步提取均匀分布的地面控制点，然后计算提取的每个地面控制点的坡度，并与相邻的在预设的第一区域内的地面控制点进行坡度比较，保留坡度较低的达到预设的第一数量地面控制点，若该第一区域内地面控制点坡度相同，则随机选择，获取最终的地面控制点。

10.根据权利要求1所述的一种基于永久散射体的地面控制点自动选取方法，其特征在于，选取的所述地面控制点用于生成数字高程模型的过程中进行轨道精炼和重去平。

Images