CN113514829A - 面向InSAR的初始DSM的区域网平差方法 - Google Patents
面向InSAR的初始DSM的区域网平差方法 Download PDFInfo
- Publication number
- CN113514829A CN113514829A CN202110783009.4A CN202110783009A CN113514829A CN 113514829 A CN113514829 A CN 113514829A CN 202110783009 A CN202110783009 A CN 202110783009A CN 113514829 A CN113514829 A CN 113514829A
- Authority
- CN
- China
- Prior art keywords
- initial
- dsm
- initial dsm
- coordinates
- east
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 68
- 239000011159 matrix material Substances 0.000 claims abstract description 35
- 238000003384 imaging method Methods 0.000 claims description 42
- 230000000977 initiatory effect Effects 0.000 claims description 12
- 238000013441 quality evaluation Methods 0.000 claims description 7
- 238000004458 analytical method Methods 0.000 claims description 6
- 230000009466 transformation Effects 0.000 claims description 6
- 239000000126 substance Substances 0.000 claims description 3
- 238000006243 chemical reaction Methods 0.000 abstract description 13
- 238000013507 mapping Methods 0.000 abstract description 11
- 239000013598 vector Substances 0.000 description 20
- 238000004364 calculation method Methods 0.000 description 16
- 238000005516 engineering process Methods 0.000 description 9
- 238000010587 phase diagram Methods 0.000 description 7
- 238000012545 processing Methods 0.000 description 7
- 238000004590 computer program Methods 0.000 description 6
- 238000012937 correction Methods 0.000 description 6
- 230000000694 effects Effects 0.000 description 4
- 238000003860 storage Methods 0.000 description 4
- 230000001174 ascending effect Effects 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 238000001914 filtration Methods 0.000 description 3
- 238000012952 Resampling Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 238000005305 interferometry Methods 0.000 description 2
- 238000012805 post-processing Methods 0.000 description 2
- 238000001303 quality assessment method Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000013459 approach Methods 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000008030 elimination Effects 0.000 description 1
- 238000003379 elimination reaction Methods 0.000 description 1
- 238000005111 flow chemistry technique Methods 0.000 description 1
- 238000009432 framing Methods 0.000 description 1
- 230000002452 interceptive effect Effects 0.000 description 1
- 238000012886 linear function Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 230000009897 systematic effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/88—Radar or analogous systems specially adapted for specific applications
- G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
- G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
- G01S13/9094—Theoretical aspects
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- Physics & Mathematics (AREA)
- Electromagnetism (AREA)
- Computer Networks & Wireless Communication (AREA)
- General Physics & Mathematics (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种面向InSAR的初始DSM的区域网平差方法，该方法包括以下步骤：步骤S1，生成初始DSM及附加数据；步骤S2，计算控制点在主影像中的像素坐标；步骤S3，确定各个初始DSM的邻接关系矩阵；步骤S4，均匀抽取各个初始DSM重叠处的连接点，并计算所述连接点在所述初始DSM中的像素坐标以及与所述初始DSM相邻的影像中的像素坐标；步骤S5，确定所述连接点的经相位高程转换之后的高程值；步骤S6，构建最小二乘平差方程，对所述初始DSM进行区域网平差。本发明能够支持业务化InSAR地形测绘过程中，流程化生产的初始DSM产品的精度的进一步提升。
Description
技术领域
本发明涉及图像处理技术领域，特别涉及一种面向InSAR的初始DSM的区域网平差方法。
背景技术
合成孔径雷达干涉(InSAR)技术已经成为了最为有效的全球测图手段之一。2000年，SRTM提供了全球56°S～60°N范围内的数字表面模型(DSM)数据，其标称精度达到了16m，格网大小30米，在全球尺度的科研以及工程应用中发挥了不可估量的作用。在2016年，德宇航宣布使用TanDEM-X双星绕飞星座完成了全球DSM的测绘，90％地区的精度为3.5米，如果不考虑极低精度区域，全球90％地区的精度则达到了0.8米，数据已经可以用于精细的军事化定位和军事目标识别。InSAR地形测绘的巨大潜力，使得InSAR技术在工程应用上迎来了一个新的时代。然而我国SAR卫星技术发展缓慢，现阶段尚没有民用的业务化InSAR地形测绘卫星，这导致了我国的应用技术发展迟缓。InSAR地形测绘技术，特别是面向大规模全球尺度的地形测绘技术，是我国将来发展的方向，也是现在技术的短板。
InSAR地形测绘是一个高度自动化的过程。在自动化处理过程中，为了精度的控制和提升，需要使用原始数据进行InSAR的流程化处理，完成数据的标准分幅、配准、重采样、干涉、去平、基线精化、滤波、解缠、相高转换、地理编码，生成具有地理坐标的DSM数据。这种处理方法是现阶段常用的处理方法，其中存在两个较为严重的工程化应用难题。第一，基线精化过程中，需要使用外部的DSM数据或其它高精度三维坐标，由于国外地面控制的缺乏，使得这种方法的全球可用性极低。第二，这种方法生成的DSM中包含了流程化处理过程中的各类误差，从其严密成像几何模型可以发现，这些误差主要包括了几何定位相关的误差源，即卫星成像时刻的位置、速度、斜距、多普勒等参数，以及高程解算相关的误差源，即干涉基线和相位等参数，这些误差会在1m到10m之间波动，造成DSM的产品精度稳定性较差。
本发明将前述自动化处理得到的DSM称为初始DSM，针对初始DSM误差较大，精度不统一，存在明显的接边误差的情况，引入地面高程控制点进行平差处理，最终解算得到区域网平差的各平差系数，为初始DSM的精度改善提供了修正参数。
发明内容
针对上述问题，本发明提供了一种面向InSAR的初始DSM的区域网平差方法，可充分地消除由于现有技术的限制和缺陷导致的一个或多个问题。
本发明提供了一种面向InSAR的初始DSM的区域网平差方法，所述方包括以下步骤：
步骤S1，利用合成孔径雷达干涉(InSAR)生成初始DSM及附加数据，所述附加数据包括主影像头文件、从影像头文件、控制点质量评估文件、初始DSM对应的SRTM和初始DSM对应的SRTM误差；
步骤S2，构建斜距-多普勒方程和椭球方程，并利用构建的所述斜距-多普勒方程和椭球方程，计算所述控制点质量评估文件中的控制点在主影像中的像素坐标；
步骤S3，根据各个所述初始DSM的经纬度范围，确定各个所述初始DSM的邻接关系矩阵；
步骤S4，基于所述邻接关系矩阵，均匀抽取各个所述初始DSM重叠处的连接点，并计算所述连接点在所述初始DSM中的像素坐标以及与所述初始DSM相邻的影像中的像素坐标；
步骤S5，基于所述连接点在所述初始DSM中的像素坐标，确定所述连接点的经相位高程转换之后的高程值；
步骤S6，基于所述控制点的像素坐标和所述连接点在所述初始DSM中的像素坐标以及所述连接点的高程，构建最小二乘平差方程，对所述初始DSM进行区域网平差，确定所述初始DSM的误差多项式系数。
优选的，步骤S3具体包括以下子步骤：
步骤S3.1，构建各个所述初始DSM的初始邻接关系矩阵，所述初始邻接关系矩阵F表达式为：
其中，f(i,j)为第i列第j行的初始DSM，g为东西向的初始DSM数量，k为南北向的初始DSM数量；
步骤S3.2，对步骤S3.1中构建的初始邻接关系矩阵中的所有初始DSM逐一进行迭代分析，确定各个初始DSM的邻接关系矩阵。
优选的，步骤S3.2具体包括以下子步骤：
步骤S3.2.1，确定位于最西北角的初始DSM，并将该DSM记为f(0,0)；
步骤S3.2.2，分别确定与所述最西北角的初始DSM的东向邻接的初始DSM和南向邻接的初始DSM；
步骤S3.2.3，依次判断除最西北角初始DSM的其他初始DSM的邻接初始DSM，以确定所有初始DSM的邻接初始DSM，并最终确定初始DSM的邻接关系矩阵。
优选的，将所述最西北角的初始DSMf(0,0)的经度范围记为[lons1,lone1]，纬度范围记为[lats1,late1]；将任一其他初始DSMf(i,j)的经度范围记为[lons2,lone2]，纬度范围记为[lats2,late2]；
其中，确定与所述最西北角的初始DSM的东向邻接的初始DSM具体包括以下子步骤：
(1)判断任一其他初始DSMf(i,j)与所述最西北角的初始DSMf(0,0)在经度方向是否有重叠，具体包括：
如果lons2>lons1，且lons2<lone1，则判定f(i,j)与f(0,0)在经度方向有重叠，并继续执行步骤(2)，否则判定无重叠，即，判定f(i,j)与f(0,0)在东西向无邻接关系，不需要执行下述步骤；
(2)计算f(i,j)与f(0,0)的纬度重叠范围；
纬度重叠范围表达为
lat_com＝min(late2,late1)-max(lats2,lats1)
其中，lat_com为纬度方向的重叠范围，min()为求取最小值，max()为求取最大值；
(3)判断f(i,j)与f(0,0)是否具备东向邻接关系；
具体的，如果lat_com>(late1-lats1)*0.7或者lat_com>(late2-lats2)*0.7，则判定f(i,j)与f(0,0)具备东向邻接关系，并赋值i＝0+1，j＝0，以寻找f(1,0)的东向和南向邻接DSM；否则，判定f(i,j)与f(0,0)不具备东向邻接关系，并确定与所述最西北角的初始DSM的南向邻接的初始DSM。
优选的，确定与所述最西北角的初始DSM的南向邻接的初始DSM具体包括以下子步骤：
(1)判断f(i,j)与f(0,0)在纬度方向是否有重叠，具体包括：
如果lats2>lats1，且lats2<late1，则判定f(i,j)与f(0,0)在纬度方向有重叠，并继续执行步骤(2)，否则判定无重叠，即，判定f(i,j)与f(0,0)在南北向无邻接关系，不需要执行下述步骤；
(2)计算f(i,j)与f(0,0)的经度重叠范围。
经度重叠范围表达为
lon_com＝min(lone2,lone1)-max(lons2,lons1)
其中，lon_com为纬度方向的重叠范围，min()为求取最小值，max()为求取最大值；
(3)判断f(i,j)与f(0,0)是否具备南向邻接关系；
其中，如果lon_com>(lone1-lons1)*0.7或者lon_com>(lone2-lons2)*0.7，则判定f(i,j)与f(0,0)具备南向邻接关系，并赋值i＝0，j＝0+1；否则，判定f(i,j)与f(0,0)不具备南向邻接关系。
优选的，步骤S4具体包括以下子步骤：
步骤S4.1，根据卫星成像时刻的飞行方向确定所述初始DSM在地理编码过程中的镜像变换关系；
步骤S4.2，根据通过步骤S3确定的邻接关系矩阵以及通过步骤S4.1确定的镜像变换关系，确定所述连接点在所述初始DSM中的像素坐标，其中，所述连接点包括南向连接点和东向连接点；
步骤S4.3，将步骤S4.2获取的南向连接点和东向连接点的像素坐标转换为地理坐标；
步骤S4.4，计算步骤S4.3解算的地理坐标在邻接关系矩阵中与所述初始DSM的相邻影像中的像素坐标。
优选的，在步骤S5中，最小二乘平差方程的表达式如下：
BX＝L
其中，
X＝(a1 b1 c1…aK bK cK)′
其中，x为控制点的东西方向像素坐标，
为第k景初始DSM中的第Mk个控制点的东西方向像素坐标，
为与第一景初始DSM有邻接关系的初始DSM中第Mk+Pk个连接点的东西方向坐标，
为第L个初始DSM中第Mi+Pi个连接点的东西向坐标，
为与第L个初始DSM相邻的初始DSM中第Mi+Pi个连接点的东西向坐标，y为控制点的南北方向像素坐标，
为第k景初始DSM中的第Mk个控制点的南北方向像素坐标，
为与第一景初始DSM有邻接关系的初始DSM中第Mk+Pk个连接点的南北方向坐标，
为第L个初始DSM中第Mi+Pi个连接点的南北向坐标，
为与第L个初始DSM相邻的初始DSM中第Mi+Pi个连接点的南北向坐标，
为第L个初始DSM中第M+p个连接点的高程，
为与第L个初始DSM相邻的初始DSM中第M+p个连接点的高程，Mk为第k景初始DSM中的控制点数量，K为参与平差的初始DSM影像总数量，P为第P个连接点，Pk为第k景初始DSM影像中的连接点总数量，Ad为与处于最西北角的第一景初始DSM相邻的初始DSM，L为第L景具有连接点的初始DSM，R为第L景初始DSM的邻接初始DSM；(aK bK cK)为第K景初始DSM的待求多项式参数，
为第K景初始DSM的第M个控制点的高程误差。
优选的，上述最小二乘平差方程的最小二乘解为
X＝(BTB)-1BTL
本发明还提供了一种计算机程序，其中，当所述计算机程序由计算设备的处理器执行时，执行本发明实施例所述的方法。
本发明还提供了一种计算机可读存储介质，其中，所述计算机可读存储介质上存储计算机程序，当所述计算机程序由计算设备的处理器执行时，执行本发明实施例所述的方法。
本发明给出了初始DSM误差的消除方法，可消除初始DSM的系统误差，提高DSM的精度，使其满足最终的DSM精度需求。
本发明针对InSAR流程化处理得到的初始DSM中存在几何定位及高程解算的各类误差，从而导致初始DSM误差较大、精度不统一、存在明显的接边误差等问题，提出了初始DSM的区域网平差方法。本发明使用初始DSM及附加数据，通过控制点以及连接点的几何坐标转换，在SAR影像空间中构建区域网平差方程，并解算获取初始DSM的误差多项式系数，为初始DSM的精度改进提供了修正参数。
附图说明
图1为根据本发明实施例的InSAR初始DSM区域网平差方法的流程图。
图2为根据本发明实施例的初始DSM的邻接关系示意图。
图3为根据本发明实施例的侧视方向、飞行方向与镜像变换关系示意图。
具体实施方式
下面对本发明进行更为详细的描述，描述过程中参照了上述附图，并使用示例性实施例进行了详解。
针对背景技术中提出的问题，本发明提出了一种面向InSAR初始DSM的区域网平差方法。这种方法针对多景影像进行区域网平差算法研究，完成系统误差的建模和消除，提高最终的DSM产品精度。
本发明目的在于确保DSM产品的稳定性，因此针对传统处理方法的应用难题，一方面去除了基线精化的步骤，消除了对地面三维坐标的依赖性。另一方面，本发明专注于DSM的后处理，即，使用ICESat数据进行DSM的区域网平差，消除前述各类误差源带来的高程误差。本发明为InSAR全球测绘提供了技术基础。
总体来讲，本发明主要包括以下步骤：
步骤1，准备区域网平差过程中需要的各类数据，所有数据都是InSAR流程化生产DSM之后生成的，具体的，利用合成孔径雷达干涉(InSAR)生成初始DSM及附加数据，所述附加数据包括主影像头文件、从影像头文件、控制点质量评估文件、初始DSM对应的SRTM和初始DSM对应的SRTM误差。本发明中，这些数据均为输入数据。
步骤S2，由于区域网平差的过程是在SAR影像空间进行的，因此所有的坐标均需转换为像素坐标，本步骤将控制点的地理坐标转换为在主影像中的像素坐标；具体的，构建斜距-多普勒方程和椭球方程，并利用构建的所述斜距-多普勒方程和椭球方程，计算所述控制点质量评估文件中的控制点在主影像中的像素坐标。
步骤S3，进行DSM邻接性判断。基于所提供的DSM的地理坐标，对输入的DSM进行东南西北的重新定位。具体的，根据各个所述初始DSM的经纬度范围，确定各个所述初始DSM的邻接关系矩阵。
步骤S4，基于所述邻接关系矩阵，均匀抽取各个所述初始DSM重叠处的连接点，并计算所述连接点在所述初始DSM中的像素坐标以及与所述初始DSM相邻的影像中的像素坐标。其中，基于DSM的邻接关系，确定连接区域，并从连接区域抽取连接点。对于东西侧邻接的数据来说，抽取西侧DSM的东边界数据作为连接点，对于南北侧邻接的数据来说，抽取北侧DSM的南边界数据作为连接点。
步骤S5，基于所述连接点在所述初始DSM中的像素坐标，确定所述连接点的经相位高程转换之后的高程值。由于连接点的坐标已知，因此需要联合步骤1中的输入数据，进行高程信息的抽取。
步骤S6，基于所述控制点的像素坐标和所述连接点在所述初始DSM中的像素坐标以及所述连接点的高程，构建最小二乘平差方程，对所述初始DSM进行区域网平差，确定所述初始DSM的误差多项式系数。其中，区域网平差的观测方程为距离向坐标x和方位向坐标y的一次函数，每景影像仅包含三个多项式参数，理论上来说，三个控制点即可解算对应参数，为了结果的准确性，实际过程中将使用所有的控制点和连接点数据进行参数解算。
下面结合具体实施例，对本发明的方案做详细说明。
如图1所示，本发明提出的面向InSAR的初始DSM的区域网平差方法主要包括以下步骤：
步骤S1，利用合成孔径雷达干涉(InSAR)生成初始DSM及附加数据，所述附加数据包括主影像头文件、从影像头文件、控制点质量评估文件、初始DSM对应的SRTM和初始DSM对应的SRTM误差。
本发明是在利用InSAR技术生成的初始DSM的基础上进行的，所以首先需要利用InSAR技术生成初始DSM数据。步骤S1主要包括以下子步骤：
步骤S1.1，采集同一地区的两幅SAR影像，将其中的一副SAR影像作为主影像，另一副作为从影像，对所述主影像和从影像进行配准。具体的，找到从影像与主影像的同名点，并拟合多项式，建立从影像到主影像的空间对应关系，对从影像进行重采样，使从影像和主影像的像素一一对应，生成重采样之后的从影像及对应的头文件。
步骤S1.2，对主影像和从影像进行干涉，得到干涉相位图。具体的，将主影像与步骤S1.1重采样获取的从影像进行逐像素共轭相乘，得到干涉相位图。
步骤S1.3，计算平地相位并从所述干涉相位图中去除所述平地相位，得到去平相位图。具体的，基于本地平均高程，计算步骤S1.2对应的影像空间内的平地相位，并从步骤S1.2获取的干涉相位图中，减掉平地相位，得到去平相位图。
步骤S1.4，使用干涉图滤波方法对步骤S1.3获取的去平相位图进行滤波，得到滤波相位图。
步骤S1.5，将步骤S1.4获取的滤波相位图进行相位解缠，得到解缠相位图。
步骤S1.6，计算干涉相位图在影像中心点的绝对相位。具体的，使用步骤S1.1获取的同名点信息解算绝对相位。
步骤S1.7，将所述解缠相位图与所述绝对相位相加，并使用主影像从文件以及步骤S1.1得到的重采样之后的从影像头文件，将相位转换到高程，获取相位高程转换之后的高程图。
步骤S1.8，使用所述相位高程转换之后的高程图，将高程信息投影到地理坐标空间，得到初始DSM。
以上即通过InSAR技术生成初始DSM及附加数据的过程，其中，生成DSM附加数据包括初始DSM对应的高程以及高程误差文件，这些文件以及文件的说明如下表1所示。
表1.InSAR生产的初始DSM以及附加数据的说明。
根据本发明的优选实施例，选取某地区的TanDEM-X卫星星座数据作为数据源，采用InSAR技术进行数据的流程化处理之后，生成的初始DSM共有4景，具体包括西北向的初始DSM、东北向的初始DSM、西南向的初始DSM、东南向的初始DSM，其基本参数如表2所示。初始DSM均以“景”为单位进行存储，每景初始DSM均由一对TanDEM-X数据生成，这一对数据由两颗绕飞的SAR卫星提供，其中第一颗卫星在2006年发射，记为TSX，第二颗卫星在2010年发射，记为TDX。这一对TanDEM-X影像的头文件均可提供斜距-多普勒(Range-Doppler,R-D)模型的计算参数，其作用类似。R-D模型是连接初始DSM坐标与像素坐标的重要桥梁。根据本发明的优选本实施例将TDX卫星对应的影像作为主影像，将TSX卫星对应的影像作为从影像，每景初始DSM均由一对干涉影像对(主影像和从影像)生成。4景初始DSM对应的4个TDX的头文件信息表格如表3–表6所示。表格中给出了影像成像的所有参数，后续计算过程中的几何参数均将从表格中获取。
表2.本发明优选实施例用到的4景初始DSM及基本参数
表3.NW对应的TDX数据头文件信息表格
表4.NE对应的TDX数据头文件信息表格
表5.SW对应的TDX数据头文件信息表格
表6.SE对应的TDX数据头文件信息表格
步骤S2，构建斜距-多普勒方程和椭球方程，并利用构建的所述斜距-多普勒方程和椭球方程，计算所述控制点质量评估文件中的控制点在主影像中的像素坐标。
区域网平差的过程是在SAR影像空间进行的，因此所有的坐标均需转换为像素坐标。本步骤提供的控制点像素坐标，后续的步骤S4提供的连接点像素坐标，两者共同参与步骤S6的区域网平差。
具体的，步骤S2包括以下子步骤：
步骤S2.1，构建斜距-多普勒方程(R-D方程)和椭球方程；
具体的，采用步骤S1中表3-表6给出的TDX数据头文件信息，结合SAR成像过程中常用的R-D方程确定主影像四个角点的坐标。R-D方程是进行SAR影像定位的重要方程，它通过表3-表6中提供的斜距、多普勒参数，结合地面点的椭球参数，建立某一点P的像素坐标(x,y)与地面坐标(XP,YP,ZP)之间的对应关系。R-D方程共包含三个子方程，即斜距方程、多普勒方程以及椭球方程。
步骤S2.1具体包括以下步骤：
步骤S2.1.1，构建斜距方程。
斜距方程表达的是，成像时刻的卫星坐标与地面坐标之间的距离等于头文件中获取的斜距。斜距方程如下所示：
其中，R0是近地点斜距，ΔR是斜距分辨率，(XS,YS,ZS)是成像时刻t的卫星坐标，(XP,YP,ZP)是地面点的三维坐标，x为地面点在所述主影像中的距离向像素坐标。
各参数的解算具体如下。
步骤S2.1.1.1，获取斜距参数。
近地点斜距R0和斜距分辨率ΔR是已知量，分别是表3–表6中的主影像的近地点斜距和距离分辨率。
步骤S2.1.1.2，获取成像时刻参数t。
y为地面点在所述主影像中的方位向像素坐标，y与成像时刻t之间存在下述关系
t＝t0+Δt·(y-1) (2)
其中，t0为方位向成像起始时间，即步骤S1中，表3–表6给出的方位向起始时间，Δt为方位向成像时间间隔，即步骤S1中，表3–表6给出的方位向单行时间。
步骤S2.1.1.3，计算卫星三维坐标。
卫星的三维坐标(XS,YS,ZS)是采用卫星的N个状态矢量插值而成。首先使用N个状态矢量构建卫星位置与卫星成像时刻的多项式关系，优选的，采用4阶多项式，以成像时刻的卫星坐标X轴的状态矢量为例，其与卫星成像时刻之间的多项式关系可表达为：
其中，(t1 t2…tN)为N个主影像的状态矢量对应的时刻，(a1 a2 a3 a4 a5)为状态矢量的X轴分量对成像时刻的多项式拟合系数，
为卫星成像时刻的N个状态矢量的X轴分量。在优选的实施例中，表3–表6均包含了19个主影像的状态矢量，因此N＝19。那么在方程(3)中，状态矢量的X轴分量、N个状态矢量对应的时间均为已知量，仅有(a1a2 a3 a4 a5)为未知量。令
那么依据最小二乘原理，即可得到
不失一般的，将表3的19个状态矢量以及对应的成像时刻带入公式(2)–(6)中，得到
L＝(-2309389.6640-2309678.0530…-2314049.2310)， (8)
同样的，将成像时刻的卫星坐标Y轴的状态矢量带入公式(2)–(6)，可以得到状态矢量的Y轴分量对成像时刻的多项式拟合系数
以及状态矢量的Z轴分量对成像时刻的多项式拟合系数
步骤S2.1.1.4，构建斜距方程。
结合步骤S2.1.1.1–步骤S2.1.1.3的参数，构建的斜距方程为：
将表3中的参数带入上述方程中，得到斜距方程的实施例如下所示。
至此，斜距方程中未知量只包含两类，即地面坐标(XP,YP,ZP)以及像素坐标(x,y)。
步骤S2.1.2，构建多普勒方程。
多普勒方程描述的是卫星与地面的相对运动带来的多普勒效应。成像时刻的卫星速度(VXS VYS VZS)与卫星侧视方向矢量((XS-XP)(YS-YP)(ZS-ZP))之间的叉乘代表了对应的多普勒效应。
构建的多普勒方程的具体表达式为：
其中，(VXS VYS VZS)为成像时刻的卫星速度，(XS,YS,ZS)为成像时刻的卫星坐标，(XP,YP,ZP)是地面点的三维坐标，fd为多普勒效应值，λ为雷达波长，R0是近地点斜距，ΔR是斜距分辨率，x为地面点在所述主影像中的距离向像素坐标。
在星载SAR处理过程中，会通过星上偏航导引以及地面的成像，保证多普勒效应为0，称为零多普勒效应。因此在本发明中，方程右侧的值恒为0。在上述方程中，未知量仅包含成像时刻的卫星速度。
成像时刻的卫星速度解算过程如下：
步骤S2.1.2.1，计算卫星的速度矢量。
卫星速度矢量是采用卫星的N个速度矢量插值而成。首先使用N个速度矢量构建卫星速度与卫星成像时间的多项式关系，优选的，采用4阶多项式。其构建方式与步骤S2.1.1.3中计算卫星的三维坐标的过程完全一致，此处不再赘述。以NW的TDX影像为例，采用与公式(7)–(13)中成像时刻的卫星位置相同的解算方法，可得到成像时刻的卫星速度与成像时间t之间的关系：
式中，(VXS VYS VZS)为成像时刻的卫星速度，(d1 d2 d3 d4 d5)为VXs对成像时刻的多项式拟合系数，(e1 e2 e3 e4 e5)为VYs对成像时刻的多项式拟合系数，(f1 f2 f3 f4 f5)为VXs对成像时刻的多项式拟合系数，t0为方位向成像起始时间，即步骤S1中，表3–表6给出的方位向起始时间，Δt为方位向成像时间间隔，即步骤S1中，表3–表6给出的方位向单行时间。y是地面点在所述主影像中的方位向像素坐标，t是成像时刻。
与步骤S2.1.1.3类似，将表3中的参数代入公式(15)，即可得到卫星速度矢量的实施例如下。
步骤S2.1.2.2，构建多普勒方程。
式中，(VXS VYS VZS)为成像时刻的卫星速度，(XS,YS,ZS)为成像时刻的卫星坐标，(XP,YP,ZP)是地面点的三维坐标。
至此，多普勒方程中未知量只包含两类，即地面坐标(XP,YP,ZP)，以及地面点在所述主影像中的方位向和距离向像素坐标(x,y)。
步骤S2.1.3，构建椭球方程。
构建的椭球方程为：
其中，(XP,YP,ZP)是地面点的三维坐标，h为地面点的高程，Re与Rp分别为椭球的赤道半径与极半径，在WGS84坐标系统中，赤道半径为6378137.00m，极半径为6356752.31m。因此严格来讲，如果想要获取地面点的精确坐标，需确保地面点高程为已知量。在本实施例中，精确的高程信息可通过ICESat数据获取，此时椭球方程中的未知参数只有地面坐标。
步骤S2.2，利用步骤2.1.1-2.1.3构建的斜距方程、多普勒方程和椭球方程，解算控制点在主影像中的像素坐标。
步骤S1生成的初始DSM包含各种误差，这些误差在影像方位向和距离向表现为低阶多项式，因此平差的核心方程是在影像的像素坐标空间内建立的。在构建方程之前，需将控制点的地理坐标转换为像素坐标。
坐标解算过程中需要用到步骤S2.1.1、S2.1.2以及S2.1.3中的三个方程。在方程中，未知量共包含三类，即像素坐标、地面点坐标、地面点高程。在此步骤中，需要借助控制点地面经纬高等三维坐标，以解算控制点对应的像素坐标。
以NW对应的初始DSM为例，其提供的控制点质量评估文件由ICESat激光点构成，文件中包含了2161个点，每个点包含经度、纬度、高程以及高程误差。不失一般的，选用第一个ICESat点，其经纬度坐标为109.378°E，35.0299°N，高程为791.216m，使用R-D模型可解算得到对应的列数为1285.046，行数为3230.392(列数、行数即控制点的像素坐标)。表2–表5给出了本实施例中，4个初始DSM对应的ICESat坐标由经纬度转到像素坐标的转换结果。
表7.NW对应的DSM的ICESat坐标转换结果。
| 索引 | 经度(°) | 纬度(°) | 列数 | 行数 | 高程(m) | 高程误差(m) |
| 1 | 109.378 | 35.029 | 1285.046 | 3230.392 | 791.216 | 1.875 |
| 2 | 109.378 | 35.031 | 1286.824 | 3220.940 | 812.873 | -6.005 |
| …… | …… | …… | …… | …… | …… | …… |
| 2160 | 109.319 | 35.478 | 2125.049 | 475.223 | 951.999 | -8.027 |
| 2161 | 109.319 | 35.481 | 2132.862 | 456.368 | 934.159 | -11.670 |
表8.NE对应的DSM的ICESat坐标转换结果。
| 索引 | 经度(°) | 纬度(°) | 列数 | 行数 | 高程(m) | 高程误差(m) |
| 1 | 109.190 | 34.648 | 1942.185 | 2954.766 | 383.666 | -3.446 |
| 2 | 109.186 | 34.673 | 1992.822 | 2802.262 | 389.483 | -2.944 |
| …… | …… | …… | …… | …… | …… | …… |
| 2949 | 109.388 | 35.077 | 1269.884 | 95.080 | 1046.253 | -3.539 |
| 2950 | 109.386 | 35.085 | 1286.741 | 47.582 | 1038.191 | 1.302 |
表9.SW对应的DSM的ICESat坐标转换结果。
表10.SE对应的DSM的ICESat坐标转换结果。
| 索引 | 经度(°) | 纬度(°) | 列数 | 行数 | 高程(m) | 高程误差(m) |
| 1 | 109.456 | 34.574 | 1995.792 | 3539.833 | 350.129 | -1.385 |
| 2 | 109.456 | 34.575 | 1998.975 | 3529.530 | 350.420 | -2.906 |
| …… | …… | …… | …… | …… | …… | …… |
| 410 | 109.459 | 34.660 | 2079.761 | 2949.161 | 353.111 | -2.593 |
| 411 | 109.459 | 34.662 | 2082.940 | 2938.859 | 354.363 | -1.038 |
步骤S3，根据各个所述初始DSM的经纬度范围，确定各个所述初始DSM的邻接关系矩阵。
需要说明的是，步骤S2是计算得到控制点的像素坐标，步骤S3-S5是计算连接点的像素坐标，步骤S6利用控制点的像素坐标和连接点的像素坐标来平差。所以，步骤S2和步骤S3-S5可并行进行。
步骤S3依据表2给出的初始DSM的经纬度范围(初始DSM的经纬度范围为初始DSM的基本参数，属已知数据)，确定各个初始DSM的邻接关系矩阵，其邻接矩阵按照由西向东，由北向南的方式建立。
步骤S3具体包括以下子步骤：
步骤S3.1，构建各个所述初始DSM的初始邻接关系矩阵。
其中，各个所述初始DSM的初始邻接关系矩阵F表达式为：
其中，f(i,j)为第i列第j行的初始DSM，g为东西向的初始DSM数量，k为南北向的初始DSM数量。
上式中，f(i,j)为第i列第j行的初始DSM，其南北东西四个方向的初始DSM分别为f(i,j+1)、f(i,j-1)、f(i+1,j)、f(i-1,j)。
步骤S3.2，对步骤S3.1中构建的邻接关系矩阵中的所有初始DSM逐一进行迭代分析，确定各个初始DSM的邻接关系矩阵。
步骤S3.2具体包括以下子步骤：
步骤S3.2.1，确定位于最西北角的初始DSM，并将该DSM记为f(0,0)。
根据初始DSM的经纬度范围，判断经度最小、纬度最大的DSM对应的编号。在本实施例中，初始DSM的经纬度范围在表2中，其中经度最小、纬度最大的初始DSM对应表2中的NW编号。将该DSM记为f(0,0)。
步骤S3.2.2，分别确定与所述最西北角的初始DSM的东向邻接的初始DSM和南向邻接的初始DSM。
具体的，将最西北角的初始DSMf(0,0)的经度范围记为[lons1,lone1]，纬度范围记为[lats1,late1]。记任一其他初始DSMf(i,j)的经度范围为[lons2,lone2]，纬度范围为[lats2,late2]。
“确定与所述最西北角的初始DSM的东向邻接的初始DSM”具体包括以下子步骤：
(1)判断任一其他初始DSMf(i,j)与所述最西北角的初始DSMf(0,0)在经度方向是否有重叠，具体包括：
如果lons2>lons1，且lons2<lone1，则判定f(i,j)与f(0,0)在经度方向有重叠，并继续执行步骤(2)，否则判定无重叠，即，判定f(i,j)与f(0,0)在东西向无邻接关系，不需要执行下述步骤。
(2)计算f(i,j)与f(0,0)的纬度重叠范围。
纬度重叠范围表达为
lat_com＝min(late2,late1)-max(lats2,lats1) (20)
其中，lat_com为纬度方向的重叠范围，min()为求取最小值，max()为求取最大值。
(3)判断f(i,j)与f(0,0)是否具备东向邻接关系。
一般来说，要求具备东向邻接关系的影像，其南北纬度覆盖范围应大致相同。在本实施例中，要求南北纬度覆盖的70％以上是一致的，即：
如果lat_com>(late1-lats1)*0.7或者lat_com>(late2-lats2)*0.7，则判定f(i,j)与f(0,0)具备东向邻接关系，并赋值i＝0+1，j＝0，以寻找f(1,0)的东向和南向邻接DSM。否则，判定f(i,j)与f(0,0)不具备东向邻接关系，并执行步骤“确定与所述最西北角的初始DSM的南向邻接的初始DSM”。
“确定与所述最西北角的初始DSM的南向邻接的初始DSM”具体包括以下子步骤：
(1)判断f(i,j)与f(0,0)在纬度方向是否有重叠，具体包括：
如果lats2>lats1，且lats2<late1，则判定f(i,j)与f(0,0)在纬度方向有重叠，并继续执行步骤(2)，否则判定无重叠，即，判定f(i,j)与f(0,0)在南北向无邻接关系，不需要执行下述步骤。
(2)计算f(i,j)与f(0,0)的经度重叠范围。
经度重叠范围表达为
lon_com＝min(lone2,lone1)-max(lons2,lons1) (21)
其中，lon_com为纬度方向的重叠范围，min()为求取最小值，max()为求取最大值。
(3)，判断f(i,j)与f(0,0)是否具备南向邻接关系。
类似的，一般来说，要求具备南向邻接关系的影像，其东西经度覆盖范围应大致相同。在本实施例中，要求东西纬度覆盖的70％以上是一致的，即：
如果lon_com>(lone1-lons1)*0.7或者lon_com>(lone2-lons2)*0.7，则判定f(i,j)与f(0,0)具备南向邻接关系，并赋值i＝0，j＝0+1。否则，判定f(i,j)与f(0,0)不具备南向邻接关系。
步骤S3.2.3，依次判断除最西北角初始DSM的其他初始DSM的邻接初始DSM，以确定所有初始DSM的邻接初始DSM，并最终确定初始DSM的邻接关系矩阵。
其他初始DSM的邻接初始DSM的判断方式与步骤S3.2.3的判断方式类似，在此不再赘述。
在本实施例中，依据上述步骤计算西北、东北、西南、东南4景初始DSM的邻接关系矩阵。西北角的初始DSM编号为NW，表达为f(0,0)；东北角的初始DSM编号为NE，与f(0,0)在东向有邻接关系，表达为f(1,0)；西南角的初始DSM编号为SW，与f(0,0)在南向有邻接关系，表达为f(0,1)。第二次迭代时，针对f(1,0)进行判断，其初始DSM编号为NE，东向无初始DSM；东南角的初始DSM编号为SE，与f(1,0)在南向有邻接关系，表达为f(1,1)。第三次迭代时，针对f(0,1)进行判断，其初始DSM编号为SW，南向无初始DSM；东南角的初始DSM编号为SE，与f(0,1)在南向有邻接关系，表达为f(1,1)。第四次迭代时，针对f(1,1)进行判断，其初始DSM编号为SE，东向和南向均无初始DSM。至此，遍历完毕，针对西北、东北、西南、东南4景初始DSM的邻接关系矩阵，得到的最终的邻接关系矩阵为：
其中，F为邻接矩阵，f(0,0)为第0列第0行的初始DSM，f(1,0)为第1列第0行的初始DSM，f(0,1)为第0列第1行的初始DSM，f(1,1)为第1列第1行的初始DSM。NW、NE、SW、SE分别为西北、东北、西南、东南4景初始DSM。根据本发明实施例的初始DSM的邻接关系如图2所示。
步骤S4，基于所述邻接关系矩阵，均匀抽取各个所述初始DSM重叠处的连接点，并计算所述连接点在所述初始DSM中的像素坐标以及与所述初始DSM相邻的影像中的像素坐标。
步骤S4计算过程中需要使用步骤S2中的斜距-多普勒(Range-Doppler,R-D)方程。
步骤S4具体包括以下子步骤：
步骤S4.1，根据卫星成像时刻的飞行方向确定所述初始DSM在地理编码过程中的镜像变换关系。
侧视方向、飞行方向与镜像变换关系如图3所示。SAR影像为侧视成像，成像时的侧视方向与飞行方向的夹角称为方向角，在表3-表6中有相关参数。方向角为90°时，为右侧视，方向角为-90°时，为左侧视。此外，航向角接近0°或者360°时，为升轨飞行，航向角接近180°时，为降轨飞行。按照表3–表6提供的SAR影像的成像参数，结合图3的镜像变换关系，可得到4景初始DSM的镜像情况如下所示。
表11.影像地理编码过程中的镜像变换关系
步骤S4.2，根据通过步骤S3确定的邻接关系矩阵以及通过步骤S4.1确定的镜像变换关系，确定所述连接点在所述初始DSM中的像素坐标，其中，所述连接点包括南向连接点和东向连接点。
处理过程中，从最西北方向开始，逐一进行连接点选取。假设要处理的初始DSM为f(i,j)，即第i列第j行的初始DSM。对每个初始DSM进行南向和东向连接点分析。步骤S4.2具体包括以下子步骤：
步骤S4.2.1，根据步骤S3提供的邻接关系，以及S4.1提供的镜像关系，对初始DSMf(i,j)进行南向连接点分析。步骤S4.2.1具体包括以下子步骤：
步骤S4.2.1.1，如果要处理的初始DSMf(i,j)为右侧视，降轨，则应从TDX影像坐标空间中选择靠近影像坐标空间最后一行的点为连接点，连接点行坐标表达式为：
yi,j＝l-b (23)
其中，l为初始DSM对应的TDX影像行数，b为方位向无效像素行数，在本实施例中，b优选为10，yi,j为第i列，第j行的初始DSM中连接点的行坐标。
第i列，第j行的初始DSM中，连接点列坐标表达式为
xi,j＝c×Δm+bs×s(c＝0,1,2,...,nx-1) (24)
其中，c为从0到nx-1的整数，nx为行方向的连接点数，在本实施例中，nx优选100，bs距离向无效像素列数比例，s为TDX的影像列数，Δm为连接点列坐标间隔，表达为：
其中，bs为距离向无效像素列数比例，用来避免可能出现的无效黑色边框，在本实施例中，bs优选为0.1。s为TDX的影像列数。
步骤S4.2.1.2，如果f(i,j)对应的影像为右侧视，升轨影像，则应从原始数据中选择靠近第一行的点为连接点，第i列，第j行的初始DSM中连接点行坐标表达式为
yi,j＝b (26)
其中，b为方位向无效像素行数，在本实施例中，优选为10。
连接点的列坐标表达式与公式(23)相同。
步骤S4.2.1.3，如果f(i,j)对应的影像为左侧视，降轨影像，则应从原始数据中选择靠近最后一行的点为连接点，连接点行坐标的表达式与公式(22)相同。列坐标的表达式与公式(23)相同。
步骤S4.2.1.4，如果f(i,j)对应的影像为左侧视，升轨影像，则应从原始数据中选择靠近第一行的点为连接点，连接点行坐标的表达式与公式(25)相同。列坐标的表达式与公式(23)相同。
在本实施例中，依据上述步骤，代入表3–表6中各影像的参数，可得到选取的连接点坐标如下表所示。
表12.南向连接点的坐标
| 初始DSM编号 | 列坐标 | 行坐标 |
| NW | 214,231……，1915 | 3223 |
| NE | 208,225……1860 | 3533 |
| SW | 无 | 无 |
| SE | 无 | 无 |
步骤S4.2.2，根据步骤S3提供的邻接关系，以及S4.1提供的镜像关系，进行东向的连接点分析。步骤S4.2.2具体包括以下子步骤：
步骤S4.2.2.1，如果f(i,j)对应的影像为右侧视，降轨影像，则应从原始数据中选择靠近第一列的点为连接点，连接点列坐标表达式为
xi,j＝b (27)
其中，b为可能出现的距离向无效像素列数，在本实施例中，优选为10。
连接点行坐标表达式为
yi,j＝c×Δm+bl×l(c＝0,1,2,...,ny) (28)
其中，c为从0到nx-1的整数，ny为行方向的点数，在本实施例中，ny优选100，Δm为连接点行坐标间隔，表达为
bl为可能出现的距离向无效像素列数比例，用来避免可能出现的无效黑色边框，在本实施例中，bl优选为0.1，l为初始DSM对应的TDX的影像行数。
步骤S4.2.2.2，如果f(i,j)对应的影像为右侧视，升轨影像，则应从原始数据中选择靠近最后一列的点为连接点，连接点列坐标表达式为
xi,j＝s-b (30)
其中，b为距离向无效像素列数，在本实施例中，优选为10。s为初始DSM对应的TDX影像列数。
连接点的行坐标表达式与公式(28)相同。
步骤S4.2.2.3，如果f(i,j)对应的影像为左侧视，降轨影像，则应从原始数据中选择靠近最后一列的点为连接点，连接点列坐标的表达式与公式(30)相同，行坐标的表达式与公式(28)相同。
步骤S4.2.2.4，如果f(i,j)对应的影像为左侧视，升轨影像，则应从原始数据中选择靠近第一列的点为连接点，连接点列坐标的表达式与公式(27)相同。行坐标的表达式与公式(28)相同。
在本实施例中，依据步骤S4.2.2.1–4.2.2.3，代入表3–表6中各影像的航向角、方位角、行列数等参数，可得到选取的连接点坐标如下表所示。
表13.东向连接点的坐标
| 初始DSM编号 | 列坐标 | 行坐标 |
| NW | 10 | 323,349，……，2883 |
| NE | 无 | 无 |
| SW | 10 | 323,349，……，2883 |
| SE | 无 | 无 |
步骤S4.3，将步骤S4.2获取的南向连接点和东向连接点的像素坐标转换为地理坐标。
地理坐标的计算依然是使用的步骤S2.1中提供的R-D方程。方程中，已知参数为像素坐标，未知参数为地面点的坐标，R-D方程为三个方程，其中包含三个未知量，可通过牛顿迭代法完成。经过计算之后的连接点地理坐标分别如表14所示。
步骤S4.4，计算步骤S4.3解算的地理坐标在邻接关系矩阵中与所述初始DSM的相邻影像中的像素坐标。
使用相邻影像参数计算连接点对应的像素坐标，其计算过程中的核心方程依然是步骤S2.1中提供的R-D方程。此时方程中的已知参数为地面点坐标，未知参数为像素坐标(x,y)，由于R-D方程包含三个方程，仅有两个未知参数，因此方程可解，解算结果如表14和表15所示。
表14.南向连接点的坐标映射关系
表15.东向连接点的坐标映射关系
步骤S5，基于所述连接点在所述初始DSM中的像素坐标，确定所述连接点的经相位高程转换之后的高程值。
步骤S4仅获取了东向连接点和南向连接点的像素坐标，而没有获取连接点对应的高程，本步骤用以获取高程值，提供步骤S6的输入参数。
表1中的附加数据包含相位高程转换之后的高程图。根据东向连接点和南向连接点的像素坐标，获取相位高程转换之后的高程图中的高程值。考虑到像素坐标一般为整数，在计算过程中，采用向下取整的方式，获取坐标值。下表中给出了对应的高程数据，其中，NW南侧连接点高程来源于表14第一行，NW东侧连接点高程来源于表15第一行，NE西侧连接点高程由本步骤S5计算得来，NE南侧连接点高程来源于表14第二行，SW东侧连接点高程来源于表15第三行，SW北侧连接点高程由本步骤S5计算得来，SE西侧连接点高程由本步骤S5计算得来，SE北侧连接点高程连接点高程由本步骤S5计算得来。
表16.连接点的高程信息。
步骤S6，基于所述控制点的像素坐标和所述连接点在所述初始DSM中的像素坐标以及所述连接点的高程，构建最小二乘平差方程，对所述初始DSM进行区域网平差，确定所述初始DSM的误差多项式系数。
使用步骤S1中表1提及的控制点的高程、步骤S2提供的控制点在主影像中的像素坐标、步骤S4提供的连接点在所述初始DSM中的像素坐标，以及步骤S5获取的经相位高程转换之后的连接点高程，构建最小二乘平差方程，获取平差系数。最小二乘平差方程的表达式如下：
BX＝L (31)
其中，
X＝(a1 b1 c1…aK bK cK)′ (33)
其中，x为控制点的东西方向像素坐标，
为第k景初始DSM中的第Mk个控制点的东西方向像素坐标，
为与第一景初始DSM有邻接关系的初始DSM中第Mk+Pk个连接点的东西方向坐标，
为第L个初始DSM中第Mi+Pi个连接点的东西向坐标，
为与第L个初始DSM相邻的初始DSM中第Mi+Pi个连接点的东西向坐标，y为控制点的南北方向像素坐标，
为第k景初始DSM中的第Mk个控制点的南北方向像素坐标，
为与第一景初始DSM有邻接关系的初始DSM中第Mk+Pk个连接点的南北方向坐标，
为第L个初始DSM中第Mi+Pi个连接点的南北向坐标，
为与第L个初始DSM相邻的初始DSM中第Mi+Pi个连接点的南北向坐标，
为第L个初始DSM中第M+p个连接点的高程，
为与第L个初始DSM相邻的初始DSM中第M+p个连接点的高程，Mk为第k景初始DSM中的控制点数量，K为参与平差的初始DSM影像总数量，P为第P个连接点，Pk为第k景初始DSM影像中的连接点总数量，Ad为与处于最西北角的第一景初始DSM相邻的初始DSM，L为第L景具有连接点的初始DSM，R为第L景初始DSM的邻接初始DSM；(aK bK cK)为第K景初始DSM的待求多项式参数，
为第K景初始DSM的第M个控制点的高程误差。
上述最小二乘平差方程的最小二乘解为
X＝(BTB)-1BTL (35)
公式(32)包含了两类主要信息，首先，基线误差以及各类其他误差在影像中表现为像素坐标x与y的一阶函数，这消除了各类低频误差的误差量，确保剩余的误差中仅保留了无法消除的随机误差；其次，每个初始DSM的东向和南向邻接初始DSM的连接点高程是相等的，这确保了初始DSM拼接过程中，不会出现拼接误差。
在本实施例中，带入表7-表10对应的控制点坐标、表14-表15对应的连接点坐标，可获得X的解算结果见下表。
表17.经区域网平差计算得到的初始DSM修正系数。
| 索引 | A | B | C |
| NW | 0.000530123 | -0.000213244 | -0.000729342 |
| NE | -0.000623581 | -0.000238376 | 0.000302065 |
| SW | 0.000434656 | -0.000303871 | -0.000554629 |
| SE | -0.000896193 | -0.000342665 | 0.0000 |
DSM产品的质量是由卫星与地面应用共同保障的。其中卫星方确保指标的合理性，使得DSM的产品精度能够达到10m量级，从而为地面应用的后处理提供可靠的输入条件。地面应用负责的是卫星无控精度到有控精度的提升。提升过程中需要使用大量的地面控制或参考点。在外部数据参与平差的过程中，平面位置被认为是极为精确的。在这种前提下，仅需要考虑高程误差。值得说明的是，此时的修正对象已经不同于传统测量检校过程中的成像参数修正，而是直接针对产品进行修正。这种修正的前提是，10m的初始DSM产品误差仅仅包含低阶系统性误差，而不包含高阶误差。在满足上述条件的情况下，本发明能够将高程精度提升至1m左右。
本发明还提供了一种计算机程序，其中，当所述计算机程序由计算设备的处理器执行时，执行本文实施例所述的方法。
本发明还提供了一种计算机可读存储介质，其中，所述计算机可读存储介质上存储计算机程序，当所述计算机程序由计算设备的处理器执行时，执行本文实施例所述的方法。
以上内容仅为本发明的较佳实施例，对于本领域的普通技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，本说明书内容不应理解为对本发明的限制。
Claims (8)
1.一种面向InSAR的初始DSM的区域网平差方法，其特征在于，所述方包括以下步骤：
步骤S1，利用合成孔径雷达干涉(InSAR)生成初始DSM及附加数据，所述附加数据包括主影像头文件、从影像头文件、控制点质量评估文件、初始DSM对应的SRTM和初始DSM对应的SRTM误差；
步骤S2，构建斜距-多普勒方程和椭球方程，并利用构建的所述斜距-多普勒方程和椭球方程，计算所述控制点质量评估文件中的控制点在主影像中的像素坐标；
步骤S3，根据各个所述初始DSM的经纬度范围，确定各个所述初始DSM的邻接关系矩阵；
步骤S4，基于所述邻接关系矩阵，均匀抽取各个所述初始DSM重叠处的连接点，并计算所述连接点在所述初始DSM中的像素坐标以及与所述初始DSM相邻的影像中的像素坐标；
步骤S5，基于所述连接点在所述初始DSM中的像素坐标，确定所述连接点的经相位高程转换之后的高程值；
步骤S6，基于所述控制点的像素坐标和所述连接点在所述初始DSM中的像素坐标以及所述连接点的高程，构建最小二乘平差方程，对所述初始DSM进行区域网平差，确定所述初始DSM的误差多项式系数。
2.根据权利要求1所述的方法，其特征在于，步骤S3具体包括以下子步骤：
步骤S3.1，构建各个所述初始DSM的初始邻接关系矩阵，所述初始邻接关系矩阵F表达式为：
其中，f(i,j)为第i列第j行的初始DSM，g为东西向的初始DSM数量，k为南北向的初始DSM数量；
步骤S3.2，对步骤S3.1中构建的初始邻接关系矩阵中的所有初始DSM逐一进行迭代分析，确定各个初始DSM的邻接关系矩阵。
3.根据权利要求2所述的方法，其特征在于，步骤S3.2具体包括以下子步骤：
步骤S3.2.1，确定位于最西北角的初始DSM，并将该DSM记为f(0,0)；
步骤S3.2.2，分别确定与所述最西北角的初始DSM的东向邻接的初始DSM和南向邻接的初始DSM；
步骤S3.2.3，依次判断除最西北角初始DSM的其他初始DSM的邻接初始DSM，以确定所有初始DSM的邻接初始DSM，并最终确定初始DSM的邻接关系矩阵。
4.根据权利要求3所述的方法，其特征在于，将所述最西北角的初始DSM f(0,0)的经度范围记为[lons1,lone1]，纬度范围记为[lats1,late1]；将任一其他初始DSM f(i,j)的经度范围记为[lons2,lone2]，纬度范围记为[lats2,late2]；
其中，确定与所述最西北角的初始DSM的东向邻接的初始DSM具体包括以下子步骤：
(1)判断任一其他初始DSMf(i,j)与所述最西北角的初始DSMf(0,0)在经度方向是否有重叠，具体包括：
如果lons2>lons1，且lons2<lone1，则判定f(i,j)与f(0,0)在经度方向有重叠，并继续执行步骤(2)，否则判定无重叠，即，判定f(i,j)与f(0,0)在东西向无邻接关系，不需要执行下述步骤；
(2)计算f(i,j)与f(0,0)的纬度重叠范围；
纬度重叠范围表达为
lat_com＝min(late2,late1)-max(lats2,lats1)
其中，lat_com为纬度方向的重叠范围，min()为求取最小值，max()为求取最大值；
(3)判断f(i,j)与f(0,0)是否具备东向邻接关系；
具体的，如果lat_com>(late1-lats1)*0.7或者lat_com>(late2-lats2)*0.7，则判定f(i,j)与f(0,0)具备东向邻接关系，并赋值i＝0+1，j＝0，以寻找f(1,0)的东向和南向邻接DSM；否则，判定f(i,j)与f(0,0)不具备东向邻接关系，并确定与所述最西北角的初始DSM的南向邻接的初始DSM。
5.根据权利要求4所述的方法，其特征在于，确定与所述最西北角的初始DSM的南向邻接的初始DSM具体包括以下子步骤：
(1)判断f(i,j)与f(0,0)在纬度方向是否有重叠，具体包括：
如果lats2>lats1，且lats2<late1，则判定f(i,j)与f(0,0)在纬度方向有重叠，并继续执行步骤(2)，否则判定无重叠，即，判定f(i,j)与f(0,0)在南北向无邻接关系，不需要执行下述步骤；
(2)计算f(i,j)与f(0,0)的经度重叠范围。
经度重叠范围表达为
lon_com＝min(lone2,lone1)-max(lons2,lons1)
其中，lon_com为纬度方向的重叠范围，min()为求取最小值，max()为求取最大值；
(3)判断f(i,j)与f(0,0)是否具备南向邻接关系；
其中，如果lon_com>(lone1-lons1)*0.7或者lon_com>(lone2-lons2)*0.7，则判定f(i,j)与f(0,0)具备南向邻接关系，并赋值i＝0，j＝0+1；否则，判定f(i,j)与f(0,0)不具备南向邻接关系。
6.根据权利要求1所述的方法，其特征在于，步骤S4具体包括以下子步骤：
步骤S4.1，根据卫星成像时刻的飞行方向确定所述初始DSM在地理编码过程中的镜像变换关系；
步骤S4.2，根据通过步骤S3确定的邻接关系矩阵以及通过步骤S4.1确定的镜像变换关系，确定所述连接点在所述初始DSM中的像素坐标，其中，所述连接点包括南向连接点和东向连接点；
步骤S4.3，将步骤S4.2获取的南向连接点和东向连接点的像素坐标转换为地理坐标；
步骤S4.4，计算步骤S4.3解算的地理坐标在邻接关系矩阵中与所述初始DSM的相邻影像中的像素坐标。
7.根据权利要求1所述的方法，其特征在于，在步骤S5中，最小二乘平差方程的表达式如下：
BX＝L
其中，
X＝(a1 b1 c1…aK bK cK)′
其中，x为控制点的东西方向像素坐标，
为第k景初始DSM中的第Mk个控制点的东西方向像素坐标，
为与第一景初始DSM有邻接关系的初始DSM中第Mk+Pk个连接点的东西方向坐标，
为第L个初始DSM中第Mi+Pi个连接点的东西向坐标，
为与第L个初始DSM相邻的初始DSM中第Mi+Pi个连接点的东西向坐标，y为控制点的南北方向像素坐标，
为第k景初始DSM中的第Mk个控制点的南北方向像素坐标，
为与第一景初始DSM有邻接关系的初始DSM中第Mk+Pk个连接点的南北方向坐标，
为第L个初始DSM中第Mi+Pi个连接点的南北向坐标，
为与第L个初始DSM相邻的初始DSM中第Mi+Pi个连接点的南北向坐标，
为第L个初始DSM中第M+p个连接点的高程，
为与第L个初始DSM相邻的初始DSM中第M+p个连接点的高程，Mk为第k景初始DSM中的控制点数量，K为参与平差的初始DSM影像总数量，P为第P个连接点，Pk为第k景初始DSM影像中的连接点总数量，Ad为与处于最西北角的第一景初始DSM相邻的初始DSM，L为第L景具有连接点的初始DSM，R为第L景初始DSM的邻接初始DSM；(aK bK cK)为第K景初始DSM的待求多项式参数，
为第K景初始DSM的第M个控制点的高程误差。
8.根据权利要求7所述的方法，其特征在于，上述最小二乘平差方程的最小二乘解为
X＝(BTB)-1BTL。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110783009.4A CN113514829B (zh) | 2021-07-12 | 面向InSAR的初始DSM的区域网平差方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110783009.4A CN113514829B (zh) | 2021-07-12 | 面向InSAR的初始DSM的区域网平差方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113514829A true CN113514829A (zh) | 2021-10-19 |
| CN113514829B CN113514829B (zh) | 2024-06-25 |
Family
ID=
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114662059A (zh) * | 2022-05-25 | 2022-06-24 | 深圳市海伊石油技术有限公司 | 一种海上卫星大地坐标的高程拟合方法与装置 |
| CN114689015A (zh) * | 2021-11-29 | 2022-07-01 | 成都理工大学 | 一种提高光学卫星立体影像dsm高程精度的方法 |
| CN116188372A (zh) * | 2022-12-14 | 2023-05-30 | 自然资源部黑龙江基础地理信息中心(自然资源部黑龙江测绘资料档案馆) | 一种基于立体测绘的阻水构筑物识别方法 |
| CN117237565A (zh) * | 2023-09-27 | 2023-12-15 | 自然资源部国土卫星遥感应用中心 | 基于高分辨率卫星立体影像的建筑物白模制作方法 |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114689015A (zh) * | 2021-11-29 | 2022-07-01 | 成都理工大学 | 一种提高光学卫星立体影像dsm高程精度的方法 |
| CN114689015B (zh) * | 2021-11-29 | 2023-01-17 | 成都理工大学 | 一种提高光学卫星立体影像dsm高程精度的方法 |
| CN114662059A (zh) * | 2022-05-25 | 2022-06-24 | 深圳市海伊石油技术有限公司 | 一种海上卫星大地坐标的高程拟合方法与装置 |
| CN116188372A (zh) * | 2022-12-14 | 2023-05-30 | 自然资源部黑龙江基础地理信息中心(自然资源部黑龙江测绘资料档案馆) | 一种基于立体测绘的阻水构筑物识别方法 |
| CN117237565A (zh) * | 2023-09-27 | 2023-12-15 | 自然资源部国土卫星遥感应用中心 | 基于高分辨率卫星立体影像的建筑物白模制作方法 |
| CN117237565B (zh) * | 2023-09-27 | 2024-02-13 | 自然资源部国土卫星遥感应用中心 | 基于高分辨率卫星立体影像的建筑物白模制作方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109903352B (zh) | 一种卫星遥感影像大区域无缝正射影像制作方法 | |
| Hu et al. | Understanding the rational function model: methods and applications | |
| US9378585B2 (en) | System and method for automatic geometric correction using RPC | |
| CN106127697B (zh) | 无人机机载成像高光谱几何校正方法 | |
| US8427505B2 (en) | Geospatial modeling system for images and related methods | |
| CN102654576B (zh) | 基于sar图像和dem数据的图像配准方法 | |
| CN107144293A (zh) | 一种视频卫星面阵相机的几何定标方法 | |
| CN110046563B (zh) | 一种基于无人机点云的输电线路断面高程修正方法 | |
| CN103673995A (zh) | 一种线阵推扫式相机在轨光学畸变参数标定方法 | |
| CN111473802A (zh) | 一种基于线阵推扫的光学传感器内方位元素定标方法 | |
| CN112050725A (zh) | 一种融合InSAR和GPS的地表形变监测方法 | |
| Wu et al. | Co-registration of lunar topographic models derived from Chang’E-1, SELENE, and LRO laser altimeter data based on a novel surface matchingmethod | |
| CN108562900B (zh) | 一种基于高程校正的sar图像几何配准方法 | |
| CN114241064B (zh) | 一种遥感卫星内外方位元素实时几何定标方法 | |
| CN109029379B (zh) | 一种高精度小基高比立体测绘方法 | |
| CN110660099B (zh) | 基于神经网络的遥感影像处理的有理函数模型拟合方法 | |
| CN111156969A (zh) | 一种宽幅遥感影像立体测绘方法及系统 | |
| CN109188483B (zh) | 一种时序化高精度外方位元素自动定标方法 | |
| CN108876829B (zh) | 基于非线性尺度空间及径向基函数的sar高精度配准方法 | |
| El-Ashmawy | A comparison study between collinearity condition, coplanarity condition, and direct linear transformation (DLT) method for camera exterior orientation parameters determination | |
| CN114280608B (zh) | 一种DInSAR高程相关大气效应去除方法及系统 | |
| CN115343709B (zh) | 一种适用于分布式星载sar系统的地形反演方法 | |
| Chen et al. | Large-scale block bundle adjustment of LROC NAC images for Lunar South Pole mapping based on topographic constraint | |
| CN113514829A (zh) | 面向InSAR的初始DSM的区域网平差方法 | |
| CN111044076B (zh) | 基于参考底图的高分一号b卫星几何检校方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant |