CN110133653B - 一种基于dsm数据的星载sar图像快速间接定位方法 - Google Patents

一种基于dsm数据的星载sar图像快速间接定位方法 Download PDF

Info

Publication number
CN110133653B
CN110133653B CN201910456759.3A CN201910456759A CN110133653B CN 110133653 B CN110133653 B CN 110133653B CN 201910456759 A CN201910456759 A CN 201910456759A CN 110133653 B CN110133653 B CN 110133653B
Authority
CN
China
Prior art keywords
point
azimuth
ground
satellite
distance
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.)
Active
Application number
CN201910456759.3A
Other languages
English (en)
Other versions
CN110133653A (zh
Inventor
刘敏
肖鹏
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Academy of Space Technology CAST
Original Assignee
China Academy of Space Technology CAST
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by China Academy of Space Technology CAST filed Critical China Academy of Space Technology CAST
Priority to CN201910456759.3A priority Critical patent/CN110133653B/zh
Publication of CN110133653A publication Critical patent/CN110133653A/zh
Application granted granted Critical
Publication of CN110133653B publication Critical patent/CN110133653B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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/9004—SAR image acquisition techniques
    • 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/9004—SAR image acquisition techniques
    • G01S13/9017—SAR image acquisition techniques with time domain processing of the SAR signals in azimuth

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)
  • Signal Processing (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明涉及一种基于DSM数据的星载SAR图像快速间接定位方法、装置和计算机可读存储介质，该方法包括：获取星载SAR图像的系统参数和图像参数；拟合卫星轨道曲线；由RD模型确定SAR图像定位在地面上的四个角点坐标；四个角点坐标构成一个四边形，获取这个四边形内的DSM数据构建三维格网；对于四边形内每一个地面像素点，通过快速估计方位时刻的方式确定其对应的方位位置和距离位置，实现间接定位；该方法能够有效提升地面像素点确定多普勒中心方位时刻的效率，并且确定多普勒中心方位时刻和遍历确定的多普勒中心方位时刻是相同的，准确性好。

Description

一种基于DSM数据的星载SAR图像快速间接定位方法
技术领域
本发明涉及星载合成孔径雷达图像定位技术领域，尤其涉及一种基于DSM数据的星载SAR图像快速间接定位方法、装置和计算机可读存储介质。
背景技术
星载合成孔径雷达(Synthetic Aperture Radar，SAR)是一种搭载在卫星上并工作在微波波段的主动遥感系统。由于星载SAR不受日照、天气、地理等因素限制，可以实现全天时、全天候的对地观测。因此，星载SAR在灾害监测、环境监测、海洋观测、资源勘探、农作物调查、森林勘察、测绘和军事等领域的应用具有独到的优势。
SAR图像是一种斜距图像，在方位向上应用了多普勒处理进行了图像的方位向合成，在距离向上利用了脉冲压缩技术构建了卫星和地面目标之间几何位置关系。SAR图像中，每个像点都对应了地面上的一个或多个(地面高度不同)目标空间位置。已知地面区域的高程信息，就可以确定像点对应的目标空间位置(一个或者多个目标)。
间接定位方法可概括为：已知地面目标的三维坐标，通过计算目标点的多普勒中心时刻确定方位向的位置，计算多普勒中心时刻的卫星位置和地面目标的斜距确定在距离向的位置。间接定位的计算复杂度主要在于确定目标点的多普勒中心时刻，通常的处理步骤是把每个方位时刻都遍历后寻找一个绝对值最小的时刻作为多普勒中心时刻，但这样运算量极大，非常耗时，效率很低，不利于应用。
因此，针对以上不足，需要提供一种能够快速实现星载SAR图像间接定位的方法。
发明内容
(一)要解决的技术问题
本发明要解决的技术问题是解决现有技术星载SAR图像间接定位法确定目标点的多普勒中心时刻非常耗时，效率很低的问题。
(二)技术方案
为了解决上述技术问题，本发明提供了一种基于DSM数据的星载SAR图像快速间接定位方法，包括如下步骤：
S1、获取星载SAR图像的系统参数和图像参数；
S2、根据系统参数和图像参数中多个时刻的卫星位置和卫星速度，拟合卫星轨道曲线；
S3、根据拟合的所述卫星轨道曲线，由RD模型确定SAR图像定位在地面上的四个角点坐标，四个角点坐标为SAR图像对应的第一个方位时刻和最后一个方位时刻分别到最近斜距和最远斜距的WGS84坐标；
S4、利用所得的四个角点坐标构成一个四边形，获取这个四边形内的DSM数据构建以经度、纬度和高度三个维度构成的三维格网，根据DSM数据精度设定地面像素分辨率，每个三维格网点为一个地面像素点；
S5、对于四边形内每一个地面像素点，通过快速估计方位时刻的方式确定其对应的方位位置和距离位置；
通过所述快速估计方位时刻的方式包括利用四边形内任一地面像素点的投影点和四个角点，以及三角形几何关系确定多普勒中心的估计方位时刻；通过估计方位时刻的多普勒值和场景中心点的多普勒斜率，确定实际多普勒中心方位时刻，得到该地面像素点对应的方位位置；根据多普勒中心方位时刻的卫星位置和地面像素点确定斜距，得到该地面像素点对应的距离位置。
优选地，所述步骤S5包括：
S5-1、确定第一个方位时刻到最后一个方位时刻的平均经度变化率和平均纬度变化率；
S5-2、选取在四边形内的任意一个地面像素点，将该点和其地面投影点，以及四个角点中的第一个方位时刻和最后一个方位时刻到最近斜距的第一点和第三点由WGS84坐标分别转换至笛卡尔坐标；
S5-3、在第一点与第三点构成的直线上，通过平均经度变化率和平均纬度变化率估计出平均变化条件下距第一点的距离为地面投影点到第一点距离的估计点，并转换至笛卡尔坐标；
S5-4、由步骤S5-3估计点、第一点以及地面投影点构成的三角形，根据三角形余弦定理，计算第一点分别与步骤S5-3估计点、地面投影点构成的直线的夹角；
S5-5、由步骤S5-4确定的夹角，根据三角形几何关系，计算地面投影点到第一点的距离在第一点与步骤S5-3估计点构成的直线方向上的投影长度；
S5-6、通过比例关系确定多普勒中心估计方位时刻；
S5-7、根据拟合卫星轨道曲线确定对应的卫星位置和卫星速度，计算估计方位时刻的多普勒频率值；
S5-8、根据场景中心点的第一个方位时刻和最后一个方位时刻的多普勒频率得到多普勒斜率；
S5-9、由估计方位时刻和场景中心点的多普勒斜率计算多普勒中心方位时刻；
S5-10、由多普勒中心方位时刻的卫星位置和地面像素点计算二者之间斜距。
优选地，所述步骤S1中获取的系统参数和图像参数包括：
光速c，地球模型[Re,Rp]，其中Re为地球赤道半径，Rp为地球椭球极半径，方位向和距离向采样点数[Na,Nr]，最近斜距Rmin，距离采样频率fs，多普勒中心频率f0，SAR图像的第一个方位时刻和最后一个方位时刻分别为t1和
Figure GDA0002712246410000031
脉冲重复频率prf，SAR图像每行成像时间间隔tline，波长λ，距离分辨单元rgnres，方位分辨单元azires，图像产品第i个方位时刻ti的对应卫星位置Rs(ti)＝[Rsxi,Rsyi,Rszi]和卫星速度矢量Vs(ti)＝[Vsxi,Vsyi,Vszi]。
优选地，所述步骤S2中通过多项式拟合公式拟合卫星轨道曲线。
优选地，所述步骤S2中拟合卫星轨道曲线的多项式拟合公式为：
卫星轨道曲线的位置向量：
Figure GDA0002712246410000041
卫星轨道曲线的速度向量：
Figure GDA0002712246410000042
N为不小于3的正整数，a1n(n＝1,...,N)，b1n(n＝1,...,N)，c1n(n＝1,...,N)分别表示多项式拟合求解的Xs、Ys、Zs的和Xv、Yv、Zv各项系数；时间t的取值为ti＝t1+(i-1)·tline，其中i＝1,...,Na；Xs、Ys、Zs和Xv、Yv、Zv分别表示卫星位置和速度的三个分量。
优选地，N的取值范围为3～6。
优选地，所述步骤S3包括：
选择SAR图像对应的第一个方位时刻t1和最后一个方位时刻
Figure GDA0002712246410000049
利用RD模型进行定位迭代，计算出第一个方位时刻t1到最近斜距Rmin的WGS84坐标，坐标记为第一点C1＝(B11,L11,0)；
利用RD模型进行定位迭代，计算出第一个方位时刻t1到最远斜距Rmin+(Nr-1)·rgnres的WGS84坐标，记为第二点
Figure GDA0002712246410000043
利用RD模型进行定位迭代，计算出最后一个方位时刻
Figure GDA0002712246410000044
到最近斜距Rmin的WGS84坐标，记为第三点
Figure GDA0002712246410000045
利用RD模型进行定位迭代，计算出最后一个方位时刻
Figure GDA0002712246410000046
到最远斜距Rmin+(Nr-1)·rgnres的WGS84坐标，记为第四点
Figure GDA0002712246410000047
优选地，所述步骤S5中：
步骤S5-1中确定第一个方位时刻t1到最后一个方位时刻
Figure GDA0002712246410000048
的平均经度变化率和平均纬度变化率包括：
计算在最近斜距Rmin处，第一个方位时刻t1到最后一个方位时刻
Figure GDA00027122464100000515
的经度变化率KL1和纬度变化率KB1；
计算在最远斜距Rmin+(Nr-1)·rgnres处，第一个方位时刻t1到最后一个方位时刻
Figure GDA0002712246410000051
的经度变化率
Figure GDA0002712246410000052
和纬度变化率
Figure GDA0002712246410000053
取求得结果的平均值，得到平均经度变化率和平均纬度变化率，表达式为：
平均经度变化率：
Figure GDA0002712246410000054
平均纬度变化率：
Figure GDA0002712246410000055
步骤S5-2中选取在四边形内的任意一个地面像素点T＝(Bij,Lij,Hij)，将该点和其地面投影点Tground＝(Bij,Lij,0)，以及四个角点中的第一个方位时刻t1和最后一个方位时刻
Figure GDA0002712246410000056
到最近斜距Rmin的第一点C1＝(B11,L11,0)和第三点
Figure GDA0002712246410000057
由WGS84坐标分别转换至笛卡尔坐标包括：
转换地面像素点T＝(Bij,Lij,Hij)至笛卡尔坐标为RT＝[dx-ij,dy-ij,dz-ij]；
转换地面投影点Tground＝(Bij,Lij,0)至笛卡尔坐标为
Figure GDA0002712246410000058
转换第一点C1＝(B11,L11,0)至笛卡尔坐标为RC1＝[dx-11,dy-11,dz-11]；
转换第三点
Figure GDA0002712246410000059
至笛卡尔坐标为
Figure GDA00027122464100000510
得到第一点C1＝(B11,L11,0)到第三点
Figure GDA00027122464100000511
的距离为：
Figure GDA00027122464100000512
得到地面投影点Tground＝(Bij,Lij,0)到第一点C1＝(B11,L11,0)的距离为：
Figure GDA00027122464100000513
步骤S5-3中在第一点C1＝(B11,L11,0)与第三点
Figure GDA00027122464100000514
构成的直线上，通过平均经度变化率和平均纬度变化率估计出平均变化条件下距第一点的距离为地面投影点到第一点距离的估计点Cg＝(Bg,Lg,0)，并转换至笛卡尔坐标包括：
估计点Cg＝(Bg,Lg,0)满足如下关系式：
Figure GDA0002712246410000061
Figure GDA0002712246410000062
转换估计点Cg＝(Bg,Lg,0)至笛卡尔坐标为
Figure GDA0002712246410000063
得到第一点C1＝(B11,L11,0)到Cg＝(Bg,Lg,0)的距离为：
Figure GDA0002712246410000064
得到地面投影点Tground＝(Bij,Lij,0)到Cg＝(Bg,Lg,0)的距离为：
Figure GDA0002712246410000065
步骤S5-4中由步骤S5-3估计点Cg＝(Bg,Lg,0)、第一点C1＝(B11,L11,0)以及地面投影点Tground＝(Bij,Lij,0)构成的三角形，根据三角形余弦定理，计算第一点C1＝(B11,L11,0)分别与步骤S5-3估计点Cg＝(Bg,Lg,0)、地面投影点Tground＝(Bij,Lij,0)构成的直线的夹角时，三角形TgroundC1Cg的∠TgroundC1Cg记为α，
Figure GDA0002712246410000066
步骤S5-5中由步骤S5-4确定的夹角，根据三角形几何关系，计算地面投影点Tground＝(Bij,Lij,0)到第一点C1＝(B11,L11,0)的距离
Figure GDA0002712246410000067
在第一点C1＝(B11,L11,0)与步骤S5-3估计点Cg＝(Bg,Lg,0)构成的直线C1Cg方向上的投影长度时，
Figure GDA0002712246410000068
投影到直线C1Cg方向上的长度记为
Figure GDA0002712246410000069
表达式为
Figure GDA00027122464100000610
步骤S5-6中通过比例关系确定多普勒中心估计方位时刻ti_est时，表达式为：
Figure GDA00027122464100000611
其中，
Figure GDA00027122464100000612
表示对x向下取整数；
步骤S5-7中根据拟合卫星轨道曲线确定对应的卫星位置和卫星速度，计算估计方位时刻ti_est的多普勒频率值fD(ti_est)时，表达式为：
Figure GDA0002712246410000071
其中，Vs(ti_est)＝[Vsxi_est,Vsyi_est,Vszi_est]，Rs(ti_est)＝[Rsxi_est,Rsyi_est,Rszi_est]，RT＝[dx-ij,dy-ij,dz-ij]，R＝|Rs(ti_est)-RT|；
步骤S5-8中根据场景中心点的第一个方位时刻t1和最后一个方位时刻
Figure GDA0002712246410000078
的多普勒频率得到多普勒斜率时，第一个方位时刻t1和最后一个方位时刻
Figure GDA0002712246410000079
的多普勒频率分别记为fD_center(t1)、
Figure GDA00027122464100000710
多普勒斜率Kf表达式为：
Figure GDA0002712246410000072
步骤S5-9中由估计方位时刻ti_est和场景中心点的多普勒斜率计算多普勒中心方位时刻tI时，表达式为：
Figure GDA0002712246410000073
其中，
Figure GDA0002712246410000074
表示对x向下取整数；I表示SAR图像像元的行号；
步骤S5-10中由多普勒中心方位时刻的卫星位置Rs_I＝[RsxI,RsyI,RszI]和地面像素点T＝(dx-ij,dy-ij,dz-ij)计算二者之间斜距RJ时，表达式为：
Figure GDA0002712246410000075
Figure GDA0002712246410000076
其中，
Figure GDA0002712246410000077
表示对x向下取整数；J表示SAR图像像元的列号。
本发明还提供了一种基于DSM数据的星载SAR图像快速间接定位装置，所述装置包括：
系统参数和图像参数获取模块，用于获取星载SAR图像的系统参数和图像参数；
轨道拟合模块，用于根据系统参数和图像参数中多个时刻的卫星位置和卫星速度，拟合卫星轨道曲线；
角点确定模块，用于根据拟合的所述卫星轨道曲线，由RD模型确定SAR图像定位在地面上的四个角点坐标，四个角点坐标为SAR图像对应的第一个方位时刻和最后一个方位时刻分别到最近斜距和最远斜距的WGS84坐标；
三维格网模块，用于利用所得的四个角点坐标构成一个四边形，获取这个四边形内的DSM数据构建以经度、纬度和高度三个维度构成的三维格网，根据DSM数据精度设定地面像素分辨率，每个三维格网点为一个地面像素点；
间接定位模块，用于对于四边形内每一个地面像素点，通过快速估计方位时刻的方式确定其对应的方位位置和距离位置；通过所述快速估计方位时刻的方式包括利用四边形内任一地面像素点的投影点和四个角点，以及三角形几何关系确定多普勒中心的估计方位时刻；通过估计方位时刻的多普勒值和场景中心点的多普勒斜率，确定实际多普勒中心方位时刻，得到该地面像素点对应的方位位置；根据多普勒中心方位时刻的卫星位置和地面像素点确定斜距，得到该地面像素点对应的距离位置。
本发明还提供了一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现上述任一项所述的方法的步骤。
(三)有益效果
本发明的上述技术方案具有如下优点：本发明提供了一种基于DSM数据的星载SAR图像快速间接定位方法，该方法首先通过正向定位方法获取星载合成孔径雷达图像的四个角点位置(高程为0)；然后在四个角点构成的四边形内获取DSM数据建立以经度、纬度和高度三个元素构成的三维格网；通过几何关系估算确定估计方位时刻；再通过估计方位时刻的多普勒值和场景中心点的多普勒斜率计算多普勒中心方位时刻，对每个三维格网点计算其对应的方位位置(多普勒中心时刻)和距离位置，从而实现了对星载合成孔径雷达图像的快速间接定位。该方法能够有效提升地面像素点确定多普勒中心方位时刻的效率，并且确定多普勒中心方位时刻和遍历确定的多普勒中心方位时刻是相同的，准确性好。
附图说明
图1是本发明实施例中的星载SAR图像快速间接定位方法流程图；
图2是本发明实施例中采用的北京西山附近区域的TanDEM-X卫星SAR图像；
图3是本发明实施例中采用的北京西山附近区域的DSM；
图4(a)是本发明实施例中北京西山附近区域的TanDEM-X卫星SAR图像定位区域的DSM；
图4(b)是本发明实施例中北京西山附近区域的TanDEM-X卫星SAR图像定位结果立体图；
图4(c)是本发明实施例中北京西山附近区域的TanDEM-X卫星SAR图像定位结果平面图；
图5(a)是本发明实施例中北京西山附近区域的TanDEM-X卫星SAR图像定位结果立体图(西山半山亭)；
图5(b)是本发明实施例中北京西山附近区域的TanDEM-X卫星SAR图像定位结果平面图(西山半山亭)；
图5(c)是图5(b)中白色方框区域放大图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
如图1所示，本发明实施例提供的基于DSM数据的星载SAR图像快速间接定位方法，具体包括如下步骤：
S1、获取星载SAR图像的系统参数和图像参数。
优选地，步骤S1中获取的SAR图像的系统参数和图像参数包括：光速c，地球模型[Re,Rp]，其中Re为地球赤道半径，Rp为地球椭球极半径，方位向和距离向采样点数[Na,Nr]，最近斜距Rmin，距离采样频率fs，多普勒中心频率f0，SAR图像的第一个方位时刻和最后一个方位时刻分别为t1和
Figure GDA0002712246410000101
脉冲重复频率prf，SAR图像每行成像时间间隔tline，波长λ，距离分辨单元rgnres，方位分辨单元azires，图像产品第i个方位时刻ti的对应卫星位置Rs(ti)＝[Rsxi,Rsyi,Rszi]和卫星速度矢量Vs(ti)＝[Vsxi,Vsyi,Vszi]。
S2、根据SAR图像的系统参数和图像参数中多个时刻的卫星位置和卫星速度，拟合卫星轨道曲线，得到与SAR图像对应的方位采样时刻的卫星位置和速度。
优选地，步骤S2中通过多项式拟合公式拟合卫星轨道曲线以确定任意方位时刻的卫星位置和速度。进一步地，拟合卫星轨道曲线的多项式拟合公式为：
卫星轨道曲线的位置向量：
Figure GDA0002712246410000102
卫星轨道曲线的速度向量：
Figure GDA0002712246410000103
N为不小于3的正整数。一般情况N优选为3～6，当然也可根据位置和速度的数据精度选择合适的拟合次数。a1n(n＝1,...,N)，b1n(n＝1,...,N)，c1n(n＝1,...,N)分别表示多项式拟合求解的Xs、Ys、Zs和Xv、Yv、Zv各项系数；时间t的取值为ti＝t1+(i-1)·tline，其中i＝1,...,Na；Xs、Ys、Zs和Xv、Yv、Zv分别表示卫星位置和速度的三个分量。
S3、根据拟合的所述卫星轨道曲线，设地面平坦(高程H＝0)，由RD(距离多普勒)模型确定SAR图像定位在地面上的四个角点坐标，四个角点坐标为SAR图像对应的第一个方位时刻t1和最后一个方位时刻
Figure GDA0002712246410000111
分别到最近斜距和最远斜距的WGS84坐标(World geodeticsystem-1984Coordinate System，世界大地测量系统，简称WGS)，即四个角点坐标分别为SAR图像对应的第一个方位时刻到最近斜距的第一点和到最远斜距的第二点的WGS84坐标，以及SAR图像对应的最后一个方位时刻到最近斜距的第三点和到最远斜距的第四点的WGS84坐标。
优选地，步骤S3确定四个角点坐标具体包括：
选择SAR图像对应的第一个方位时刻t1和最后一个方位时刻
Figure GDA0002712246410000112
利用RD模型进行定位迭代，计算出第一个方位时刻t1到最近斜距Rmin的WGS84坐标，坐标记为第一点C1＝(B11,L11,0)；利用RD模型进行定位迭代，计算出第一个方位时刻t1到最远斜距Rmin+(Nr-1)·rgnres的WGS84坐标，记为第二点
Figure GDA0002712246410000113
利用RD模型进行定位迭代，计算出最后一个方位时刻
Figure GDA0002712246410000114
到最近斜距Rmin的WGS84坐标，记为第三点
Figure GDA0002712246410000115
利用RD模型进行定位迭代，计算出最后一个方位时刻
Figure GDA0002712246410000116
到最远斜距Rmin+(Nr-1)·rgnres的WGS84坐标，记为第四点
Figure GDA0002712246410000117
S4、利用步骤S3所得的四个角点坐标构成一个四边形，获取这个四边形内的DSM数据构建以经度、纬度和高度三个维度构成的三维格网，根据DSM数据精度设定地面像素分辨率，每个三维格网点为一个地面像素点。
S5、对四边形内的任意地面像素进行间接快速定位：对于四边形内每一个地面像素点，通过快速估计方位时刻的方式确定其对应的方位位置和距离位置，进而实现快速间接定位。
其中，通过所述快速估计方位时刻的方式包括：
利用四边形内任一地面像素点的投影点和四个角点，以及三角形几何关系确定多普勒中心估计方位时刻；
通过估计方位时刻的多普勒值和场景中心点的多普勒斜率，确定实际多普勒中心方位时刻，得到该地面像素点对应的方位位置；
根据多普勒中心方位时刻的卫星位置和地面像素点确定斜距，得到该地面像素点对应的距离位置。
优选地，步骤S5包括：
S5-1、确定第一个方位时刻t1到最后一个方位时刻
Figure GDA00027122464100001211
的平均经度变化率和平均纬度变化率。
优选地，首先计算在最近斜距Rmin处，第一个方位时刻t1到最后一个方位时刻
Figure GDA00027122464100001210
的经度变化率KL1和纬度变化率KB1，相应的表达式为：
经度变化率：
Figure GDA0002712246410000121
纬度变化率：
Figure GDA0002712246410000122
其次计算在最远斜距Rmin+(Nr-1)·rgnres处，第一个方位时刻t1到最后一个方位时刻
Figure GDA0002712246410000123
的经度变化率
Figure GDA0002712246410000124
和纬度变化率
Figure GDA0002712246410000125
相应的表达式为：
经度变化率：
Figure GDA0002712246410000126
纬度变化率：
Figure GDA0002712246410000127
最后取求得结果的平均值，得到平均经度变化率和平均纬度变化率，相应的表达式为：
平均经度变化率：
Figure GDA0002712246410000128
平均纬度变化率：
Figure GDA0002712246410000129
S5-2、选取在四边形内的任意一个地面像素点T＝(Bij,Lij,Hij)，将该地面像素点T＝(Bij,Lij,Hij)和其对应的地面投影点(该地面像素点对应的地面投影点WGS84坐标记为Tground＝(Bij,Lij,0))，以及四个角点中的第一个方位时刻t1和最后一个方位时刻
Figure GDA0002712246410000131
到最近斜距Rmin的第一点和第三点(四个角点中的第一个方位时刻t1到最近斜距Rmin的WGS84坐标记为第一点C1＝(B11,L11,0)，和最后一个方位时刻
Figure GDA0002712246410000132
到最近斜距Rmin的WGS84坐标记为第三点
Figure GDA0002712246410000133
由WGS84坐标分别转换至笛卡尔坐标。
优选地，计算地面投影点Tground＝(Bij,Lij,0)到C1＝(B11,L11,0)的距离
Figure GDA0002712246410000134
计算C1＝(B11,L11,0)到
Figure GDA0002712246410000135
的距离DC1-C3。
优选地，WGS84坐标(B,L,H)转换为笛卡尔坐标(dx,dy,dz)的关系式为：
卯酉圈半径：
Figure GDA0002712246410000136
其中，偏心率
Figure GDA0002712246410000137
Figure GDA0002712246410000138
将地面像素点T＝(Bij,Lij,Hij)、其地面投影点Tground＝(Bij,Lij,0)、第一点C1＝(B11,L11,0)和第三点
Figure GDA0002712246410000139
的WGS84球坐标分别转换为笛卡尔坐标RT、
Figure GDA00027122464100001310
RC1和RC3，那么有：
转换地面像素点T＝(Bij,Lij,Hij)点坐标至笛卡尔坐标为RT＝[dx-ij,dy-ij,dz-ij]，则有：
Figure GDA00027122464100001311
Figure GDA00027122464100001312
转换地面投影点Tground＝(Bij,Lij,0)点坐标至笛卡尔坐标为
Figure GDA0002712246410000141
则有：
Figure GDA0002712246410000142
Figure GDA0002712246410000143
转换第一点C1＝(B11,L11,0)点坐标至笛卡尔坐标为RC1＝[dx-11,dy-11,dz-11]，则有：
Figure GDA0002712246410000144
Figure GDA0002712246410000145
转换第三点
Figure GDA0002712246410000146
点坐标至笛卡尔坐标为
Figure GDA0002712246410000147
则有：
Figure GDA0002712246410000148
Figure GDA0002712246410000149
进而可得到第一点C1＝(B11,L11,0)到第三点
Figure GDA00027122464100001410
的距离为：
Figure GDA00027122464100001411
得到地面投影点Tground＝(Bij,Lij,0)到第一点C1＝(B11,L11,0)的距离为：
Figure GDA00027122464100001412
S5-3、在第一点C1＝(B11,L11,0)与第三点
Figure GDA00027122464100001413
构成的直线上，通过步骤S5-1求出的平均经度变化率和平均纬度变化率估计出平均变化条件下距第一点的距离为地面投影点到第一点距离的估计点Cg＝(Bg,Lg,0)，并转换至笛卡尔坐标。
优选地，在直线C1C3上，估计点Cg＝(Bg,Lg,0)满足如下关系式：
Figure GDA0002712246410000151
Figure GDA0002712246410000152
转换Cg＝(Bg,Lg,0)点坐标至笛卡尔坐标为
Figure GDA0002712246410000153
则有：
Figure GDA0002712246410000154
Figure GDA0002712246410000155
得到第一点C1＝(B11,L11,0)到Cg＝(Bg,Lg,0)的距离为：
Figure GDA0002712246410000156
得到地面投影点Tground＝(Bij,Lij,0)到Cg＝(Bg,Lg,0)的距离为：
Figure GDA0002712246410000157
S5-4、由步骤S5-3估计点Cg＝(Bg,Lg,0)、第一点C1＝(B11,L11,0)以及地面投影点Tgroud＝(Bij,Lij,0)构成的三角形，根据三角形余弦定理，计算第一点C1＝(B11,L11,0)分别与步骤S5-3估计点Cg＝(Bg,Lg,0)、地面投影点Tground＝(Bij,Lij,0)构成的直线的夹角。
优选地，根据三角形的余弦定理，计算三角形TgroundC1Cg的∠TgroundC1Cg记为α，则表达式为：
Figure GDA0002712246410000158
S5-5、由步骤S5-4确定的夹角，根据三角形几何关系，计算地面投影点Tground＝(Bij,Lij,0)到第一点C1＝(B11,L11,0)的距离
Figure GDA0002712246410000159
在第一点C1＝(B11,L11,0)与步骤S5-3估计点Cg＝(Bg,Lg,0)构成的直线C1Cg方向上的投影长度。
优选地，地面投影点Tground＝(Bij,Lij,0)到第一点C1＝(B11,L11,0)的距离
Figure GDA0002712246410000161
投影到直线C1Cg方向上的长度记为
Figure GDA0002712246410000162
表达式为：
Figure GDA0002712246410000163
S5-6、通过比例关系确定多普勒中心估计方位时刻。
估计地面像素点的方位时刻，记为ti_est，表达式：
Figure GDA0002712246410000164
这里，
Figure GDA0002712246410000165
表示对x向下取整数。
以上步骤S5-1至S5-6对应通过所述快速估计方位时刻的方式包括利用四边形内任一地面像素点的投影点和四个角点，以及三角形几何关系确定多普勒中心的估计方位时刻。
S5-7、根据拟合卫星轨道曲线确定对应的卫星位置和卫星速度，计算估计方位时刻ti_est的多普勒频率值。
估计方位时刻ti_est的多普勒频率记为fD(ti_est)，表达式为：
Figure GDA0002712246410000166
这里有，Vs(ti_est)＝[Vsxi_est,Vsyi_est,Vszi_est]，Rs(ti_est)＝[Rsxi_est,Rsyi_est,Rszi_est]，RT＝[dx-ij,dy-ij,dz-ij]，R＝|Rs(ti_est)-RT|。
S5-8、根据场景中心点的第一个方位时刻t1和最后一个方位时刻
Figure GDA0002712246410000167
的多普勒频率得到多普勒斜率。
类似地，可以得到场景中心点的第一个方位时刻t1和最后一个方位时刻
Figure GDA0002712246410000168
的多普勒频率，分别记为fD_center(t1)、
Figure GDA0002712246410000169
那么多普勒斜率可以记为：
Figure GDA00027122464100001610
以上步骤S5-7至S5-8对应通过估计方位时刻的多普勒值和场景中心点的多普勒斜率，确定实际多普勒中心方位时刻，得到该地面像素点对应的方位位置。
S5-9、由估计方位时刻ti_est和场景中心点的多普勒斜率计算多普勒中心方位时刻tI。
Figure GDA0002712246410000171
这里，
Figure GDA0002712246410000172
表示对x向下取整数。tI即为多普勒中心方位时刻，I记为SAR图像像元的行号(方位向)。
那么，多普勒中心方位时刻tI与估计方位时刻ti_est之间的误差dti＝ti_est-tI就被称作估计误差。
S5-10、由多普勒中心方位时刻tI的卫星位置Rs_I＝[RsxI,RsyI,RszI]和地面像素点T＝(dx-ij,dy-ij,dz-ij)计算它们之间距离DT-R，即为斜距RJ。那么有：
Figure GDA0002712246410000173
Figure GDA0002712246410000174
这里，
Figure GDA0002712246410000175
表示对x向下取整数。J表示SAR图像像元的列号(距离向)。
以上步骤S5-9至S5-10对应根据多普勒中心方位时刻的卫星位置和地面像素点确定斜距，得到该地面像素点对应的距离位置。
对四边形内的所有地面像素按照上述快速估计方位时刻的方式进行计算，即可实现间接快速定位。
在星载合成孔径雷达图像定位方法中，应用本发明的方法具有如下优点：
(1)正向定位确定地面四个角点，对四角点构成的四边形内的地面像点，利用几何关系和三角形余弦定理等直接确定估计方位时刻，再通过估计方位时刻的多普勒值和场景中心点的多普勒斜率计算多普勒中心方位时刻，提升地面像点确定多普勒中心方位时刻的效率。
(2)与遍历所有方位时刻相比，本发明确定多普勒中心方位时刻和遍历确定的多普勒中心方位时刻是相同的。
(3)在入射角度为θ时，高程数据的高度误差为H。采用直接定位的方法处理后，高程误差带来的水平定位的误差为H·cot θ；而用间接定位方法处理后，高程误差带来的水平定位误差约为H·cosθ。通过比较，间接定位误差约为直接定位误差的sinθ。
进一步地，在一个具体的实施方式中，为了说明本发明所提供的方法的有效性，本发明还采用了现有技术中较为成熟的星载SAR图像数据产品进行处理，相同之处不再赘述，不同之处在于：采用德国的雷达卫星TanDEM-X卫星在北京西山附近区域的SAR影像以及该地区的DSM数据对SAR图像中所有像点进行快速间接定位。其中，TanDEM-X卫星的SAR影像参数见如下表1。
表1 TanDEM-X卫星的SAR影像参数
参数名称 数值
SAR图像大小 6235×10808
方位向间隔 0.8705m
距离向间隔 0.4547m
第一个像素的斜距 706315.7084m
脉冲重复频率 3720.4126Hz
SAR图像每行成像时间间隔 1.234566E-4s
雷达频率 9649.998153MHz
卫星状态矢量 12
图像获取时间(UTC) Dec-20-2018
获取模式 聚束
升轨/降轨 升轨
视线方向 右视
请参阅图2至图5(c)，其中：图2是北京西山附近区域的TanDEM-X卫星SAR图像。图3是北京西山附近区域的DSM，其网格大小为0.69m×0.69m。图4(a)至图4(c)是北京西山附近区域的TanDEM-X卫星SAR图像定位结果。其中，图4(a)是SAR图像定位区域的DSM；图4(b)是SAR图像定位结果立体图；图4(c)是SAR图像定位结果平面图。图5(a)至图5(c)是北京西山附近区域的TanDEM-X卫星SAR图像定位结果(西山半山亭)。其中图5(a)是SAR图像定位结果立体图(西山半山亭)；图5(b)是SAR图像定位结果平面图(西山半山亭)；图5(c)是SAR图像定位结果平面图(西山半山亭)的白色方框区域放大图。
表2给出了北京西山附近区域的TanDEM-X卫星SAR图像定位结果(西山半山亭)，这里采用千寻位置公司的FindCM服务对西山半山亭的八个角点(P1-P8)位置进行多次实测，如表2所示对检查点P1的千寻实测数据标记为P1-QX，并认为实测数据为真值。M1表示直接采用获取的DSM数据进行快速间接定位处理的结果。用一个点测出一个水平面的实际高度，将其与DSM该位置的高度进行对比，这个两个水平高度差值即为校正高度值，该地区的校正DSM高度为8.5米，这里M2就表示DSM数据校正高度后快速间接定位处理的结果。M3表示DSM数据校正高度和校正电离层延迟和对流层延迟后的快速间接定位处理的结果。从表2可以看出，八个检查点M1处理结果的平均三维定位误差约为14.4米，八个检查点M2处理结果的平均三维定位误差约为4.8米，八个检查点M3处理结果的平均三维定位误差约为0.57米。
表2 北京西山附近区域的TanDEM-X卫星SAR图像定位结果(西山半山亭)
Figure GDA0002712246410000191
Figure GDA0002712246410000201
综上所述，本发明提出了一种基于DSM数据的星载合成孔径雷达图像的快速间接定位方法，该方法利用几何关系确定估计方位时刻，通过估计方位时刻快速确定多普勒中心方位时刻，从而实现了对星载合成孔径雷达图像的快速间接定位。该方法可有效解决遍历每个方位时刻以寻找多普勒中心时刻耗时长、效率低的问题。
在一些优选的实施方式中，提供了一种基于DSM数据的星载SAR图像快速间接定位装置，装置包括：
系统参数和图像参数获取模块，用于获取星载SAR图像的系统参数和图像参数；
轨道拟合模块，用于根据系统参数和图像参数中多个时刻的卫星位置和卫星速度，拟合卫星轨道曲线；
角点确定模块，用于根据拟合的所述卫星轨道曲线，由RD模型确定SAR图像定位在地面上的四个角点坐标，四个角点坐标为SAR图像对应的第一个方位时刻和最后一个方位时刻分别到最近斜距和最远斜距的WGS84坐标；
三维格网模块，用于利用所得的四个角点坐标构成一个四边形，获取这个四边形内的DSM数据构建以经度、纬度和高度三个维度构成的三维格网，根据DSM数据精度设定地面像素分辨率，每个三维格网点为一个地面像素点；
间接定位模块，用于对于四边形内每一个地面像素点，通过快速估计方位时刻的方式确定其对应的方位位置和距离位置；通过所述快速估计方位时刻的方式包括利用四边形内任一地面像素点的投影点和四个角点，以及三角形几何关系确定多普勒中心的估计方位时刻；通过估计方位时刻的多普勒值和场景中心点的多普勒斜率，确定实际多普勒中心方位时刻，得到该地面像素点对应的方位位置；根据多普勒中心方位时刻的卫星位置和地面像素点确定斜距，得到该地面像素点对应的距离位置。
本领域的技术人员可以清楚地了解到，为描述的方便和简洁，仅以上述各功能模块的划分进行举例说明，实际应用中，可以根据需要而将上述功能分配由不同的功能模块完成，即将装置的内部结构划分成不同的功能模块，以完成以上描述的全部或者部分功能。上述描述功能模块的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。
在一些优选的实施方式中，提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述任一实施方式中所述的基于DSM数据的星载SAR图像快速间接定位方法。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (10)

1.一种基于DSM数据的星载SAR图像快速间接定位方法，其特征在于，包括如下步骤：
S1、获取星载SAR图像的系统参数和图像参数；
S2、根据系统参数和图像参数中多个时刻的卫星位置和卫星速度，拟合卫星轨道曲线；
S3、根据拟合的所述卫星轨道曲线，由RD模型确定SAR图像定位在地面上的四个角点坐标，四个角点坐标为SAR图像对应的第一个方位时刻和最后一个方位时刻分别到最近斜距和最远斜距的WGS84坐标；
S4、利用所得的四个角点坐标构成一个四边形，获取这个四边形内的DSM数据构建以经度、纬度和高度三个维度构成的三维格网，根据DSM数据精度设定地面像素分辨率，每个三维格网点为一个地面像素点；
S5、对于四边形内每一个地面像素点，通过快速估计方位时刻的方式确定其对应的方位位置和距离位置；
通过所述快速估计方位时刻的方式包括利用四边形内任一地面像素点的投影点和四个角点，以及三角形几何关系确定多普勒中心的估计方位时刻；通过估计方位时刻的多普勒值和场景中心点的多普勒斜率，确定实际多普勒中心方位时刻，得到该地面像素点对应的方位位置；根据多普勒中心方位时刻的卫星位置和地面像素点确定斜距，得到该地面像素点对应的距离位置。
2.根据权利要求1所述的基于DSM数据的星载SAR图像快速间接定位方法，其特征在于，所述步骤S5包括：
S5-1、确定第一个方位时刻到最后一个方位时刻的平均经度变化率和平均纬度变化率；
S5-2、选取在四边形内的任意一个地面像素点，将该点和其地面投影点，以及四个角点中的第一个方位时刻和最后一个方位时刻到最近斜距的第一点和第三点由WGS84坐标分别转换至笛卡尔坐标；
S5-3、在第一点与第三点构成的直线上，通过平均经度变化率和平均纬度变化率估计出平均变化条件下距第一点的距离为地面投影点到第一点距离的估计点，并转换至笛卡尔坐标；
S5-4、由步骤S5-3估计点、第一点以及地面投影点构成的三角形，根据三角形余弦定理，计算第一点分别与步骤S5-3估计点、地面投影点构成的直线的夹角；
S5-5、由步骤S5-4确定的夹角，根据三角形几何关系，计算地面投影点到第一点的距离在第一点与步骤S5-3估计点构成的直线方向上的投影长度；
S5-6、通过比例关系确定多普勒中心估计方位时刻；
S5-7、根据拟合卫星轨道曲线确定对应的卫星位置和卫星速度，计算估计方位时刻的多普勒频率值；
S5-8、根据场景中心点的第一个方位时刻和最后一个方位时刻的多普勒频率得到多普勒斜率；
S5-9、由估计方位时刻和场景中心点的多普勒斜率计算多普勒中心方位时刻；
S5-10、由多普勒中心方位时刻的卫星位置和地面像素点计算二者之间斜距。
3.根据权利要求2所述的基于DSM数据的星载SAR图像快速间接定位方法，其特征在于，所述步骤S1中获取的系统参数和图像参数包括：
光速c，地球模型[Re,Rp]，其中Re为地球赤道半径，Rp为地球椭球极半径，方位向和距离向采样点数[Na,Nr]，最近斜距Rmin，距离采样频率fs，多普勒中心频率f0，SAR图像的第一个方位时刻和最后一个方位时刻分别为t1和
Figure FDA0002670640010000021
脉冲重复频率prf，SAR图像每行成像时间间隔tline，波长λ，距离分辨单元rgnres，方位分辨单元azires，图像产品第i个方位时刻ti的对应卫星位置Rs(ti)＝[Rsxi,Rsyi,Rszi]和卫星速度矢量Vs(ti)＝[Vsxi,Vsyi,Vszi]。
4.根据权利要求1所述的基于DSM数据的星载SAR图像快速间接定位方法，其特征在于：所述步骤S2中通过多项式拟合公式拟合卫星轨道曲线。
5.根据权利要求4所述的基于DSM数据的星载SAR图像快速间接定位方法，其特征在于：所述步骤S2中拟合卫星轨道曲线的多项式拟合公式为：
卫星轨道曲线的位置向量：
Figure FDA0002670640010000031
卫星轨道曲线的速度向量：
Figure FDA0002670640010000032
N为不小于3的正整数，a1n(n＝1,...,N)，b1n(n＝1,...,N)，c1n(n＝1,...,N)分别表示多项式拟合求解的Xs、Ys、Zs的和Xv、Yv、Zv各项系数；时间t的取值为ti＝t1+(i-1)·tline，其中i＝1,...,Na；Xs、Ys、Zs和Xv、Yv、Zv分别表示卫星位置和速度的三个分量，Na为方位向采样点数，SAR图像的第一个方位时刻为t1，第i个方位时刻为ti，SAR图像每行成像时间间隔tline。
6.根据权利要求5所述的基于DSM数据的星载SAR图像快速间接定位方法，其特征在于：N的取值范围为3～6。
7.根据权利要求3所述的基于DSM数据的星载SAR图像快速间接定位方法，其特征在于，所述步骤S3包括：
选择SAR图像对应的第一个方位时刻t1和最后一个方位时刻
Figure FDA0002670640010000033
利用RD模型进行定位迭代，计算出第一个方位时刻t1到最近斜距Rmin的WGS84坐标，坐标记为第一点C1＝(B11,L11,0)；
利用RD模型进行定位迭代，计算出第一个方位时刻t1到最远斜距Rmin+(Nr-1)·rgnres的WGS84坐标，记为第二点
Figure FDA0002670640010000041
利用RD模型进行定位迭代，计算出最后一个方位时刻
Figure FDA0002670640010000042
到最近斜距Rmin的WGS84坐标，记为第三点
Figure FDA0002670640010000043
利用RD模型进行定位迭代，计算出最后一个方位时刻
Figure FDA0002670640010000044
到最远斜距Rmin+(Nr-1)·rgnres的WGS84坐标，记为第四点
Figure FDA0002670640010000045
8.根据权利要求7所述的基于DSM数据的星载SAR图像快速间接定位方法，其特征在于，所述步骤S5中：
步骤S5-1中确定第一个方位时刻t1到最后一个方位时刻
Figure FDA0002670640010000046
的平均经度变化率和平均纬度变化率包括：
计算在最近斜距Rmin处，第一个方位时刻t1到最后一个方位时刻
Figure FDA0002670640010000047
的经度变化率KL1和纬度变化率KB1；
计算在最远斜距Rmin+(Nr-1)·rgnres处，第一个方位时刻t1到最后一个方位时刻
Figure FDA0002670640010000048
的经度变化率
Figure FDA0002670640010000049
和纬度变化率
Figure FDA00026706400100000410
取求得结果的平均值，得到平均经度变化率和平均纬度变化率，表达式为：
平均经度变化率：
Figure FDA00026706400100000411
平均纬度变化率：
Figure FDA00026706400100000412
步骤S5-2中选取在四边形内的任意一个地面像素点T＝(Bij,Lij,Hij)，将该点和其地面投影点Tground＝(Bij,Lij,0)，以及四个角点中的第一个方位时刻t1和最后一个方位时刻
Figure FDA00026706400100000413
到最近斜距Rmin的第一点C1＝(B11,L11,0)和第三点
Figure FDA00026706400100000414
由WGS84坐标分别转换至笛卡尔坐标包括：
转换地面像素点T＝(Bij,Lij,Hij)至笛卡尔坐标为RT＝[dx-ij,dy-ij,dz-ij]；
转换地面投影点Tground＝(Bij,Lij,0)至笛卡尔坐标为
Figure FDA00026706400100000415
转换第一点C1＝(B11,L11,0)至笛卡尔坐标为RC1＝[dx-11,dy-11,dz-11]；
转换第三点
Figure FDA0002670640010000051
至笛卡尔坐标为
Figure FDA0002670640010000052
得到第一点C1＝(B11,L11,0)到第三点
Figure FDA0002670640010000053
的距离为：
Figure FDA0002670640010000054
得到地面投影点Tground＝(Bij,Lij,0)到第一点C1＝(B11,L11,0)的距离为：
Figure FDA0002670640010000055
步骤S5-3中在第一点C1＝(B11,L11,0)与第三点
Figure FDA0002670640010000056
构成的直线上，通过平均经度变化率和平均纬度变化率估计出平均变化条件下距第一点的距离为地面投影点到第一点距离的估计点Cg＝(Bg,Lg,0)，并转换至笛卡尔坐标包括：
估计点Cg＝(Bg,Lg,0)满足如下关系式：
Figure FDA0002670640010000057
Figure FDA0002670640010000058
转换估计点Cg＝(Bg,Lg,0)至笛卡尔坐标为
Figure FDA0002670640010000059
得到第一点C1＝(B11,L11,0)到Cg＝(Bg,Lg,0)的距离为：
Figure FDA00026706400100000510
得到地面投影点Tground＝(Bij,Lij,0)到Cg＝(Bg,Lg,0)的距离为：
Figure FDA00026706400100000511
步骤S5-4中由步骤S5-3估计点Cg＝(Bg,Lg,0)、第一点C1＝(B11,L11,0)以及地面投影点Tground＝(Bij,Lij,0)构成的三角形，根据三角形余弦定理，计算第一点C1＝(B11,L11,0)分别与步骤S5-3估计点Cg＝(Bg,Lg,0)、地面投影点Tground＝(Bij,Lij,0)构成的直线的夹角时，三角形TgroundC1Cg的∠TgroundC1Cg记为α，
Figure FDA00026706400100000512
步骤S5-5中由步骤S5-4确定的夹角，根据三角形几何关系，计算地面投影点Tground＝(Bij,Lij,0)到第一点C1＝(B11,L11,0)的距离
Figure FDA0002670640010000061
在第一点C1＝(B11,L11,0)与步骤S5-3估计点Cg＝(Bg,Lg,0)构成的直线C1Cg方向上的投影长度时，
Figure FDA0002670640010000062
投影到直线C1Cg方向上的长度记为
Figure FDA0002670640010000063
表达式为
Figure FDA0002670640010000064
步骤S5-6中通过比例关系确定多普勒中心估计方位时刻ti_est时，表达式为：
Figure FDA0002670640010000065
其中，
Figure FDA0002670640010000066
表示对x向下取整数；
步骤S5-7中根据拟合卫星轨道曲线确定对应的卫星位置和卫星速度，计算估计方位时刻ti_est的多普勒频率值fD(ti_est)时，表达式为：
Figure FDA0002670640010000067
其中，Vs(ti_est)＝[Vsxi_est,Vsyi_est,Vszi_est]，Rs(ti_est)＝[Rsxi_est,Rsyi_est,Rszi_est]，RT＝[dx-ij,dy-ij,dz-ij]，R＝|Rs(ti_est)-RT|；
步骤S5-8中根据场景中心点的第一个方位时刻t1和最后一个方位时刻
Figure FDA0002670640010000068
的多普勒频率得到多普勒斜率时，第一个方位时刻t1和最后一个方位时刻
Figure FDA0002670640010000069
的多普勒频率分别记为fD_center(t1)、
Figure FDA00026706400100000610
多普勒斜率Kf表达式为：
Figure FDA00026706400100000611
步骤S5-9中由估计方位时刻ti_est和场景中心点的多普勒斜率计算多普勒中心方位时刻tI时，表达式为：
Figure FDA00026706400100000612
其中，
Figure FDA00026706400100000613
表示对x向下取整数；I表示SAR图像像元的行号；
步骤S5-10中由多普勒中心方位时刻的卫星位置Rs_I＝[RsxI,RsyI,RszI]和地面像素点T＝(dx-ij,dy-ij,dz-ij)计算二者之间斜距RJ时，表达式为：
Figure FDA0002670640010000071
Figure FDA0002670640010000072
其中，DT-R表示多普勒中心方位时刻tI的卫星位置Rs_I＝[RsxI,RsyI,RszI]和地面像素点T＝(dx-ij,dy-ij,dz-ij)之间距离，
Figure FDA0002670640010000073
表示对x向下取整数；J表示SAR图像像元的列号。
9.一种基于DSM数据的星载SAR图像快速间接定位装置，其特征在于，所述装置包括：
系统参数和图像参数获取模块，用于获取星载SAR图像的系统参数和图像参数；
轨道拟合模块，用于根据系统参数和图像参数中多个时刻的卫星位置和卫星速度，拟合卫星轨道曲线；
角点确定模块，用于根据拟合的所述卫星轨道曲线，由RD模型确定SAR图像定位在地面上的四个角点坐标，四个角点坐标为SAR图像对应的第一个方位时刻和最后一个方位时刻分别到最近斜距和最远斜距的WGS84坐标；
三维格网模块，用于利用所得的四个角点坐标构成一个四边形，获取这个四边形内的DSM数据构建以经度、纬度和高度三个维度构成的三维格网，根据DSM数据精度设定地面像素分辨率，每个三维格网点为一个地面像素点；
间接定位模块，用于对于四边形内每一个地面像素点，通过快速估计方位时刻的方式确定其对应的方位位置和距离位置；通过所述快速估计方位时刻的方式包括利用四边形内任一地面像素点的投影点和四个角点，以及三角形几何关系确定多普勒中心的估计方位时刻；通过估计方位时刻的多普勒值和场景中心点的多普勒斜率，确定实际多普勒中心方位时刻，得到该地面像素点对应的方位位置；根据多普勒中心方位时刻的卫星位置和地面像素点确定斜距，得到该地面像素点对应的距离位置。
10.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1至8中任一项所述的方法的步骤。
CN201910456759.3A 2019-05-29 2019-05-29 一种基于dsm数据的星载sar图像快速间接定位方法 Active CN110133653B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910456759.3A CN110133653B (zh) 2019-05-29 2019-05-29 一种基于dsm数据的星载sar图像快速间接定位方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910456759.3A CN110133653B (zh) 2019-05-29 2019-05-29 一种基于dsm数据的星载sar图像快速间接定位方法

Publications (2)

Publication Number Publication Date
CN110133653A CN110133653A (zh) 2019-08-16
CN110133653B true CN110133653B (zh) 2020-12-08

Family

ID=67582758

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910456759.3A Active CN110133653B (zh) 2019-05-29 2019-05-29 一种基于dsm数据的星载sar图像快速间接定位方法

Country Status (1)

Country Link
CN (1) CN110133653B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111289977A (zh) * 2020-03-02 2020-06-16 中国科学院电子学研究所 一种定位方法及装置
CN118068331B (zh) * 2024-04-22 2024-06-25 中国科学院空天信息创新研究院 星载合成孔径雷达数据流式分景方法、装置、设备及介质

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102243074A (zh) * 2010-05-13 2011-11-16 中国科学院遥感应用研究所 基于光线追踪技术的航空遥感成像几何变形仿真方法
CN102565797A (zh) * 2011-12-21 2012-07-11 北京航空航天大学 一种针对聚束模式星载sar图像的几何校正方法
US8242948B1 (en) * 2010-02-26 2012-08-14 Lockheed Martin Corporation High fidelity simulation of synthetic aperture radar
EP2446298B1 (de) * 2009-06-25 2014-12-03 Airbus Defence and Space GmbH Verfahren zur bestimmung der geographischen koordinaten von bildpunkten in sar bildern
CN105069843A (zh) * 2015-08-22 2015-11-18 浙江中测新图地理信息技术有限公司 一种面向城市三维建模的密集点云的快速提取方法
CN105510913A (zh) * 2015-11-11 2016-04-20 湖北工业大学 基于类光学像方改正的异源光学和sar遥感影像联合定位方法
CN106226768A (zh) * 2016-08-09 2016-12-14 北京空间飞行器总体设计部 超高分辨率敏捷sar卫星滑动聚束模式系统参数设计方法
CN107341778A (zh) * 2017-07-10 2017-11-10 国家测绘地理信息局卫星测绘应用中心 基于卫星控制点库和dem的sar影像正射纠正方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8842036B2 (en) * 2011-04-27 2014-09-23 Lockheed Martin Corporation Automated registration of synthetic aperture radar imagery with high resolution digital elevation models

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2446298B1 (de) * 2009-06-25 2014-12-03 Airbus Defence and Space GmbH Verfahren zur bestimmung der geographischen koordinaten von bildpunkten in sar bildern
US8242948B1 (en) * 2010-02-26 2012-08-14 Lockheed Martin Corporation High fidelity simulation of synthetic aperture radar
CN102243074A (zh) * 2010-05-13 2011-11-16 中国科学院遥感应用研究所 基于光线追踪技术的航空遥感成像几何变形仿真方法
CN102565797A (zh) * 2011-12-21 2012-07-11 北京航空航天大学 一种针对聚束模式星载sar图像的几何校正方法
CN105069843A (zh) * 2015-08-22 2015-11-18 浙江中测新图地理信息技术有限公司 一种面向城市三维建模的密集点云的快速提取方法
CN105510913A (zh) * 2015-11-11 2016-04-20 湖北工业大学 基于类光学像方改正的异源光学和sar遥感影像联合定位方法
CN106226768A (zh) * 2016-08-09 2016-12-14 北京空间飞行器总体设计部 超高分辨率敏捷sar卫星滑动聚束模式系统参数设计方法
CN107341778A (zh) * 2017-07-10 2017-11-10 国家测绘地理信息局卫星测绘应用中心 基于卫星控制点库和dem的sar影像正射纠正方法

Also Published As

Publication number Publication date
CN110133653A (zh) 2019-08-16

Similar Documents

Publication Publication Date Title
CA2705809C (en) Method and apparatus of taking aerial surveys
CN104268935A (zh) 一种基于特征的机载激光点云与影像数据融合系统及方法
CN111679288B (zh) 一种点云数据的空间分布度量方法
CN110133653B (zh) 一种基于dsm数据的星载sar图像快速间接定位方法
CN105444778B (zh) 一种基于成像几何反演的星敏感器在轨定姿误差获取方法
CN105004354A (zh) 大斜视角下无人机可见光和红外图像目标定位方法
CN109631863A (zh) 一种空地结合的潮间带一体化测绘方法
CN109239710A (zh) 雷达高程信息的获取方法及装置、计算机可读存储介质
CN111597496B (zh) 基于有理多项式系数单幅星载遥感影像目标高度计算方法
CN111487621B (zh) 一种基于雷达图像的海表流场反演方法及电子设备
CN111311659B (zh) 一种基于斜交平面镜三维成像的校准方法
CN113189551A (zh) 基于场景DEM的GB-InSAR重轨误差补偿方法
Molnár Making a georeferenced mosaic of historical map series using constrained polynomial fit
Spore et al. Collection, processing, and accuracy of mobile terrestrial lidar survey data in the coastal environment
CN109946682A (zh) 基于ICESat/GLAS的GF3数据基线估计方法
Zhang Photogrammetric processing of low altitude image sequences by unmanned airship
CN111681299B (zh) 基于InSAR解缠相位生成数字表面模型的方法及装置
CN115018973A (zh) 一种低空无人机点云建模精度的无靶标评估方法
CN113238202A (zh) 光子激光三维成像系统的坐标系点云计算方法及其应用
CN111611929A (zh) 基于LiDAR和InSAR技术的河道行洪风险点识别方法、装置、服务器及存储介质
CN117745779B (zh) 一种光学和sar共孔径一致性成像方法
Wu et al. Improvement of LiDAR data accuracy using 12-parameter affine transformation
Guastaferro et al. Rectification of spot 5 satellite imagery for marine geographic information systems
Meneghini et al. Use of Mercator cartographic representation for Landsat 8 imageries
AU2013260677B2 (en) Method and apparatus of taking aerial surveys

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
GR01 Patent grant