CN107341778A - 基于卫星控制点库和dem的sar影像正射纠正方法 - Google Patents
基于卫星控制点库和dem的sar影像正射纠正方法 Download PDFInfo
- Publication number
- CN107341778A CN107341778A CN201710554588.9A CN201710554588A CN107341778A CN 107341778 A CN107341778 A CN 107341778A CN 201710554588 A CN201710554588 A CN 201710554588A CN 107341778 A CN107341778 A CN 107341778A
- Authority
- CN
- China
- Prior art keywords
- dem
- sar
- coordinate
- data
- refdem
- 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 47
- 238000004088 simulation Methods 0.000 claims abstract description 64
- 238000012937 correction Methods 0.000 claims abstract description 63
- 241001269238 Data Species 0.000 claims abstract description 34
- 230000000149 penetrating effect Effects 0.000 claims abstract description 22
- 238000003384 imaging method Methods 0.000 claims description 73
- 239000011159 matrix material Substances 0.000 claims description 38
- 239000013598 vector Substances 0.000 claims description 37
- 238000013507 mapping Methods 0.000 claims description 26
- 238000012545 processing Methods 0.000 claims description 21
- 230000008569 process Effects 0.000 claims description 13
- 238000004364 calculation method Methods 0.000 claims description 8
- 238000005314 correlation function Methods 0.000 claims description 3
- 230000009467 reduction Effects 0.000 abstract description 9
- 238000005516 engineering process Methods 0.000 abstract description 4
- 230000000694 effects Effects 0.000 abstract description 2
- 238000006243 chemical reaction Methods 0.000 description 8
- 241000209094 Oryza Species 0.000 description 3
- 235000007164 Oryza sativa Nutrition 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 230000014509 gene expression Effects 0.000 description 3
- 238000005259 measurement Methods 0.000 description 3
- 230000003287 optical effect Effects 0.000 description 3
- 235000009566 rice Nutrition 0.000 description 3
- FNMKZDDKPDBYJM-UHFFFAOYSA-N 3-(1,3-benzodioxol-5-yl)-7-(3-methylbut-2-enoxy)chromen-4-one Chemical compound C1=C2OCOC2=CC(C2=COC=3C(C2=O)=CC=C(C=3)OCC=C(C)C)=C1 FNMKZDDKPDBYJM-UHFFFAOYSA-N 0.000 description 2
- 230000008901 benefit Effects 0.000 description 2
- 238000004422 calculation algorithm Methods 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 238000012821 model calculation Methods 0.000 description 2
- 238000003860 storage Methods 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 241000208340 Araliaceae Species 0.000 description 1
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 description 1
- 235000003140 Panax quinquefolius Nutrition 0.000 description 1
- 238000000205 computational method Methods 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000004069 differentiation Effects 0.000 description 1
- 235000008434 ginseng Nutrition 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 230000004807 localization Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000004321 preservation Methods 0.000 description 1
- 239000013643 reference control Substances 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000012876 topography Methods 0.000 description 1
Classifications
-
- G06T5/80—
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
- G06T2207/10044—Radar image
Landscapes
- Radar Systems Or Details Thereof (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种基于卫星控制点库和DEM的SAR影像正射纠正方法,其利用外部已有DEM数据和资源三号控制点数据,通过SAR影像模拟技术和影像精确匹配,得到资源三号控制点在SAR影像中的像点坐标,随后通过高精度的资源三号控制点参与SAR正射纠正,得到SAR正射纠正后的DOM影像。本发明可以自动地利用资源三号控制点库参与SAR影像正射纠正,减少人工干预带来的效率低下或人工选点带来的刺点误差;同时能有效用于SAR的正射纠正,降低获取控制点的成本,SAR影像纠正效果好。
Description
技术领域
本发明涉及雷达图像几何处理技术领域,特别涉及一种基于卫星控制点库和DEM的SAR影像正射纠正方法。
背景技术
资源三号测绘卫星于2012年1月9号成功发射,是我国第一颗民用三线阵立体测图卫星,主要用于1:50000立体测图及更大比例尺基础地理产品的生产和更新以及开展国土资源调查与监测。为满足国产资源三号测绘卫星影像数据处理、应急测绘保障等规模化应用,首次建立了覆盖全国范围的满足1:50000精度的影像控制点数据库,包含361025个高精度控制点数据,它们已全面应用于我国1:50000比例尺光学数字正射影像生产和1:25000比例尺的全国正射影像更新。由于光学影像与SAR影像的成像几何不同,散射机理不同,导致纹理和强度信息相差极大,因此使直接将资源三号卫星影像控制点转刺到SAR影像中,会带来数个像素的转刺误差。而如果直接采用严密成像模型直接定位,则会由于成像几何误差的影响带来更大的定位误差,直接应用于SAR影像正射纠正很难满足我国测绘比例尺的精度要求。
一般在无实测地面控制点数据时,SAR影像的正射纠正是利用外部参考DEM数据模拟SAR影像进行几何校正,其原理是:根据真实SAR影像的轨道状态矢量和成像参数以及外部参考DEM数据,通过SAR模拟技术模拟SAR影像的成像几何,得到模拟SAR影像。由于模拟SAR影像与真实SAR影像使用相同的轨道状态矢量和成像参数,它们之间的几何差异可以通过影像匹配确定,进而得到模拟SAR影像与真实SAR影像一一对应的像点映射关系。而模拟SAR影像与参考DEM之间同样存在一一映射的像点关系,可以在DEM上选取一部分点作为外部参考控制点数据,用于SAR影像的正射纠正处理。一般工程应用均使用的是WGS84参考椭球系统,WGS84椭球高与EGM96正高系统之间需要通过大地水准面差距进行转换。此外,国内学者对SRTM DEM的高程精度进行验证,利用三峡、华北、云贵地区总计206个GPS观测数据进行DEM高程精度验证,结果表明高程均方根误差为10.31m。资源三号卫星影像控制点数据的平面和高程精度优于厘米级,其高程精度要高于从SRTM DEM上选取的控制点数据,即使用资源三号卫星控制点数据可以进一步提高SAR影像正射纠正精度。现有技术中一般对SAR影像采用严密的成像模型进行定位处理和模型参数优化,采用的模型参数较多且计算复杂。
因此,需要提出一种能够克服上述缺点的使用资源三号卫星控制点数据和DEM辅助下的用于SAR影像正射纠正的方法。
发明内容
本发明提供了一种基于卫星控制点库和DEM的SAR影像正射纠正方法,其特征在于,所述方法具体包括以下步骤:
步骤S1,采集原始数据,所述原始数据优选的包括资源三号卫星影像控制点数据、外部参考DEM数据、原始SAR影像以及原始SAR影像成像参数;
步骤S2,对所述原始SAR影像进行多视处理,得到多视SAR影像;
步骤S3,基于所述外部参考DEM数据和多视SAR影像,计算外接矩形数字高程模型(RefDEM)数据、DEM模拟SAR影像、DEM模拟SAR影像和多视SAR影像雷达坐标空间一一映射的坐标关系;
步骤S4,根据步骤S3.3获取的DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系和步骤S3.4获取的DEM模拟SAR影像,计算得到与多视SAR影像雷达坐标空间一致的SAR模拟影像;
步骤S5,对所述SAR模拟影像和所述多视SAR影像进行灰度匹配,计算得到新的DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系;
步骤S6,基于所述新的DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系,计算资源三号卫星影像控制点数据在多视SAR影像中的像点坐标集;
步骤S7,将所述资源三号卫星影像控制点数据在多视SAR影像中的像点坐标集转换为所述资源三号卫星影像控制点数据在原始SAR影像中的坐标集;
步骤S8,利用资源三号卫星影像控制点数据在原始SAR影像中的坐标集,进行SAR影像正射纠正处理,计算得到距离向和方位向多项式的改正系数;
步骤S9,根据所述距离向和方位向多项式的改正系数,对所述原始SAR影像进行SAR正射纠正,得到SAR正射纠正后的DOM影像。
优选的,步骤2具体包括:根据下式对所述原始SAR影像进行多视处理,得到多视SAR影像:
其中,Pix(i,j)表示多视SAR影像中的像元(i,j)对应的像元值;nAzML表示原始SAR影像的方位向多视数,nRgML表示原始SAR影像的距离向多视数;Pix(n,m)表示原始SAR影像中的像元(n,m)对应的像元值,n表示原始SAR影像中的像元(n,m)的方位向坐标,取值范围为n=1,...,nAzML,m表示原始SAR影像中的像元(n,m)的距离向坐标,取值范围为m=1,...,nRgML。
优选的,所述步骤3优选的包括以下子步骤:
步骤S3.1,基于所述外部参考DEM数据和多视SAR影像,计算外接矩形数字高程模型(RefDEM)数据,其中,所述外接矩形数字高程模型(RefDEM)数据包括所述外接矩形数字高程模型(RfDEM)在外部参考DEM数据中的宽、高、起始点地理坐标、纬度方向分辨率、经度方向分辨率和高程值;
步骤S3.2,根据步骤S3.1获取的外接矩形数字高程模型(RefDEM)的宽nWidthDEM和高nHeightDEM,创建两个大小均为(nHeightDEM,nWidthDEM)的实数空矩阵,用于分别存储外接矩形数字高程模型(RefDEM)数据和DEM模拟SAR影像数据;创建一个大小为(nHeightDEM,nWidthDEM)的复数空矩阵,用于存储DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系数据;
步骤S3.3,基于所述外部参考DEM数据和多视SAR影像,根据SAR严密成像模型(R-D模型),计算外接矩形数字高程模型(RefDEM)数据与多视SAR影像雷达坐标空间一一映射的坐标关系;
步骤S3.4,根据SAR影像模拟公式,根据外接矩形数字高程模型(RefDEM)数据模拟SAR像元(i,j)处的影像值PixSim(i,j),最终所有模拟SAR像元值PixSim(i,j)的集合构成DEM模拟SAR影像。
优选的,所述步骤S3.1可以包括:
步骤S3.1.1,利用外部参考DEM(ExtDEM)数据和多视SAR影像的轨道状态、成像参数信息,根据SAR严密成像模型(R-D模型)计算得到多视SAR影像的四个角点坐标在外部参考DEM(ExtDEM)数据中的像素坐标;
步骤S3.1.2,根据步骤S3.1.1获取的多视SAR影像的四个角点在外部参考DEM(ExtDEM)数据中的像素坐标(iDEM,jDEM),计算外接矩形数字高程模型(RefDEM)的四个角点坐标在外部参考DEM(ExtDEM)数据中的坐标,以及外接矩形数字高程模型(RefDEM)在外部参考DEM数据中的宽nWidthDEM和高nHeightDEM;
步骤S3.1.3,根据步骤S3.1.2得到的外接矩形数字高程模型(RefDEM)的左上角Clt(i,j)角点坐标,结合外部参考DEM(ExtDEM)数据提供的起点地理坐标(Bext,0,Lext,0)、纬度方向的分辨率Bext,res、经度方向的分辨率Lext,res,根据下式计算外接矩形数字高程模型(RefDEM)的起始点地理坐标(B1,L1);
B1=Bext,0+Bext,res*i
L1=Lext,0+Lext,res*j
其中,(B1,L1)为外接矩形数字高程模型(RefDEM)的起始点地理坐标;(Bext,0,Lext,0)为外部参考DEM(ExtDEM)数据提供的起点地理坐标,Bext,res为外部参考DEM(ExtDEM)数据的纬度方向的分辨率,Lext,res为外部参考DEM(ExtDEM)数据的经度方向的分辨率;i为外接矩形数字高程模型(RefDEM)在外部参考DEM(ExtDEM)数据中的左上角的行方向坐标,j为外接矩形数字高程模型(RefDEM)在外部参考DEM(ExtDEM)数据中的左上角的列方向坐标;
步骤S3.1.4,根据步骤S3.1.3获取的外接矩形数字高程模型(RefDEM)的起始点地理坐标(B1,L1),计算外接矩形数字高程模型(RefDEM)的高程值。
优选的,所述步骤S3.1可以包括以下子步骤:
步骤S3.1.1.1,根据多视SAR影像的角点的方位向坐标iAZ,计算卫星方位向成像时间tT,并根据卫星方位向成像时间tT计算卫星位置矢量(Xs,Ys,Zs)和速度矢量
步骤S3.1.1.2,对多视SAR影像的中心点经纬度进行坐标转换,将其从大地坐标系转换到地心空间直角坐标系,作为待求点的地心坐标初值(XT 0,YT 0,ZT 0);
步骤S3.1.1.3,基于步骤S3.1.1.1求得的卫星位置矢量(Xs,Ys,Zs)和速度矢量和步骤S3.1.1.2求得的待求点的地心坐标初值(XT 0,YT 0,ZT 0),根据SAR严密成像模型,解算所述待求点的改正数;
步骤S3.1.1.4,根据下式,利用待求点的改正数更新待求点的地心坐标:
XT 1=XT 0+ΔXT,YT 1=YT 0+ΔYT,ZT 1=ZT 0+ΔZT;
其中,是更新后的待求点的地心坐标,是待求点的地心坐标初值,ΔXT,ΔYT,ΔZT是待求点的改正坐标;
步骤S3.1.1.5,利用更新后的待求点的地心坐标(XT 1,YT 1,ZT 1)进行迭代解算,重复S3.1.1.3-S3.1.1.4,直至改正数(ΔXT,ΔYT,ΔZT)的绝对值都小于设定阈值ε迭代结束,确定待求点的地心坐标(XT,YT,ZT);
步骤S3.1.1.6,对通过步骤S3.1.1.5获取的待求点的地心坐标(XT,YT,ZT)进行坐标转换,将其从地心空间直角坐标系转换到大地坐标系,得到并输出待求点的大地坐标(BT,LT,HT);
步骤S3.1.1.7,根据下述公式,计算得到多视SAR影像的角点坐标在外部参考DEM(ExtDEM)数据中的像素坐标(iDEM,jDEM):
iDEM=(Bext,0-Bi)/Bext,res
jDEM=(Lext,0-Li)/Lext,res
其中,(iDEM,jDEM)为多视SAR影像的角点坐标在外部参考DEM(ExtDEM)数据中的像素坐标;(Bext,0,Lext,0)为外部参考DEM(ExtDEM)数据的起点地理坐标;Bext,res为纬度方向分辨率,Lext,res为经度方向分辨率;Bi为多视SAR影像的角点坐标对应的大地坐标纬度值,Li为多视SAR影像角点坐标对应的大地坐标经度值。
优选的,步骤S3.1.1.3优选的包括:
所述SAR严密成像模型为:
其中,VxS·(XS-XP)+VyS·(YS-YP)+VzS·(ZS-ZP)=-fd·λ·(R0+ΔR·x)/2为多普勒方程,为距离方程,(XS,YS,ZS)为卫星位置矢量,(XP,YP,ZP)为角点的地心空间直角坐标,(VXS,VYS,VZS)为卫星速度矢量,x为角点的距离向坐标,R0为近地点斜距,ΔR为距离向分辨率,fd为多普勒中心频率,λ为雷达波长;
根据上述公式(2),构建观测方程,计算多视SAR影像的角点的误差方程系数,形成误差方程系数矩阵,再根据误差方程系数矩阵建立法方程:BTBX-BTL=0,根据下式解算待求点的改正数:
X=(BTB)-1BTL
其中,X=(ΔXT,ΔYT,ΔZT),ΔXT,ΔYT,ΔZT是待求点的改正数坐标;B为误差方程系数矩阵;L为常数项矩阵;;
其中,,在步骤S3.1.4中,根据下述公式计算外接矩形数字高程模型(RefDEM)的高程值:
BRefDEM=B1+Bref,res*m,m=0,1,…,nHeightDEM-1
LRefDEM=L1+Lref,res*n,n=0,1,…,n WidthDEM-1
iDEM=(Bext,0-BRefDEM)/Bext,res
jDEM=(Lext,0-LRefDEM)/Lext,res
HDEM=PixDEM(iDEM,jDEM)
其中,(BRefDEM,LRefDEM)为外接矩形数字高程模型(RefDEM)的某个点的经纬度坐标,BRefDEM为外接矩形数字高程模型(RefDEM)的某个点的纬度坐标,LRefDEM为外接矩形数字高程模型(RefDEM)的某个点的经度坐标;(B1,L1)为外接矩形数字高程模型(RefDEM)的起始点地理坐标,Bref,res为外接矩形数字高程模型(RefDEM)的纬度方向的分辨率,Lref,res为外接矩形数字高程模型(RefDEM)的经度方向的分辨率,nWidthDEM和nHeightDEM分别为外接矩形数字高程模型(RefDEM)的宽和高,m为高度方向的坐标索引,其中m=0,1,…,nHeightDEM-1;n为宽度方向的坐标索引,其中n=0,1,…,nWidthDEM-1;(Bext,0,Lext,0)为外部参考DEM(ExtDEM)数据的起点地理坐标,Bext,res为外部参考DEM(ExtDEM)在纬度方向的分辨率,Lext,res为外部参考DEM(ExtDEM)在经度方向的分辨率;(iDEM,jDEM)为外接矩形数字高程模型(RefDEM)的某个点的经纬度坐标在外部参考DEM(ExtDEM)数据中的行列坐标;PixDEM(iDEM,jDEM)表示行列坐标为(iDEM,jDEM)的点在外部参考DEM(ExtDEM)数据中的高程值,HDEM表示外接矩形数字高程模型(RefDEM)的高程值。
优选的,步骤S3.3优选的包括以下子步骤:
步骤S3.3.1,基于公式(2)中的多普勒方程迭代计算外接矩形数字高程模型(RefDEM)数据中的每个像素点的初始方位向坐标y;
S3.3.2,根据公式(2)的距离方程计算外接矩形数字高程模型(RefDEM)数据中的每个像素点的初始距离向坐标x;
S3.3.3,组合步骤S3.3.1和步骤S3.3.2得到的外接矩形数字高程模型(RefDEM)数据中的每个像素点的初始方位向坐标y和初始距离向坐标x,得到外接矩形数字高程模型(RefDEM)数据中的每个像素点在多视SAR影像雷达坐标空间中的像点坐标(x,y),所述像点坐标(x,y)为外接矩形数字高程模型(RefDEM)数据与多视SAR影像雷达坐标空间一一映射的坐标关系。
优选的,
步骤S4具体包括以下三个步骤:
步骤S4.1,获得DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系(即坐标映射数据1)的每一个像元pix(i,j)的复数像元值(iValue,jValue),如满足如下条件式,则保留该像元pix(i,j),参与下一步计算:
其中,ivalue为某一个像元pix(i,j)对应的复数像元值的高度方向坐标,jvalue为某一个像元pix(i,j)对应的复数像元值的宽度方向坐标,nMLHeight为多视SAR影像的高,nMLWidth为多视SAR影像的宽;
步骤S4.2,根据步骤S4.1获取的像元pix(i,j)和对应的复数像元值(iValue,jValue),在DEM模拟SAR影像中取(i,j)像元值填充到复数像元值(iValue,jValue)对应的SAR模拟影像坐标中;
步骤S4.3,对DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系(即坐标映射数据1)的每一个像元均采用步骤S4.1和步骤S4.2,得到SAR模拟影像;
步骤S5具体包括以下三个子步骤:
步骤S5.1,采用匹配方法,在SAR模拟影像上规则选取一定数量的点,对搜索窗口数为nWin范围内的像素进行匹配,得到窗口内每个像元的匹配测度ρ,取最大的匹配测度ρ,以及对应的方位向偏移量记为off_Az,距离向偏移量记为off_Rg。匹配方法采用常用的NCC匹配方法(归一化互相关系数),其公式如下:
其中,ri表示多视SAR影像(参考影像)窗口内的像素值;mi表示SAR模拟影像(待匹配影像)窗口内的像素值;表示多视SAR影像窗口内的像素均值;表示模拟SAR影像窗口内的像素均值;σr表示多视SAR影像窗口内的像素标准差;σm表示模拟SAR影像窗口内的像素标准差;n表示匹配窗口内的像素个数,n=nWin2,nWin为搜索窗口大小;
步骤S5.2,对步骤S5.1获取的方位向的偏移量off_Az、距离向的偏移量off_Rg,进行多项式拟合,得到方位向和距离向的偏移量拟合系数;多项式拟合公式为:
y(xAz)=a0+a1xAz
y(xRg)=b0+b1xRg
其中,y(xAz)、y(xRg)为方位向和距离向的偏移量矩阵,矩阵大小为k个元素;xAz、xRg为对应的元素索引,取值为xAz=1,2,...,k,xRg=1,2,...,k;a0、a1为方位向的多项式常数项和一次项系数;b0、b1为距离向的多项式常数项和一次项系数;
步骤S5.3,根据步骤S5.2获取的方位向和距离向的偏移量拟合系数,更新坐标映射数据1得到新的坐标映射数据2。更新的公式如下:
对坐标映射数据1的每一个像元(i,j),得到其复数像元值(iValue,1,jValue,1),根据方位向和距离向的偏移量,有
iValue,2=iValue,1+a0+a1·iValue,1
jValue,2=jValue,1+b0+b1·jValue,1
其中,iValue,2、jValue,2为更新后的坐标映射数据2的方位向和距离向像元值;iValue,1、jValue,1为坐标映射数据1的方位向和距离向像元值;a0、a1为方位向的多项式常数项和一次项系数;b0、b1为距离向的多项式常数项和一次项系数。
优选的,步骤S6具体包括以下子步骤:
步骤S6.1,计算ZY3数据在外接矩形数字高程模型(RefDEM)中的坐标(iDEM,jDEM),计算公式如下所示:
iDEM=(Bi-Bref,1)/Bref,res
jDEM=(Li-Lref,1)/Lref,res
其中,(Bref,1,Lref,1)为外接矩形数字高程模型(RefDEM)的起始点地理坐标;(Bi,Li)为资源三号卫星影像控制点数据ZY3的经纬度地理坐标;Bref,res为纬度方向外接矩形数字高程模型(RefDEM)分辨率,Lref,res为经度方向外接矩形数字高程模型(RefDEM)分辨率;(iDEM,jDEM)为ZY3数据在外接矩形RefDEM中的坐标;
步骤S6.2,计算(iDEM,jDEM)在SAR模拟影像中的像素坐标(iMSAR,jMSAR);
步骤S6.3,将所有ZY3数据重复步骤S6.1和S6.2,获得ZY3在SAR模拟影像中的像点坐标(iMSAR,jMSAR),并基于所述像点坐标(iMSAR,jMSAR)和ZY3的地理坐标和高程构成ZY3控制点数据在多视SAR影像中的坐标集(Bi,Li,Hi,iMSAR,jMSAR);其中,(Bi,Li,Hi)为ZY3的地理坐标和高程;(iMSAR,jMSAR)为ZY3在多视SAR影像中的像点坐标。
优选的,步骤S9具体包括以下子步骤:
步骤S9.1,根据原始SAR影像的四个角点像素坐标、原始SAR成像参数和输入的DOM数据的纬度方向分辨率、经度方向分辨率,确定待校正的DOM数据的宽和高,以及DOM起始点地理坐标;
步骤S9.2,根据步骤S9.1的地理参数,利用步骤S3.4构建的RD模型,逐像素对待校正的DOM影像进行计算,得到待校正的DOM像素在原始SAR影像中的像素坐标(x',y');
步骤S9.3,对步骤S9.2得到的像素坐标(x',y')进行改正,得到改正后的原始SAR像素坐标(x,y),获取原始SAR像素坐标(x,y)对应的像元值PixDOM(x,y),将所述像元值PixDOM(x,y)保存到步骤S9.2所述的待校正的DOM像素坐标pix(xDOM,yDOM)中,遍历整个待校正的DOM像素,获得DOM正射纠正影像。
相对于现有技术,本发明的主要优点在于:
1、通过本方法可以自动地利用资源三号控制点库参与SAR影像正射纠正,减少人工干预带来的效率低下或人工选点带来的刺点误差;
2、充分利用了通过光学手段获取的高精度资源三号控制点库,使其能有效用于SAR的正射纠正,降低获取控制点的成本;
3、通过外部DEM模拟SAR影像并与多视SAR影像进行匹配,确定了高精度资源三号控制点库在SAR影像上的坐标集,比直接在外部DEM上选取的点参与SAR影像纠正效果要好。
附图说明
图1为基于卫星控制点库和DEM的SAR正射纠正方法流程图。
图2(a)是外部参考DEM,图2(b)是多视SAR影像。
图3(a)为外接矩形RefDEM;图3(b)为坐标映射数据;图3(c)为DEM模拟SAR影像。
图4(a)为SAR模拟影像;图4(b)为多视SAR影像。
具体实施方式
下面参照附图对本发明进行更全面的描述,其中说明本发明的示例性实施例。
本发明的主要目的是利用外部已有DEM数据和资源三号控制点数据,通过SAR影像模拟技术和影像精确匹配,得到资源三号控制点在SAR影像中的像点坐标,随后通过高精度的资源三号控制点(此时包含像点坐标和经纬度高程坐标)参与SAR正射纠正,得到SAR正射纠正后的DOM影像。
如图1所示,本发明提出的基于卫星控制点库和DEM(Digital Elevation Model,DEM;数字高程模型)的SAR正射纠正方法,利用资源三号卫星影像控制点数据和外部参考DEM数据,首先对原始SAR影像进行多视处理得到多视SAR影像,其次通过SAR模拟技术,模拟一副SAR影像,运用匹配策略对SAR模拟影像和多视SAR影像进行匹配,确定它们之间的几何偏移关系,随后通过几何偏移关系计算得到高精度资源三号卫星影像控制点数据在多视SAR影像中的像点坐标,进而得到资源三号卫星影像控制点数据在原始SAR影像中的像点坐标,最后利用资源三号卫星影像控制点的高程数据替换外部参考DEM对应的高程数据,得到高精度的资源三号控制点(此时包含像点坐标和经纬度高程坐标),参与SAR影像正射纠正,其可实现资源三号卫星影像控制点数据应用于SAR影像正射纠正处理。该方法具体包括以下步骤:
步骤S1,采集原始数据,所述原始数据包括资源三号卫星影像控制点数据、外部参考DEM数据、原始SAR影像以及原始SAR影像成像参数。优选的,外部参考DEM数据可以是SRTM(Shuttle Radar Topography Mission,SRTM;美国航天飞机雷达测图任务)DEM数据。
具体数据如下:1)资源三号卫星影像控制点数据采用高精度的GPS设备进行外业测量,从资源三号卫星影像控制点数据库中查找并选取实验区控制点数据,如表1所示的资源三号卫星影像控制点数据,以ZY3表示,控制点坐标包括:第一列表示控制点的名称,第二列表示纬度,第三列表示经度,第四列表示大地高程;2)外部参考DEM数字高程模型数据,以ExtDEM表示,它的地理范围是完全包括SAR影像的地理范围;3)原始SAR影像,为了描述方便,将步骤1获取的SAR影像称为“原始SAR影像”;4)原始SAR影像对应的成像参数如表2所示。
控制点名称 | 纬度/度 | 经度/度 | 大地高程/米 |
F03 | 39.80890 | 116.09411 | 75.428 |
F18 | 39.74416 | 116.01044 | 50.419 |
F05 | 39.79243 | 116.16860 | 38.754 |
C24 | 39.88796 | 116.01093 | 220.683 |
C22 | 39.89176 | 116.02806 | 230.632 |
W15 | 40.12384 | 116.22403 | 30.164 |
C02 | 40.12662 | 116.22181 | 33.304 |
W17 | 40.06218 | 116.28779 | 36.213 |
表1 资源三号卫星影像控制点数据ZY3
表2 原始SAR影像对应的成像参数
由于资源三号卫星影像控制点数据精度更高,因此本发明用资源三号卫星影像控制点数据替换具有相同地理坐标的外部参考DEM高程数据,参与SAR正射纠正处理,能得到更好的SAR影像正射纠正精度,表3为资源三号卫星影像控制点数据与SRTM DEM高程精度对比结果。
控制点名称 | 纬度/度 | 经度/度 | 实测高程/米 | SRTM高程/米 | 误差值 |
F03 | 39.80890 | 116.09411 | 75.428 | 76.727 | -1.299 |
F18 | 39.74416 | 116.01044 | 50.419 | 53.834 | -3.415 |
F05 | 39.79243 | 116.16860 | 38.754 | 37.279 | 1.475 |
C24 | 39.88796 | 116.01093 | 220.683 | 234.741 | -14.058 |
C22 | 39.89176 | 116.02806 | 230.632 | 234.092 | -3.460 |
W15 | 40.12384 | 116.22403 | 30.164 | 37.105 | -6.941 |
C02 | 40.12662 | 116.22181 | 33.304 | 34.631 | -1.327 |
W17 | 40.06218 | 116.28779 | 36.213 | 38.467 | -2.254 |
表3 资源三号卫星影像控制点数据与SRTM DEM高程精度比较
步骤S2,对所述原始SAR影像进行多视处理,得到多视SAR影像。
具体的,在步骤S2中,利用步骤S1获得的原始SAR影像和原始SAR影像成像参数,对原始SAR影像进行多视处理,得到多视SAR影像。
考虑到原始SAR影像受噪声影响比较大,为了增加SAR模拟影像与SAR影像匹配的成功率,首先对原始SAR影像进行多视处理,达到降低噪声对匹配结果的影响。将原始SAR影像的方位向(行方向)多视数设为nAzML,距离向(列方向)多视数设为nRgML。假设一幅原始SAR影像的宽为nWidth,高为nHeight,那么进行多视处理之后,其多视SAR影像的宽nMLWidth为nWidth/nRgML,多视SAR影像的高nMLHeight为nHeight/nAzML。根据下式对所述原始SAR影像进行多视处理,得到多视SAR影像:
其中,Pix(i,j)表示多视SAR影像中的像元(i,j)对应的像元值;nAzML表示原始SAR影像的方位向多视数,nRgML表示原始SAR影像的距离向多视数;Pix(n,m)表示原始SAR影像中的像元(n,m)对应的像元值,n表示原始SAR影像中的像元(n,m)的方位向坐标,取值范围为n=1,...,nAzML,m表示原始SAR影像中的像元(n,m)的距离向坐标,取值范围为m=1,...,nRgML。对所有像元均采用公式(1)进行计算,得到多视SAR影像。
优选的,方位向多视数nAzML=4,距离向多视数nRgML=4;多视SAR影像的宽nMLWidth=3461,高nMLHeight=7176。
步骤S3,基于所述外部参考DEM数据和多视SAR影像,计算外接矩形数字高程模型(RefDEM)数据、DEM模拟SAR影像、DEM模拟SAR影像和多视SAR影像雷达坐标空间一一映射的坐标关系。
具体的,利用步骤S1提供的外部参考DEM数据和步骤S2提供的多视SAR影像的轨道状态和成像参数信息,计算得到外接矩形数字高程模型(RefDEM)数据,通过SAR模拟技术得到DEM模拟SAR影像,以及计算得到DEM模拟SAR影像和多视SAR影像之间一一映射的坐标关系。
SAR模拟技术,首先包括几何模拟,它是以数字高程模型(DEM)、卫星状态矢量参数和成像参数为基本输入信息,基于距离-多普勒方程(如公式2所示)建立几何构像模型,模拟出影像的几何位置;其次,后向散射模拟,本发明基于半经验后向散射模型进行SAR影像的灰度模拟。完成模拟SAR影像之后,采用影像匹配技术,自动完成SAR模拟影像与多视SAR影像上同名目标点的寻找,是SAR正射纠正处理中常用的一种方法。
RefDEM(即,上文所述的外接矩形数字高程模型)数据是模拟SAR影像对应的外接矩形数字高程模型数据;DEM模拟SAR影像是基于外接矩形数字高程模型(RefDEM)数据模拟计算得到的模拟SAR影像,其大小和影像方向与外接矩形数字高程模型(RefDEM)数据相一致。
步骤S3包括:首先根据外部参考DEM数据和多视SAR影像的轨道状态、成像参数信息,结合多视SAR影像的四个角点像元坐标,计算外接矩形数字高程模型(RefDEM)在外部参考DEM数据中的宽nWidthDEM和高nHeightDEM,以及相应外接矩形数字高程模型(RefDEM)数据的起始点地理坐标(B1,L1);其次,创建两个大小为(nHeightDEM,nWidthDEM)的实数空矩阵,用于存储外接矩形数字高程模型(RefDEM)数据和DEM模拟SAR影像数据,创建一个大小为(nHeightDEM,nWidthDEM)的复数空矩阵,用于存储DEM模拟SAR影像和多视SAR影像雷达坐标空间一一映射的坐标关系数据;随后,根据外接矩形数字高程模型(RefDEM)数据的起始点地理坐标(B1,L1),在外部参考DEM数据中取高程值,数据值存储于外接矩形数字高程模型(RefDEM)数据中;根据SAR严密成像模型(R-D模型),结合外接矩形数字高程模型(RefDEM)数据的起始点地理坐标(B1,L1)和相应高程数据,计算其在多视SAR影像雷达坐标空间中的像点坐标,存储于DEM模拟SAR影像和多视SAR影像雷达坐标空间一一映射的坐标关系数据中;随后,根据SAR影像模拟公式(即半经验后向散射模型),模拟SAR影像灰度值,存储于DEM模拟SAR数据中;最后输出三个数据文件,参与后续影像匹配、资源三号控制点坐标获取和SAR正射纠正处理。
具体来说,步骤S3具体包括以下子步骤:
步骤S3.1,基于所述外部参考DEM数据和多视SAR影像,计算外接矩形数字高程模型(RefDEM)数据,其中,所述外接矩形数字高程模型(RefDEM)数据包括外接矩形数字高程模型(RefDEM)在外部参考DEM数据中的宽、高、起始点地理坐标、纬度方向分辨率、经度方向分辨率和高程值。
具体的,在步骤S3.1中,利用外部参考DEM(ExtDEM)数据和多视SAR影像的轨道状态矢量、成像参数信息,计算外接矩形数字高程模型(RefDEM)在外部参考DEM数据中的宽nWidthDEM和高nHeightDEM,以及外接矩形数字高程模型(RefDEM)的起始点地理坐标(B1,L1)。
如图2(a)所示,其中,21表示外部参考DEM数据,22表示多视SAR影像的实际范围,23表示计算得到的外接矩形数字高程模型数据。图2(a)是图2(b)所使用到的外部参考DEM数据,图2(a)的范围要完全包含图2(b)的范围,图2(b)对应的是图2(a)中的“多视SAR影像的实际范围22”所示的框线。
外部参考DEM数据21的范围很大,所以要根据多视SAR影像参数文件(SAR影像的实际范围如图2(a)的22所示),进行裁剪,只保留与SAR影像范围一致的DEM数据,即计算得到外接矩形数字高程模型(RefDEM)(如图2(a)的23所示)。
如图2(a)所示,要确定外接矩形DEM数据在外部参考DEM数据中的范围,首先是根据图2(a)中多视SAR影像的实际范围22(“斜着的白色”矩形)的4个角点坐标,反算这4个角点坐标在外部参考DEM中的4个新坐标,进而根据这4个新坐标确定外接矩形数字高程模型23(“正着的白色”矩形)的四个新角点坐标,最后根据四个新角点坐标得到宽nWidthDEM和高nHeightDEM以及外接矩形数字高程模型23的左上角起始点地理坐标(B1,L1)。
步骤S3.1具体包括以下子步骤:
步骤S3.1.1,利用外部参考DEM(ExtDEM)数据和多视SAR影像的轨道状态矢量、成像参数信息,根据SAR严密成像模型(R-D模型)计算得到多视SAR影像的四个角点坐标在外部参考DEM(ExtDEM)数据中的像素坐标;
具体的,利用外部参考DEM(ExtDEM)数据和多视SAR影像的轨道状态矢量、成像参数信息,根据SAR严密成像模型(R-D模型)进行直接地理定位得到多视SAR影像数据的每个角点的地理坐标,所述SAR严密成像模型公式如下式(2)所示,计算多视SAR影像的四个角点坐标在外部参考DEM(ExtDEM)数据中的坐标,得到多视SAR影像数据的四个角点坐标在外部参考DEM(ExtDEM)数据中的四个角点坐标为左上角坐标Clt'(i,j)、右上角坐标Crt'(i,j)、左下角坐标Cld'(i,j)、右下角坐标Crd'(i,j);
其中,VxS·(XS-XP)+VyS·(YS-YP)+VzS·(ZS-ZP)=-fd·λ·(R0+ΔR·x)/2为多普勒方程,为距离方程,(XS,YS,ZS)为卫星位置矢量,(XP,YP,ZP)为角点的地心空间直角坐标,(VXS,VYS,VZS)为卫星速度矢量,x为角点的距离向坐标,R0为近地点斜距,ΔR为距离向分辨率,fd为多普勒中心频率,λ为雷达波长。在求解的过程中,均是在WGS84椭球系下进行。
步骤S3.1.1具体包括以下子步骤:
步骤S3.1.1.1,根据多视SAR影像的角点(所述角点是指多视SAR影像的左上角点、右上角点、左下角点、右下角点这四个角点中的一个,需要说明的是,这四个角点坐标的计算方法相同)的方位向坐标iAZ,计算卫星方位向成像时间tT,并根据卫星方位向成像时间tT计算卫星位置矢量(Xs,Ys,Zs)和速度矢量
具体的,根据下式计算卫星方位向成像时间tT:tT=t0+tAzi·iAz,其中,iAz表示方位向坐标,t0表示方位向起始时间,tAzi表示方位向时间间隔。
其中,方位向起始时间为t0,方位向时间间隔tAzi(即,多视SAR影像的成像参数信息),均由表2成像参数信息提供。优选的,令多视SAR影像的左上角的方位向和距离向坐标为(0,0),右上角为(0,nMLWidth),左下角为(nMLHeight,0),右下角(nMLHeight,nMLWidth)。以多视SAR影像的左上角点(0,0)为例进行说明:多视SAR影像的四个角点方位向和距离向坐标(iAz,jRg)分别为左上角坐标(0,0)、右上角坐标(0,3461)、左下角坐标(7176,0)、右下角坐标(7176,3461)。
根据卫星方位向成像时间tT计算卫星位置矢量(Xs,Ys,Zs)具体包括:
拟合卫星位置矢量与所述卫星方位向成像时间tT之间的二阶多项式,所述二阶多项式如下公式所示:
其中,(XS,YS,ZS)为待拟合的卫星位置矢量,为位置拟合二阶多项式系数矩阵,tT为卫星方位向成像时间;
以X方向为例,使用表2提供的19个状态位置矢量(即,上文提到的多视SAR影像的轨道状态信息)可构建19个方程,即
X=(a0 a1 a2)T
其中,X=(-2628711.7465 -2629649.8296 -2630583.9876 -2631514.2213 -2632440.5295 -2633362.9117 -2634281.3685 -2635195.8984 -2636106.5024 -2637013.1796 -2637915.9289 -2638814.7513 -2639709.6459 -2640600.6118 -2641487.6501 -2642370.7589 -2643249.9397 -2644125.1910 -2644996.5124)';T=(t0t1 … t12)′=(80787.910594 80788.910594 80789.910594 80790.910594 80791.91059480792.910594 80793.910594 80794.910594 80795.910594 80796.910594 80797.91059480798.910594 80799.910594 80800.910594 80801.910594 80802.910594 80803.91059480804.910594 80805.910594),状态矢量的时间间隔为1s。解算可获得多项式系数,
(ax0 ax1 ax2)=(-2627769.7393 -943.9692 1.9619)
以上计算得到了多项式系数,进而根据多项式系数和卫星方位向成像时间计算任意位置的卫星位置矢量。
根据卫星方位向成像时间tT计算卫星速度矢量具体包括:
拟合卫星速度矢量与所述卫星方位向成像时间tT之间的二阶多项式,所述二阶多项式如公式(4)所示:
其中,为待拟合的卫星速度矢量,tT为卫星方位向成像时间,为速度拟合二阶多项式系数矩阵。
步骤S3.1.1.2,对多视SAR影像的中心点经纬度进行坐标转换,将其从大地坐标系转换到地心空间直角坐标系,作为待求点的地心坐标初值(XT 0,YT 0,ZT 0),所述待求点是指多视SAR影像的左上角点、右上角点、左下角点、右下角点四个角点中的一个(需要说明的是,这四个角点坐标的计算方法相同);
步骤S3.1.1.3,基于步骤S3.1.1.1求得的卫星位置矢量(Xs,Ys,Zs)和速度矢量和步骤S3.1.1.2求得的待求点的地心坐标初值(XT 0,YT 0,ZT 0),根据SAR严密成像模型,解算所述待求点的改正数;
具体的,如上所述,所述SAR严密成像模型为:
其中,VxS·(XS-XP)+VyS·(YS-YP)+VzS·(ZS-ZP)=-fd·λ·(R0+ΔR·x)/2为多普勒方程,为距离方程,(XS,YS,ZS)为卫星位置矢量,(XP,YP,ZP)为角点的地心空间直角坐标,为卫星速度矢量,x为角点的距离向坐标,R0为近地点斜距,ΔR为距离向分辨率,fd为多普勒中心频率,λ为雷达波长。
根据上述公式(2),构建观测方程,计算多视SAR影像的角点的误差方程系数,形成误差方程系数矩阵,再根据误差方程系数矩阵建立法方程(该过程为本领域公知的技术,在此不再赘述):BTBX-BTL=0,根据下式解算待求点的改正数:
X=(BTB)-1BTL
其中,X=(ΔXT,ΔYT,ΔZT),ΔXT,ΔYT,ΔZT是待求点的改正数坐标;B为误差方程系数矩阵(在该步骤的计算中,上述公式(2)中所涉及的卫星位置矢量和速度矢量由步骤S3.1.1.1计算得到);L为常数项矩阵,由步骤S3.1.1.2计算得到。
步骤S3.1.1.4,根据下式,利用待求点的改正数更新待求点的地心坐标:
XT 1=XT 0+ΔXT,YT 1=YT 0+ΔYT,ZT 1=ZT 0+ΔZT;
其中,是更新后的待求点的地心坐标,是待求点的地心坐标初值,ΔXT,ΔYT,ΔZT是待求点的改正坐标。
步骤S3.1.1.5,利用更新后的待求点的地心坐标(XT 1,YT 1,ZT 1)进行迭代解算,重复S3.1.1.3-S3.1.1.4,直至改正数(ΔXT,ΔYT,ΔZT)的绝对值都小于设定阈值ε迭代结束,确定待求点的地心坐标(XT,YT,ZT);根据本发明的优选实施例,阈值ε为10-8。
步骤S3.1.1.6,对通过步骤S3.1.1.5获取的待求点的地心坐标(XT,YT,ZT)进行坐标转换,将其从地心空间直角坐标系转换到大地坐标系,得到并输出待求点的大地坐标(BT,LT,HT);
该大地坐标(BT,LT,HT)是多视SAR影像的直接定位结果,作为多视SAR影像待求点的大地坐标(即经纬度坐标)。
步骤S3.1.1.7,根据下述公式(5),计算得到多视SAR影像的角点坐标在外部参考DEM(ExtDEM)数据中的像素坐标(iDEM,jDEM);
具体的,根据外部参考DEM(ExtDEM)数据提供的起点地理坐标(Bext,0,Lext,0)、纬度方向的分辨率Bext,res、经度方向的分辨率Lext,res,再结合步骤S3.1.1.6输出的多视SAR影像的角点的大地坐标,利用下述公式(5)进行计算,得到多视SAR影像角点坐标在外部参考DEM(ExtDEM)数据中的像素坐标(iDEM,jDEM);
其中,(iDEM,jDEM)为多视SAR影像的角点坐标在外部参考DEM(ExtDEM)数据中的像素坐标;(Bext,0,Lext,0)为外部参考DEM(ExtDEM)数据的起点地理坐标;Bext,res为外部参考DEM(ExtDEM)数据的纬度方向分辨率,Lext,res为外部参考DEM(ExtDEM)数据的经度方向分辨率;Bi为多视SAR影像的角点坐标对应的大地坐标纬度值,Li为多视SAR影像角点坐标对应的大地坐标经度值,由步骤S3.1.1.6获得。
需要说明的是,计算多视SAR影像的左上角点、右上角点、左下角点、右下角点这四个角点的大地坐标的过程相同,所以通过步骤S3.1.1.1-S3.1.1.7可以计算得到多视SAR影像的左上角点、右上角点、左下角点、右下角点这四个角点的大地坐标。
根据本发明的具体实施例,外部参考DEM(ExtDEM)数据的起点地理坐标(Bext,0,Lext,0)=(41.00042,114.999583)、纬度方向分辨率为Bext,res=-8.3333333e-5、经度方向分辨率为Lext,res=8.3333333e-5。由此计算得到的多视SAR影像的四个角点坐标在外部参考DEM(ExtDEM)数据中的坐标分别为左上角Clt'(i,j)=(9315,13865)、右上角Crt'(i,j)=(12374,19341)、左下角Cld'(i,j)=(17310,9605)、右下角Crd'(i,j)=(19620,16874)。
步骤S3.1.2,根据步骤S3.1.1获取的多视SAR影像的四个角点在外部参考DEM(ExtDEM)数据中的像素坐标(iDEM,jDEM),计算外接矩形数字高程模型
(RefDEM)的四个角点(即,外接矩形数字高程模型的左上角点、右上角点、左下角点、右下角点)坐标在外部参考DEM(ExtDEM)数据中的坐标,以及外接矩形数字高程模型(RefDEM)在外部参考DEM数据中的宽nWidthDEM和高nHeightDEM;
在步骤S3.1.1中得到多视SAR影像四个角点在外部参考DEM(ExtDEM)数据中的像素坐标(iDEM,jDEM),即图2(a)中的附图标记22所示的框图的四个角点在外部参考DEM(ExtDEM)数据中的(iDEM,jDEM)。可以看到附图标记22所示的框不是一个规则的矩形,为了计算方便,需要寻找(iDEM,jDEM)中的行和列的最大最小值,最终确定图2中的附图标记23所示的框图的像素坐标和该框图的宽、高。
要确定外接矩形数字高程模型(RefDEM)在外部参考DEM(ExtDEM)数据中的边界范围,必须对多视SAR影像的四个角点在外部参考DEM(ExtDEM)数据中的像素坐标(iDEM,jDEM)进行判断,寻找它们中的行方向的最大、最小值,列方向的最大、最小值,而行方向的最大值与最小值之差即为宽nWidthDEM,列方向的最大值与最小值之差即为高nHeightDEM,进而得到外接矩形数字高程模型(RefDEM)在外部参考DEM数据中的宽nWidthDEM和高nHeightDEM;另外,在对多视SAR影像的四个角点在外部参考DEM(ExtDEM)数据中的像素坐标(iDEM,jDEM)进行判断,寻找它们中的行方向的最大、最小值,列方向的最大、最小值后,可获得外接矩形数字高程模型(RefDEM)的四个角点坐标在外部参考DEM(ExtDEM)数据中分别为左上角坐标Clt(i,j)(即行最小,列最小)、右上角坐标Crt(i,j)(即行最小,列最大)、左下角坐标Cld(i,j)(即行最大,列最小)、右下角坐标Crd(i,j)(即行最大,列最大);
举例来说,根据步骤S3.1.1.7提供的多视SAR影像的四个角点坐标在外部参考DEM(ExtDEM)数据中的坐标分别为左上角Clt'(i,j)=(9315,13865)、右上角Crt'(i,j)=(12374,19341)、左下角Cld'(i,j)=(17310,9605)、右下角Crd'(i,j)=(19620,16874),它们中的行方向最大、最小值分别为19341和9605,它们中的列方向最大、最小值分别为19620和9315,因此外接矩形数字高程模型(RefDEM)的四个角点坐标在外部参考DEM(ExtDEM)数据中分别为左上角Clt(i,j)=(9315,9605)(即行最小,列最小)、右上角Crt(i,j)=(9315,19341)(即行最小,列最大)、左下角Cld(i,j)=(19620,9605)(即行最大,列最小)、右下角Crd(i,j)=(19620,19341)(即行最大,列最大)。
步骤S3.1.3,根据步骤S3.1.2得到的外接矩形数字高程模型(RefDEM)的左上角Clt(i,j)角点坐标,结合外部参考DEM(ExtDEM)数据提供的起点地理坐标(Bext,0,Lext,0)(SAR处理领域,一般以左上角为图像的起始点)、纬度方向的分辨率Bext,res、经度方向的分辨率Lext,res,由下式(6)计算外接矩形数字高程模型(RefDEM)的起始点地理坐标(B1,L1);
其中,(B1,L1)为外接矩形数字高程模型(RefDEM)的起始点地理坐标;
(Bext,0,Lext,0)为外部参考DEM(ExtDEM)数据提供的起点地理坐标,Bext,res为外部参考DEM(ExtDEM)数据的纬度方向的分辨率,Lext,res为外部参考DEM(ExtDEM)数据的经度方向的分辨率;i为外接矩形数字高程模型(RefDEM)在外部参考DEM(ExtDEM)数据中的左上角的行方向坐标,j为外接矩形数字高程模型(RefDEM)在外部参考DEM(ExtDEM)数据中的左上角的列方向坐标(该行、列坐标通过步骤S3.1.2求得)。
此外,外接矩形数字高程模型(RefDEM)的纬度方向的分辨率(Bref,res)和经度方向的分辨率(Lref,res)分别与外部参考DEM(ExtDEM)数据的纬度方向分辨率和经度方向的分辨率相同,即:
Bref,res=Bext,res,
Lref,res=Lext,res。
根据本发明的具体实施例,外接矩形数字高程模型(RefDEM)的宽和高分别为nWidthDEM=9736和nHeightDEM=10305。外接矩形数字高程模型(RefDEM)在外部参考ExtDEM数据中的范围,如图2中的23所示。相应的外接矩形数字高程模型(RefDEM)的起始点地理坐标为(B1,L1)=(40.224,115.799)。
步骤S3.1.4,根据步骤S3.1.3获取的外接矩形数字高程模型(RefDEM)的起始点地理坐标(B1,L1),计算外接矩形数字高程模型(RefDEM)的高程值;
具体的,根据步骤S3.1.3获取的外接矩形数字高程模型(RefDEM)的起始点地理坐标(B1,L1)、外接矩形数字高程模型(RefDEM)的纬度方向的分辨率Bref,res、外接矩形数字高程模型(RefDEM)的经度方向的分辨率Lref,res,在外部参考DEM(ExtDEM)数据中取高程值,并将取得的高程值存储于外接矩形数字高程模型(RefDEM)中。
具体的,根据下述公式(7)计算外接矩形数字高程模型(RefDEM)的高程值:
其中,(BRefDEM,LRefDEM)为外接矩形数字高程模型(RefDEM)的某个点的经纬度坐标,BRefDEM为外接矩形数字高程模型(RefDEM)的某个点的纬度坐标,LRefDEM为外接矩形数字高程模型(RefDEM)的某个点的经度坐标;(B1,L1)为外接矩形数字高程模型(RefDEM)的起始点地理坐标,Bref,res为外接矩形数字高程模型(RefDEM)的纬度方向的分辨率,Lref,res为外接矩形数字高程模型(RefDEM)的经度方向的分辨率,nWidthDEM和nHeightDEM分别为外接矩形数字高程模型(RefDEM)的宽和高,m为高度方向的坐标索引,其中m=0,1,…,nHeightDEM-1;n为宽度方向的坐标索引,其中n=0,1,…,nWidthDEM-1;(Bext,0,Lext,0)为外部参考DEM(ExtDEM)数据的起点地理坐标,Bext,res为外部参考DEM(ExtDEM)在纬度方向的分辨率,Lext,res为外部参考DEM(ExtDEM)在经度方向的分辨率;(iDEM,jDEM)为外接矩形数字高程模型(RefDEM)的某个点的经纬度坐标在外部参考DEM(ExtDEM)数据中的行列坐标,即,iDEM为计算得到的外接矩形数字高程模型(RefDEM)的像素坐标在外部参考DEM(ExtDEM)数据中的行坐标,jDEM为计算得到的外接矩形数字高程模型(RefDEM)的像素坐标在外部参考DEM(ExtDEM)数据中的列坐标;PixDEM(iDEM,jDEM)表示行列坐标为(iDEM,jDEM)的点在外部参考DEM(ExtDEM)数据中的高程值,HDEM表示外接矩形数字高程模型(RefDEM)的高程值。
根据本发明的具体实施例,外部参考DEM(ExtDEM)数据的起点地理坐标(Bext,0,Lext,0)=(41.00042,114.999583)、纬度方向分辨率为Bext,res=-8.3333333e-5、经度方向、分辨率为Lext,res=8.3333333e-5。
步骤S3.2,根据步骤S3.1获取的外接矩形数字高程模型(RefDEM)的宽nWidthDEM和高nHeightDEM,创建两个大小均为(nHeightDEM,nWidthDEM)的实数空矩阵,用于分别存储外接矩形数字高程模型(RefDEM)数据和DEM模拟SAR影像数据;创建一个大小为(nHeightDEM,nWidthDEM)的复数空矩阵,用于存储DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系数据。
在图3a,3b,3c中的数据对应的宽和高是根据步骤S3.1获取的;为了后续计算的方便,在步骤S3.2中先创建三个空的文件,创建两个大小均为(nHeightDEM,nWidthDEM)的实数空矩阵,分别存储外接矩形数字高程模型(RefDEM)数据(最终结果为图3a,对应的步骤为S3.1)和DEM模拟SAR影像数据(最终结果为图3c,对应的步骤为S3.4);
创建一个大小为(nHeightDEM,nWidthDEM)的复数空矩阵,用于存储DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系(即坐标映射数据1)(最终结果为图3b,对应的步骤为S3.3)。
此外还可以看到,图3c的结果为DEM模拟SAR影像,它的方向与图3a中的RefDEM的方向是一致的;而为了和原始SAR影像进行匹配,必须将图3c转换为与SAR影像一致的方向上,这个处理过程在步骤S4进行。
步骤S3.3,基于所述外部参考DEM数据和多视SAR影像,根据SAR严密成像模型(R-D模型),计算外接矩形数字高程模型(RefDEM)数据与多视SAR影像雷达坐标空间一一映射的坐标关系;
具体的,根据SAR严密成像模型(R-D模型),结合由步骤S3.1获得的外接矩形数字高程模型(RefDEM)数据的起始点地理坐标和高程,以及多视SAR影像的成像参数中的卫星状态矢量参数,计算外接矩形数字高程模型(RefDEM)数据中的每个像素点在多视SAR影像雷达坐标空间中的像点坐标(x,y),所述像点坐标(x,y)为外接矩形数字高程模型(RefDEM)数据与多视SAR影像雷达坐标空间一一映射的坐标关系,通过这个映射关系,即可在两个数据(不同坐标系)中相互找到对应关系。所述SAR严密成像模型如公式(2)所示。
DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系数据如图3(b)所示。因为,DEM模拟SAR影像是根据外接矩形数字高程模型(RefDEM)数据和SAR模拟技术得到的。这个坐标映射数据起着桥梁衔接的作用,通过坐标映射数据,可以实现DEM模拟SAR影像与多视SAR影像雷达坐标空间任何坐标的查找。
在得到DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系数据(如图3(b))之后,即可根据坐标映射数据得到SAR模拟影像(如图4(a)所示,对应步骤S4)。
步骤S3.3具体包括以下子步骤:
S3.3.1,基于公式(2)中的多普勒方程迭代计算外接矩形数字高程模型(RefDEM)数据中的每个像素点的初始方位向坐标y;
步骤S3.3.1具体包括以下步骤:
S3.3.1.1,设定外接矩形数字高程模型(RefDEM)数据中的每个像素点的像素坐标初值为多视SAR影像的中心点(icenter,jcenter)(其中jcenter=nMLWidth/2,icenter=nMLHeight/2,nMLWidth和nMLHeight分别为多视SAR影像的宽和高),根据多视SAR影像中心点的方位向坐标icenter计算卫星方位向成像时间ty,计算公式如下式所示:
ty=t0+tAzi·icenter
其中,ty为卫星方位向成像时间;icenter为视SAR影像中心点的方位向坐标;t0为卫星方位向成像起始时间,tAzi为方位向时间间隔,t0和tAzi均由表2提供。
S3.3.1.2,根据所述卫星方位向成像时间ty计算卫星位置矢量(XS,YS,ZS)(具体计算过程在步骤S3.1.1.1中已描述,不再赘述);
S3.3.1.3,根据所述卫星方位向成像时间ty计算速度矢量(具体计算过程在步骤S3.1.1.1中已描述,不再赘述);
S3.3.1.4,利用所述卫星位置矢量(XS,YS,ZS)、卫星速度矢量和外接矩形数字高程模型(RefDEM)数据中的每个像素点的地心空间直角坐标(XP,YP,ZP),根据多普勒方程计算得到多普勒值fd(多普勒方程中的雷达波长、近斜距和距离向分辨率参数由表2提供),同时计算得到多普勒变化率fd',结合多普勒测量值fdop,根据下式计算出需要调整的时间变量Δt:
Δt=-(fd-fdop)/f′d;
其中,多普勒变化率f′d是通过数值微分的方法近似求出,假设要计算当前方位向坐标y对应的时间ty的f′d,设dt=0.01s,利用多普勒方程分别计算ty和ty+dt处的多普勒值fd(ty)和fd(ty+dt),则
f′d=(fd(ty+dt)-fd(ty))/dt
此外,外接矩形数字高程模型(RefDEM)数据中的每个像素点的地心空间直角坐标(XP,YP,ZP),它是根据公式(7)计算得到的(BRefDEM,LRefDEM,HDEM),其中BRefDEM为该点的纬度,LRefDEM为该点的经度,HDEM为该点的高程。随后对(BRefDEM,LRefDEM,HDEM)进行坐标转换,转换到地心空间直角坐标系中得到(XP,YP,ZP)。
S3.3.1.5,根据下式,利用所述时间变量Δt对卫星成像时间ty进行修改,得到新的成像时刻ty+1:
ty+1=ty+Δt;
S3.3.1.6,重复步骤S3.3.1.2-S3.3.1.5,进行迭代计算,直到Δt绝对值小于阈值ε,迭代结束,输出卫星方位向成像时间ty;
S3.3.1.7,由解算的卫星方位向成像时间ty,根据下式计算每个像素点的初始方位向坐标y:y=(ty-t0)/Δty,其中,t0为卫星方位向成像起始时间,Δty为方位向时间间隔,这两个参量均由表2提供。
S3.3.2,根据公式(2)的距离方程计算外接矩形数字高程模型(RefDEM)数据中的每个像素点的初始距离向坐标x。
由距离方程结合步骤S3.3.1计算得到的每个像素点的位置矢量(XP,YP,ZP),以及由表2提供的近地点斜距R0、距离向分辨率ΔR,计算每个像素点的初始距离向坐标x。
S3.3.3,组合步骤S3.3.1和步骤S3.3.2得到的外接矩形数字高程模型(RefDEM)数据中的每个像素点的初始方位向坐标y和初始距离向坐标x,得到外接矩形数字高程模型(RefDEM)数据中的每个像素点在多视SAR影像雷达坐标空间中的像点坐标(x,y),所述像点坐标(x,y)为外接矩形数字高程模型(RefDEM)数据与多视SAR影像雷达坐标空间一一映射的坐标关系。通过这个映射关系,即可在两个数据(坐标系)中相互找到对应关系。
步骤S3.4,根据SAR影像模拟公式(即半经验后向散射模型),根据外接矩形数字高程模型(RefDEM)数据模拟SAR像元(i,j)处的影像值PixSim(i,j),最终所有模拟SAR像元值PixSim(i,j)的集合构成DEM模拟SAR影像。SAR影像模拟公式为:
其中,PixSim(i,j)为模拟SAR像元(i,j)处的影像值,σ为雷达散射截面,K为常数,η为当地入射角大小,其中,Rts为雷达入射方向矢量;为地面目标点所在表面的法线矢量。
PixSim(i,j)为模拟SAR像元(i,j)处的影像值,只是模拟SAR图像中的一个像素值,所有PixSim(i,j)组成的像素集合为DEM模拟SAR影像,如图3(C)所示。.
优选的,K=[1,10]区间内的数。
输出三个数据文件,如图3所示。图3-(a)中保存外接矩形数字高程模型(RefDEM)高程数据;因外接矩形DEM数据方向与多视SAR影像(雷达坐标)方向不一致,需要经过坐标映射才能将它们之间的像元一一对应起来,因此3-(b)中保存的是外接矩形DEM数据与多视SAR雷达坐标空间的对应关系,根据坐标映射数据即可实现相互转换;图3-(c)中保存的是基于DEM和公式(8)得到的所有像元值构成DEM模拟SAR影像,同样可以根据坐标映射数据,将DEM模拟SAR影像与多视SAR影像雷达坐标空间一一对应起来。
步骤S4,根据步骤S3.3获取的DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系和步骤S3.4获取的DEM模拟SAR影像,计算得到与多视SAR影像雷达坐标空间一致的SAR模拟影像,结果如图4a所示。
如图2(b)和3(c)所示,DEM模拟SAR影像与多视SAR影像之间的数据方向是不一致的,需要通过DEM模拟SAR影像和多视SAR影像雷达坐标空间一一映射的坐标关系(即坐标映射数据)实现数据转换,得到与多视SAR影像(雷达坐标)方向一致的SAR模拟影像。遍历坐标映射数据,取其复数像元值,若其复数像元值满足多视SAR影像的宽nMLWidth和高nMLHeight,则将对应的DEM模拟SAR影像的像元值赋值给SAR模拟影像。
步骤S4具体包括以下三个步骤:
步骤S4.1,获得DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系(即坐标映射数据1)的每一个像元pix(i,j)的复数像元值(iValue,jValue),如满足如下条件式,则保留该像元pix(i,j),参与下一步计算。
其中,ivalue为某一个像元pix(i,j)对应的复数像元值的高度方向坐标,jvalue为某一个像元pix(i,j)对应的复数像元值的宽度方向坐标,nMLHeight为多视SAR影像的高,nMLWidth为多视SAR影像的宽。
步骤S4.2,根据步骤S4.1获取的像元pix(i,j)和对应的复数像元值(iValue,jValue),在DEM模拟SAR影像中取(i,j)像元值填充到复数像元值(iValue,jValue)对应的SAR模拟影像坐标中。
步骤S4.3,对DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系(即坐标映射数据1)的每一个像元均采用步骤S4.1和步骤S4.2,得到SAR模拟影像,结果如图4-(a)所示。图4-(b)为多视SAR影像,现在它们的影像方向与大小是一致的,可以参与下一步的影像匹配。
步骤S5,对所述SAR模拟影像和所述多视SAR影像进行灰度匹配,计算得到新的DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系;
具体的,根据步骤S4得到的SAR模拟影像(简称SSAR)和多视SAR影像(简称MSAR),对SAR模拟影像与多视SAR影像在一定范围的窗口内进行灰度匹配,得到方位向和距离向的偏移量,并基于所述方位向和距离向的偏移量更新DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系(即坐标映射数据1)得到新的DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系(即坐标映射数据2)。
步骤S5,具体包括以下三个子步骤:
步骤S5.1,采用匹配方法,在SAR模拟影像上规则选取一定数量的点,对搜索窗口数为nWin范围内的像素进行匹配,得到窗口内每个像元的匹配测度ρ,取最大的匹配测度ρ,以及对应的方位向偏移量记为off_Az,距离向偏移量记为off_Rg。匹配方法采用常用的NCC匹配方法(归一化互相关系数),其公式如下:
其中,ri表示多视SAR影像(参考影像)窗口内的像素值;mi表示SAR模拟影像(待匹配影像)窗口内的像素值;表示多视SAR影像窗口内的像素均值;表示SAR模拟影像窗口内的像素均值;σr表示多视SAR影像窗口内的像素标准差;σm表示SAR模拟影像窗口内的像素标准差;n表示匹配窗口内的像素个数,n=nWin2,nWin为搜索窗口大小。
优选的,nWin=512。
步骤S5.2,对步骤S5.1获取的方位向的偏移量off_Az、距离向的偏移量off_Rg,进行多项式拟合,得到方位向和距离向的偏移量拟合系数。多项式拟合公式为:
其中,y(xAz)、y(xRg)为方位向和距离向的偏移量矩阵,矩阵大小为k个元素;xAz、xRg为对应的元素索引,取值为xAz=1,2,...,k,xRg=1,2,...,k;a0、a1为方位向的多项式常数项和一次项系数;b0、b1为距离向的多项式常数项和一次项系数。
优选的,方位向偏移量系数为常数项a0=0.77653、一次项a1=1.19188e-005;距离向偏移量系数为常数项b0=0.09955、一次项b1=2.58910e-004。
步骤S5.3,根据步骤S5.2获取的方位向和距离向的偏移量拟合系数,更新坐标映射数据1得到新的坐标映射数据2。更新的公式如下:
对坐标映射数据1的每一个像元(i,j),得到其复数像元值(iValue,1,jValue,1),根据方位向和距离向的偏移量,有
其中,iValue,2、jValue,2为更新后的坐标映射数据2的方位向和距离向像元值;iValue,1、jValue,1为坐标映射数据1的方位向和距离向像元值;a0、a1为方位向的多项式常数项和一次项系数;b0、b1为距离向的多项式常数项和一次项系数。
步骤S6,基于所述新的DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系,计算资源三号卫星影像控制点数据在多视SAR影像中的像点坐标集;
具体的,根据步骤S3.3获取的外接矩形数字高程模型(RefDEM)数据的起始角点地理坐标(Bref,1,Lref,1),结合资源三号卫星影像控制点数据ZY3,逐点计算ZY3数据在外接矩形数字高程模型(RefDEM)中的坐标(iDEM,jDEM),其次利用步骤S5.3获取的坐标映射数据2,得到ZY3数据在多视SAR影像中的像点坐标(iMSAR,jMSAR),构成ZY3控制点数据在多视SAR影像中的像点坐标集(Bi,Li,Hi,iMSAR,jMSAR)。
步骤S6具体包括以下三个子步骤:
步骤S6.1,根据步骤S3.3获取的外接矩形数字高程模型(RefDEM)数据的起始角点地理坐标(B1,L1),结合资源三号卫星影像控制点数据ZY3,见表1所示,逐点计算ZY3数据在外接矩形数字高程模型(RefDEM)中的坐标(iDEM,jDEM),计算公式如下所示:
其中,(Bref,1,Lref,1)为外接矩形数字高程模型(RefDEM)的起始点地理坐标;(Bi,Li)为资源三号卫星影像控制点数据ZY3的经纬度地理坐标,如表1所示;Bref,res为纬度方向外接矩形数字高程模型(RefDEM)分辨率,Lref,res为经度方向外接矩形数字高程模型(RefDEM)分辨率;(iDEM,jDEM)为ZY3数据在外接矩形RefDEM中的坐标。
优选的,外接矩形数字高程模型(RefDEM)的起始点地理坐标为(Bref,1,Lref,1)=(40.224,115.799)度;(Bi,Li)为表1中纬度和经度信息;外接矩形数字高程模型(RefDEM)的纬度方向分辨率为Bref,res=-8.3333e-005度;外接矩形数字高程模型(RefDEM)的经度方向分辨率为Lref,res=8.3333e-005度。
控制点名称 | 纬度/度 | 经度/度 | iDEM | jDEM |
F03 | 39.80890 | 116.09411 | 4981.22 | 3541.334 |
F18 | 39.74416 | 116.01044 | 5758.10 | 2537.29 |
F05 | 39.79243 | 116.16860 | 5178.86 | 4435.218 |
C24 | 39.88796 | 116.01093 | 4032.50 | 2543.17 |
C22 | 39.89176 | 116.02806 | 3986.90 | 2748.731 |
W15 | 40.12384 | 116.22403 | 1201.92 | 5100.38 |
C02 | 40.12662 | 116.22181 | 1168.56 | 5073.74 |
W17 | 40.06218 | 116.28779 | 1941.85 | 5865.503 |
表4 ZY3数据在外接矩形数字高程模型(RefDEM)中的坐标
步骤S6.2,利用步骤S6.1获取的ZY3数据在外接矩形数字高程模型(RefDEM)中的坐标(iDEM,jDEM),结合步骤S5.3获取的坐标映射数据2,计算(iDEM,jDEM)在SAR模拟影像中的像素坐标(iMSAR,jMSAR)。因为根据步骤S3获取的外接矩形数字高程模型(RefDEM)和坐标映射数据2,它们之间的像素关系是一一对应的,如图3所示,即仅仅需要把外接矩形数字高程模型(RefDEM)的像素坐标在坐标映射数据2中将Pix(iDEM,jDEM)的像素值输出即可,得到(iMSAR,jMSAR)。
优选的,以表1中的F03点为例,计算得到F03点在多视SAR影像中的像点坐标为(iMSAR,jMSAR)=(1285.46,5329.784)。
步骤S6.3,将所有ZY3数据重复步骤S6.1和S6.2,获得ZY3在SAR模拟影像中的像点坐标(iMSAR,jMSAR),同时由于SAR模拟影像已与多视SAR影像进行匹配,它们之间的像素关系是一一对应的,即ZY3在多视SAR影像中的像点坐标(iMSAR,jMSAR)是可以直接使用的,不需要进行任何转换,并基于所述像点坐标(iMSAR,jMSAR)和ZY3的地理坐标和高程构成ZY3控制点数据在多视SAR影像中的坐标集(Bi,Li,Hi,iMSAR,jMSAR)。其中,(Bi,Li,Hi)为ZY3的地理坐标和高程;(iMSAR,jMSAR)为ZY3在多视SAR影像中的像点坐标。
优选的,ZY3控制点数据在多视SAR影像中的坐标集(Bi,Li,Hi,iMSAR,jMSAR),如表5所示。
控制点名称 | 纬度/度 | 经度/度 | iMSAR | jMSAR |
F03 | 39.80890 | 116.09411 | 1285.46 | 5329.784 |
F18 | 39.74416 | 116.01044 | 1872.99 | 6389.175 |
F05 | 39.79243 | 116.16860 | 632.348 | 5414.823 |
C24 | 39.88796 | 116.01093 | 2134.48 | 4395.387 |
C22 | 39.89176 | 116.02806 | 1993.94 | 4309.637 |
W15 | 40.12384 | 116.22403 | 831.600 | 718.400 |
C02 | 40.12662 | 116.22181 | 856.207 | 684.1301 |
W17 | 40.06218 | 116.28779 | 180.838 | 1448.448 |
表5 ZY3数据在多视SAR影像中的坐标
步骤S7,将所述资源三号卫星影像控制点数据在多视SAR影像中的像点坐标集转换为所述资源三号卫星影像控制点数据在原始SAR影像中的坐标集;
根据步骤S2设置的SAR影像的方位向(行方向)多视数为nAzML,距离向(列方向)多视数为nRgML,需要将资源三号卫星影像控制点数据在多视SAR影像中的坐标集(Bi,Li,Hi,iMSAR,jMSAR)转换为资源三号卫星影像控制点数据在原始SAR影像(即,真实SAR影像)坐标集,即多视SAR影像的像素坐标分别乘以方位向多视数为nAzML,距离向多视数为nRgML,最后构成ZY3控制点数据在原始SAR影像中的坐标集(Bi,Li,Hi,iSAR,jSAR)。计算公式为:
其中,(iSAR,jSAR)为真实SAR像素坐标,(iMSAR,jMSAR)为多视SAR像素坐标,nAzML为方位向多视数,nRgML为距离向多视数。
优选的,多视数nAzML=4,nRgML=4;计算得到ZY3控制点数据在真实SAR影像中的坐标集(Bi,Li,Hi,iSAR,jSAR),如表6所示。
控制点名称 | 纬度/度 | 经度/度 | iSAR | jSAR |
F03 | 39.80890 | 116.09411 | 5141.87 | 21319.14 |
F18 | 39.74416 | 116.01044 | 7491.97 | 25556.70 |
F05 | 39.79243 | 116.16860 | 2529.39 | 21659.29 |
C24 | 39.88796 | 116.01093 | 8537.95 | 17581.55 |
C22 | 39.89176 | 116.02806 | 7975.77 | 17238.55 |
W15 | 40.12384 | 116.22403 | 3326.40 | 2873.600 |
C02 | 40.12662 | 116.22181 | 3424.83 | 2736.521 |
W17 | 40.06218 | 116.28779 | 723.354 | 5793.793 |
表6 ZY3数据在真实SAR影像中的坐标
步骤S8,利用资源三号卫星影像控制点数据在原始SAR影像中的坐标集,进行SAR影像正射纠正处理,计算得到距离向和方位向多项式的改正系数;
具体的,利用步骤S7获取的ZY3控制点数据集(Bi,Li,Hi,iSAR,jSAR)和步骤S1提供轨道状态信息,进行SAR影像正射纠正处理,计算得到距离向多项式的改正系数e0和e1,方位向多项式的改正系数f0和f1,参与后续计算。
SAR影像正射纠正的目的是建立目标点大地坐标与SAR影像坐标之间精确的映射关系。控制点的大地坐标(B,L,H)经过步骤S3.4构建的RD模型间接定位后得到控制点影像坐标(x',y'),这是利用定位模型得到的初步定位结果。显然对于一个精确的定位模型来说,(x',y')应该和实际测量的控制点影像坐标(x,y)是相同的。但由于直接利用SAR轨道参数信息计算得到的定位模型参数通常会存在一定的误差,实际上(x',y')与(x,y)并不相等。它们之间就建立起了目标点的初步定位结果(x',y')与其实际影像坐标(x,y)之间的关系。它们之间的几何关系表示为:
其中,(x,y)是控制点在影像上的实测坐标,(x',y')是利用RD模型由地理坐标计算得到的影像坐标,e0,e1分别为距离向多项式的常数项和一次项,f0,f1分别为方位向多项式的常数项和一次项。
步骤S8具体包括以下四个子步骤:
步骤S8.1,构建多项式误差模型,分别对SAR影像的方位向和距离向进行多项式系数的计算。考虑到SAR成像时,距离向和方位向近似垂直,两个方向上的误差相对独立,因此在影像的方位向和距离向上分别对影像坐标进行一次线性改正,表达式为:
其中,F(x)为距离向多项式函数,F(y)为方位向多项式函数,(x,y)是控制点在影像上的实测坐标,(x',y')是利用RD模型由地理坐标计算得到的影像坐标,e0,e1分别为距离向多项式的常数项和一次项,f0,f1分别为方位向多项式的常数项和一次项。
步骤S8.2,根据步骤S8.1的式(10)可以对ZY3控制点列如下线性方程,根据最小二乘平差求解影像的一次线性改正参数。解算误差方程如下:
其中,vx、vy分别为函数F(x)、F(y)的残差;为函数F(x)对e0的偏导数,为函数F(x)对e1的偏导数;为函数F(y)对f0的偏导数,为函数F(y)对f1的偏导数;Δe0为距离向多项式的常数项改正量,Δe1为距离向多项式的一次项改正量;Δf0为方位向多项式的常数项改正量,Δf1为方位向多项式的一次项改正量;Fx为已知的距离向坐标,Fy为已知的方位向坐标。
步骤S8.3,根据式(11)所列出的误差模型,列出距离向x的误差方程式
V1=A1Δx1-l1 (18)
其中,
Δx1=[Δe0 Δe1]T
其中,V1为残差项矩阵,其中vxG(i)表示第G(i)个控制点在F(x)函数上的残差,mG为ZY3控制点的个数;A1为距离向的偏导数系数矩阵,其中为函数F(x)对e0的偏导数,为函数F(x)对e1的偏导数;l1为常数项矩阵,其中Fx0G(i)表示第G(i)个控制点在F(x)函数上的模型计算值;Δx1为待求解的参数改正量矩阵,Δe0,Δe1分别为距离向多项式的常数项改正量和一次项改正量。
优选的,距离向多项式的常数项为e0=0.0741和一次项为e1=0.999985。
步骤S8.4,根据式(11)所列出的误差模型,列出方位向y的误差方程式
V2=A2Δx2-l2 (19)
其中,
Δx2=[Δf0 Δf1]T
其中,V2为残差项矩阵,其中vyG(i)表示第G(i)个控制点在F(y)函数上的残差,mG为ZY3控制点的个数;A2为距离向的偏导数系数矩阵,其中为函数F(y)对f0的偏导数,为函数F(y)对f1的偏导数;l2为常数项矩阵,其中Fy0G(i)表示第G(i)个控制点在F(y)函数上的模型计算值;Δx2为待求解的参数改正量矩阵,Δf0,Δf1分别为方位向多项式的常数项改正量和一次项改正量。
优选的,方位向多项式的常数项为f0=-0.2013和一次项为f1=1.000015。
步骤S9,根据所述距离向和方位向多项式的改正系数,对原始SAR影像进行SAR正射纠正,得到SAR正射纠正后的DOM影像。
具体的,根据步骤S8的公式(15),结合步骤S3.4构建的RD模型,对原始SAR影像的每个像素进行校正,输出数字正射影像图(Digital Orthophoto Map,DOM)影像和定位精度报告。
步骤S9具体包括以下四个子步骤:
步骤S9.1,根据原始SAR影像的四个角点像素坐标、原始SAR成像参数和输入的DOM数据的纬度方向分辨率、经度方向分辨率,确定待校正的DOM数据的宽和高,以及DOM起始点地理坐标。步骤S9.1包括以下子步骤:
步骤S9.1.1,根据原始SAR影像的四个角点像素坐标和原始SAR成像参数,利用公式(2)和步骤S3.1.1,得到原始SAR影像的四个角点像素坐标对应的四个地理坐标(计算过程和步骤S3.1.1类似,这里不再赘述)。
优选的,原始SAR影像的四个角点像素坐标对应的地理坐标分别为左上角(40.160378,116.333731)、右上角(40.216899,115.925287)、左下角(39.658992,116.212917)、右下角(39.720849,115.815612)。
步骤S9.1.2,根据四个地理坐标计算得到纬度方向的最大最小值和经度方向的最大最小值,进而得到纬度方向的纬度差值dBDOM和经度方向的经度差值dLDOM和DOM起始点地理坐标(BDOM,LDOM)。
优选的,纬度方向的纬度差值dBDOM=0.557907,经度方向的经度差值dLDOM=0.518119,DOM起始点地理坐标为纬度方向最大值和经度方向最小值得组合,即(BDOM,LDOM)=(40.216899,115.815612)。
步骤S9.1.3,根据步骤步骤S9.1.2获取的纬度方向的纬度差值dBDOM和经度方向的经度差值dLDOM,计算得到待校正的DOM数据的宽nWidthDOM和高nHeightDOM,公式如下:
nWidthDOM=dBDOM/BDOM,res
nHeightDOM=dLDOM/LDOM.res
其中,nWidthDOM为待校正的DOM数据的宽,nHeightDOM为待校正的DOM数据的高;dBDOM为纬度方向的纬度差值,dLDOM为经度方向的经度差值;BDOM,res为DOM数据的纬度方向分辨率,LDOM.res为DOM数据的经度方向分辨率。
优选的,输入的DOM数据的纬度方向分辨率和经度方向分辨率均为0.000052度;待校正的DOM数据的宽nWidthDOM=10002和高nHeightDOM=10770。
步骤S9.2,根据步骤S9.1的地理参数(包括DOM起始点地理坐标、输入的DOM数据的纬度方向分辨率、经度方向分辨率),利用步骤S3.4构建的RD模型,逐像素对待校正的DOM影像进行计算,得到待校正的DOM像素在原始SAR影像中的像素坐标(x',y')。
步骤S9.3,对步骤S9.2得到的像素坐标(x',y')进行改正,得到改正后的原始SAR像素坐标(x,y),获取原始SAR像素坐标(x,y)对应的像元值PixDOM(x,y),将所述像元值PixDOM(x,y)保存到步骤S9.2所述的待校正的DOM像素坐标pix(xDOM,yDOM)中,遍历整个待校正的DOM像素,获得DOM正射纠正影像。
具体的,根据步骤S8得到的方位向和距离向多项式系数和步骤S9.2得到的像素坐标(x',y'),利用步骤S8的公式(15)进行改正,得到改正后的原始SAR像素坐标(x,y),获取原始SAR像素坐标(x,y)对应的像元值PixDOM(x,y),将所述像元值PixDOM(x,y)保存到步骤S9.2所述的待校正的DOM像素坐标pix(xDOM,yDOM)中,遍历整个待校正的DOM像素,获得DOM正射纠正影像。
根据步骤S8得到的方位向和距离向多项式系数和步骤S9.2得到的像素坐标(x',y'),利用步骤S8的公式(15)进行改正,得到改正后的SAR像素坐标(x,y),如公式下式所示(公式15)
x=e0+e1·x'
y=f0+f1·y'
就是说,已知e0,e1分别为距离向多项式的常数项和一次项,f0,f1分别为方位向多项式的常数项和一次项,将通过步骤S9.2得到的像素坐标(x',y'),带入上述公式计算即可得到改正后的SAR像素坐标(x,y)。
优选的,步骤S9还可以包括:
步骤S9.4,对步骤S9.3得到的DOM正射纠正影像进行定位精度评价,输出校正后的精度报告。
优选的,根据表1提供的数据,对步骤S9.3输出的DOM正射纠正数据进行精度评价,统计得到纬度方向中误差为3.95m,经度方向中误差为3.01m,总的平面定位中误差为4.97m。对正射纠正影像DOM的定位精度进行统计,如表7所示。
表7 正射纠正影像DOM的定位精度统计
相对于现有技术,本发明通过外部DEM模拟SAR影像得到SAR模拟影像,并与多视SAR影像进行匹配,确定了高精度资源三号控制点库在SAR影像上的坐标集,比人工刺点精度要高;常规方法使用DEM参与SAR正射纠正精度较低,而现在通过模拟匹配的方法,使用高精度资源三号控制点库参与正射纠正,既解决了选点问题又达到了较好的纠正精度。通过匹配的方法,确定了高精度资源三号控制点库在SAR影像上的坐标集。
以上内容仅为本发明的较佳实施例,对于本领域的普通技术人员,依据本发明的思想,在具体实施方式及应用范围上均会有改变之处,本说明书内容不应理解为对本发明的限制。
Claims (10)
1.一种基于卫星控制点库和DEM的SAR影像正射纠正方法,其特征在于,所述方法具体包括以下步骤:
步骤S1,采集原始数据,所述原始数据优选的包括资源三号卫星影像控制点数据、外部参考DEM数据、原始SAR影像以及原始SAR影像成像参数;
步骤S2,对所述原始SAR影像进行多视处理,得到多视SAR影像;
步骤S3,基于所述外部参考DEM数据和多视SAR影像,计算外接矩形数字高程模型(RefDEM)数据、DEM模拟SAR影像、DEM模拟SAR影像和多视SAR影像雷达坐标空间一一映射的坐标关系;
步骤S4,根据步骤S3.3获取的DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系和步骤S3.4获取的DEM模拟SAR影像,计算得到与多视SAR影像雷达坐标空间一致的SAR模拟影像;
步骤S5,对所述SAR模拟影像和所述多视SAR影像进行灰度匹配,计算得到新的DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系;
步骤S6,基于所述新的DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系,计算资源三号卫星影像控制点数据在多视SAR影像中的像点坐标集;
步骤S7,将所述资源三号卫星影像控制点数据在多视SAR影像中的像点坐标集转换为所述资源三号卫星影像控制点数据在原始SAR影像中的坐标集;
步骤S8,利用资源三号卫星影像控制点数据在原始SAR影像中的坐标集,进行SAR影像正射纠正处理,计算得到距离向和方位向多项式的改正系数;
步骤S9,根据所述距离向和方位向多项式的改正系数,对所述原始SAR影像进行SAR正射纠正,得到SAR正射纠正后的DOM影像。
2.根据权利要求1所述的基于卫星控制点库和DEM的SAR影像正射纠正方法,其特征在于,步骤2具体包括:根据下式对所述原始SAR影像进行多视处理,得到多视SAR影像:
<mrow>
<mi>R</mi>
<mi>i</mi>
<mi>x</mi>
<mrow>
<mo>(</mo>
<mi>i</mi>
<mo>,</mo>
<mi>j</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<munderover>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>n</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mrow>
<mi>m</mi>
<mo>=</mo>
<mi>n</mi>
<mi>R</mi>
<mi>g</mi>
<mi>M</mi>
<mi>L</mi>
</mrow>
</munderover>
<mrow>
<mi>m</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mrow>
<mi>n</mi>
<mo>=</mo>
<mi>n</mi>
<mi>A</mi>
<mi>z</mi>
<mi>M</mi>
<mi>L</mi>
</mrow>
</munderover>
<mi>P</mi>
<mi>i</mi>
<mi>x</mi>
<mrow>
<mo>(</mo>
<mi>n</mi>
<mo>,</mo>
<mi>m</mi>
<mo>)</mo>
</mrow>
<mo>/</mo>
<mrow>
<mo>(</mo>
<mi>n</mi>
<mi>A</mi>
<mi>z</mi>
<mi>M</mi>
<mi>L</mi>
<mo>*</mo>
<mi>n</mi>
<mi>R</mi>
<mi>g</mi>
<mi>M</mi>
<mi>L</mi>
<mo>)</mo>
</mrow>
</mrow>
其中,Pix(i,j)表示多视SAR影像中的像元(i,j)对应的像元值;nAzML表示原始SAR影像的方位向多视数,nRgML表示原始SAR影像的距离向多视数;Pix(n,m)表示原始SAR影像中的像元(n,m)对应的像元值,n表示原始SAR影像中的像元(n,m)的方位向坐标,取值范围为n=1,...,nAzML,m表示原始SAR影像中的像元(n,m)的距离向坐标,取值范围为m=1,...,nRgML。
3.根据权利要求1所述的基于卫星控制点库和DEM的SAR影像正射纠正方法,其特征在于,所述步骤3优选的包括以下子步骤:
步骤S3.1,基于所述外部参考DEM数据和多视SAR影像,计算外接矩形数字高程模型(RefDEM)数据,其中,所述外接矩形数字高程模型(RefDEM)数据包括所述外接矩形数字高程模型(RfDEM)在外部参考DEM数据中的宽、高、起始点地理坐标、纬度方向分辨率、经度方向分辨率和高程值;
步骤S3.2,根据步骤S3.1获取的外接矩形数字高程模型(RefDEM)的宽nWidthDEM和高nHeightDEM,创建两个大小均为(nHeightDEM,nWidthDEM)的实数空矩阵,用于分别存储外接矩形数字高程模型(RefDEM)数据和DEM模拟SAR影像数据;创建一个大小为(nHeightDEM,nWidthDEM)的复数空矩阵,用于存储DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系数据;
步骤S3.3,基于所述外部参考DEM数据和多视SAR影像,根据SAR严密成像模型(R-D模型),计算外接矩形数字高程模型(RefDEM)数据与多视SAR影像雷达坐标空间一一映射的坐标关系;
步骤S3.4,根据SAR影像模拟公式,根据外接矩形数字高程模型(RefDEM)数据模拟SAR像元(i,j)处的影像值PixSim(i,j),最终所有模拟SAR像元值PixSim(i,j)的集合构成DEM模拟SAR影像。
4.根据权利要求1所述的基于卫星控制点库和DEM的SAR影像正射纠正方法,其特征在于,所述步骤S3.1可以包括:
步骤S3.1.1,利用外部参考DEM(ExtDEM)数据和多视SAR影像的轨道状态、成像参数信息,根据SAR严密成像模型(R-D模型)计算得到多视SAR影像的四个角点坐标在外部参考DEM(ExtDEM)数据中的像素坐标;
步骤S3.1.2,根据步骤S3.1.1获取的多视SAR影像的四个角点在外部参考DEM(ExtDEM)数据中的像素坐标(iDEM,jDEM),计算外接矩形数字高程模型(RefDEM)的四个角点坐标在外部参考DEM(ExtDEM)数据中的坐标,以及外接矩形数字高程模型(RefDEM)在外部参考DEM数据中的宽nWidthDEM和高nHeightDEM;
步骤S3.1.3,根据步骤S3.1.2得到的外接矩形数字高程模型(RefDEM)的左上角Clt(i,j)角点坐标,结合外部参考DEM(ExtDEM)数据提供的起点地理坐标(Bext,0,Lext,0)、纬度方向的分辨率Bext,res、经度方向的分辨率Lext,res,根据下式计算外接矩形数字高程模型(RefDEM)的起始点地理坐标(B1,L1);
B1=Bext,0+Bext,res*i
L1=Lext,0+Lext,res*j
其中,(B1,L1)为外接矩形数字高程模型(RefDEM)的起始点地理坐标;(Bext,0,Lext,0)为外部参考DEM(ExtDEM)数据提供的起点地理坐标,Bext,res为外部参考DEM(ExtDEM)数据的纬度方向的分辨率,Lext,res为外部参考DEM(ExtDEM)数据的经度方向的分辨率;i为外接矩形数字高程模型(RefDEM)在外部参考DEM(ExtDEM)数据中的左上角的行方向坐标,j为外接矩形数字高程模型(RefDEM)在外部参考DEM(ExtDEM)数据中的左上角的列方向坐标;
步骤S3.1.4,根据步骤S3.1.3获取的外接矩形数字高程模型(RefDEM)的起始点地理坐标(B1,L1),计算外接矩形数字高程模型(RefDEM)的高程值。
5.根据权利要求4所述的基于卫星控制点库和DEM的SAR影像正射纠正方法,其特征在于,所述步骤S3.1可以包括以下子步骤:
步骤S3.1.1.1,根据多视SAR影像的角点的方位向坐标iAZ,计算卫星方位向成像时间tT,并根据卫星方位向成像时间tT计算卫星位置矢量(Xs,Ys,Zs)和速度矢量
步骤S3.1.1.2,对多视SAR影像的中心点经纬度进行坐标转换,将其从大地坐标系转换到地心空间直角坐标系,作为待求点的地心坐标初值(XT 0,YT 0,ZT 0);
步骤S3.1.1.3,基于步骤S3.1.1.1求得的卫星位置矢量(Xs,Ys,Zs)和速度矢量和步骤S3.1.1.2求得的待求点的地心坐标初值(XT 0,YT 0,ZT 0),根据SAR严密成像模型,解算所述待求点的改正数;
步骤S3.1.1.4,根据下式,利用待求点的改正数更新待求点的地心坐标:
XT 1=XT 0+ΔXT,YT 1=YT 0+ΔYT,ZT 1=ZT 0+ΔZT;
其中,是更新后的待求点的地心坐标,是待求点的地心坐标初值,ΔXT,ΔYT,ΔZT是待求点的改正坐标;
步骤S3.1.1.5,利用更新后的待求点的地心坐标(XT 1,YT 1,ZT 1)进行迭代解算,重复S3.1.1.3-S3.1.1.4,直至改正数(ΔXT,ΔYT,ΔZT)的绝对值都小于设定阈值ε迭代结束,确定待求点的地心坐标(XT,YT,ZT);
步骤S3.1.1.6,对通过步骤S3.1.1.5获取的待求点的地心坐标(XT,YT,ZT)进行坐标转换,将其从地心空间直角坐标系转换到大地坐标系,得到并输出待求点的大地坐标(BT,LT,HT);
步骤S3.1.1.7,根据下述公式,计算得到多视SAR影像的角点坐标在外部参考DEM(ExtDEM)数据中的像素坐标(iDEM,jDEM):
iDEM=(Bext,0-Bi)/Bext,res
jDEM=(Lext,0-Li)/Lext,res
其中,(iDEM,jDEM)为多视SAR影像的角点坐标在外部参考DEM(ExtDEM)数据中的像素坐标;(Bext,0,Lext,0)为外部参考DEM(ExtDEM)数据的起点地理坐标;Bext,res为纬度方向分辨率,Lext,res为经度方向分辨率;Bi为多视SAR影像的角点坐标对应的大地坐标纬度值,Li为多视SAR影像角点坐标对应的大地坐标经度值。
6.根据权利要求5所述的基于卫星控制点库和DEM的SAR影像正射纠正方法,其特征在于,步骤S3.1.1.3优选的包括:
所述SAR严密成像模型为:
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<msqrt>
<mrow>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>X</mi>
<mi>S</mi>
</msub>
<mo>-</mo>
<msub>
<mi>X</mi>
<mi>P</mi>
</msub>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>+</mo>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>Y</mi>
<mi>S</mi>
</msub>
<mo>-</mo>
<msub>
<mi>Y</mi>
<mi>P</mi>
</msub>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<mo>+</mo>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>Z</mi>
<mi>S</mi>
</msub>
<mo>-</mo>
<msub>
<mi>Z</mi>
<mi>P</mi>
</msub>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
</msqrt>
<mo>=</mo>
<msub>
<mi>R</mi>
<mn>0</mn>
</msub>
<mo>+</mo>
<mi>&Delta;</mi>
<mi>R</mi>
<mo>&CenterDot;</mo>
<mi>x</mi>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>Vx</mi>
<mi>S</mi>
</msub>
<mo>&CenterDot;</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>X</mi>
<mi>S</mi>
</msub>
<mo>-</mo>
<msub>
<mi>X</mi>
<mi>P</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>+</mo>
<msub>
<mi>Vy</mi>
<mi>S</mi>
</msub>
<mo>&CenterDot;</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>Y</mi>
<mi>S</mi>
</msub>
<mo>-</mo>
<msub>
<mi>Y</mi>
<mi>P</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>+</mo>
<msub>
<mi>Vz</mi>
<mi>S</mi>
</msub>
<mo>&CenterDot;</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>Z</mi>
<mi>S</mi>
</msub>
<mo>-</mo>
<msub>
<mi>Z</mi>
<mi>P</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mo>-</mo>
<msub>
<mi>f</mi>
<mi>d</mi>
</msub>
<mo>&CenterDot;</mo>
<mi>&lambda;</mi>
<mo>&CenterDot;</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>R</mi>
<mn>0</mn>
</msub>
<mo>+</mo>
<mi>&Delta;</mi>
<mi>R</mi>
<mo>&CenterDot;</mo>
<mi>x</mi>
<mo>)</mo>
</mrow>
<mo>/</mo>
<mn>2</mn>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
其中,VxS·(XS-XP)+VyS·(YS-YP)+VzS·(ZS-ZP)=-fd·λ·(R0+ΔR·x)/2为多普勒方程,为距离方程,(XS,YS,ZS)为卫星位置矢量,(XP,YP,ZP)为角点的地心空间直角坐标,为卫星速度矢量,x为角点的距离向坐标,R0为近地点斜距,ΔR为距离向分辨率,fd为多普勒中心频率,λ为雷达波长;
根据上述公式(2),构建观测方程,计算多视SAR影像的角点的误差方程系数,形成误差方程系数矩阵,再根据误差方程系数矩阵建立法方程:BTBX-BTL=0,根据下式解算待求点的改正数:
X=(BTB)-1BTL
其中,X=(ΔXT,ΔYT,ΔZT),ΔXT,ΔYT,ΔZT是待求点的改正数坐标;B为误差方程系数矩阵;L为常数项矩阵;;
其中,,在步骤S3.1.4中,根据下述公式计算外接矩形数字高程模型(RefDEM)的高程值:
BRefDEM=B1+Bref,res*m,m=0,1,…,nHeightDEM-1
LRefDEM=L1+Lref,res*n,n=0,1,…,nWidthDEM-1
iDEM=(Bext,0-BRefDEM)/Bext,res
jDEM=(Lext,0-LRefDEM)/Lext,res
HDEM=PixDEM(iDEM,jDEM)
其中,(BRefDEM,LRefDEM)为外接矩形数字高程模型(RefDEM)的某个点的经纬度坐标,BRefDEM为外接矩形数字高程模型(RefDEM)的某个点的纬度坐标,LRefDEM为外接矩形数字高程模型(RefDEM)的某个点的经度坐标;(B1,L1)为外接矩形数字高程模型(RefDEM)的起始点地理坐标,Bref,res为外接矩形数字高程模型(RefDEM)的纬度方向的分辨率,Lref,res为外接矩形数字高程模型(RefDEM)的经度方向的分辨率,nWidthDEM和nHeightDEM分别为外接矩形数字高程模型(RefDEM)的宽和高,m为高度方向的坐标索引,其中m=0,1,…,nHeightDEM-1;n为宽度方向的坐标索引,其中n=0,1,…,nWidthDEM-1;(Bext,0,Lext,0)为外部参考DEM(ExtDEM)数据的起点地理坐标,Bext,res为外部参考DEM(ExtDEM)在纬度方向的分辨率,Lext,res为外部参考DEM(ExtDEM)在经度方向的分辨率;(iDEM,jDEM)为外接矩形数字高程模型(RefDEM)的某个点的经纬度坐标在外部参考DEM(ExtDEM)数据中的行列坐标;PixDEM(iDEM,jDEM)表示行列坐标为(iDEM,jDEM)的点在外部参考DEM(ExtDEM)数据中的高程值,HDEM表示外接矩形数字高程模型(RefDEM)的高程值。
7.根据权利要求3所述的影像正射纠正方法,其特征在于,步骤S3.3优选的包括以下子步骤:
步骤S3.3.1,基于公式(2)中的多普勒方程迭代计算外接矩形数字高程模型(RefDEM)数据中的每个像素点的初始方位向坐标y;
S3.3.2,根据公式(2)的距离方程计算外接矩形数字高程模型(RefDEM)数据中的每个像素点的初始距离向坐标x;
S3.3.3,组合步骤S3.3.1和步骤S3.3.2得到的外接矩形数字高程模型(RefDEM)数据中的每个像素点的初始方位向坐标y和初始距离向坐标x,得到外接矩形数字高程模型(RefDEM)数据中的每个像素点在多视SAR影像雷达坐标空间中的像点坐标(x,y),所述像点坐标(x,y)为外接矩形数字高程模型(RefDEM)数据与多视SAR影像雷达坐标空间一一映射的坐标关系。
8.根据权利要求1所述的影像正射纠正方法,其特征在于,
步骤S4具体包括以下三个步骤:
步骤S4.1,获得DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系(即坐标映射数据1)的每一个像元pix(i,j)的复数像元值(iValue,jValue),如满足如下条件式,则保留该像元pix(i,j),参与下一步计算:
<mrow>
<mtable>
<mtr>
<mtd>
<mrow>
<mn>0</mn>
<mo><</mo>
<mo>=</mo>
<msub>
<mi>i</mi>
<mrow>
<mi>v</mi>
<mi>a</mi>
<mi>l</mi>
<mi>u</mi>
<mi>e</mi>
</mrow>
</msub>
<mo><</mo>
<mi>n</mi>
<mi>M</mi>
<mi>L</mi>
<mi>H</mi>
<mi>e</mi>
<mi>i</mi>
<mi>g</mi>
<mi>h</mi>
<mi>t</mi>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mn>0</mn>
<mo><</mo>
<mo>=</mo>
<msub>
<mi>j</mi>
<mrow>
<mi>v</mi>
<mi>a</mi>
<mi>l</mi>
<mi>u</mi>
<mi>e</mi>
</mrow>
</msub>
<mo><</mo>
<mi>n</mi>
<mi>M</mi>
<mi>L</mi>
<mi>W</mi>
<mi>i</mi>
<mi>d</mi>
<mi>t</mi>
<mi>h</mi>
</mrow>
</mtd>
</mtr>
</mtable>
<mo>;</mo>
</mrow>
其中,ivalue为某一个像元pix(i,j)对应的复数像元值的高度方向坐标,jvalue为某一个像元pix(i,j)对应的复数像元值的宽度方向坐标,nMLHeight为多视SAR影像的高,nMLWidth为多视SAR影像的宽;
步骤S4.2,根据步骤S4.1获取的像元pix(i,j)和对应的复数像元值(iValue,jValue),在DEM模拟SAR影像中取(i,j)像元值填充到复数像元值(iValue,jValue)对应的SAR模拟影像坐标中;
步骤S4.3,对DEM模拟SAR影像与多视SAR影像雷达坐标空间一一映射的坐标关系(即坐标映射数据1)的每一个像元均采用步骤S4.1和步骤S4.2,得到SAR模拟影像;
步骤S5具体包括以下三个子步骤:
步骤S5.1,采用匹配方法,在SAR模拟影像上规则选取一定数量的点,对搜索窗口数为nWin范围内的像素进行匹配,得到窗口内每个像元的匹配测度ρ,取最大的匹配测度ρ,以及对应的方位向偏移量记为off_Az,距离向偏移量记为off_Rg。匹配方法采用常用的NCC匹配方法(归一化互相关系数),其公式如下:
<mrow>
<mi>&rho;</mi>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mi>n</mi>
</mfrac>
<mfrac>
<mrow>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>n</mi>
</munderover>
<mrow>
<mo>(</mo>
<msub>
<mi>r</mi>
<mi>i</mi>
</msub>
<mo>-</mo>
<mover>
<mi>r</mi>
<mo>&OverBar;</mo>
</mover>
<mo>)</mo>
</mrow>
<mrow>
<mo>(</mo>
<msub>
<mi>m</mi>
<mi>i</mi>
</msub>
<mo>-</mo>
<mover>
<mi>m</mi>
<mo>&OverBar;</mo>
</mover>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<msub>
<mi>&sigma;</mi>
<mi>r</mi>
</msub>
<msub>
<mi>&sigma;</mi>
<mi>m</mi>
</msub>
</mrow>
</mfrac>
</mrow>
其中,ri表示多视SAR影像(参考影像)窗口内的像素值;mi表示SAR模拟影像(待匹配影像)窗口内的像素值;表示多视SAR影像窗口内的像素均值;表示模拟SAR影像窗口内的像素均值;σr表示多视SAR影像窗口内的像素标准差;σm表示模拟SAR影像窗口内的像素标准差;n表示匹配窗口内的像素个数,n=nWin2,nWin为搜索窗口大小;
步骤S5.2,对步骤S5.1获取的方位向的偏移量off_Az、距离向的偏移量off_Rg,进行多项式拟合,得到方位向和距离向的偏移量拟合系数;多项式拟合公式为:
y(xAz)=a0+a1xAz
y(xRg)=b0+b1xRg
其中,y(xAz)、y(xRg)为方位向和距离向的偏移量矩阵,矩阵大小为k个元素;xAz、xRg为对应的元素索引,取值为xAz=1,2,...,k,xRg=1,2,...,k;a0、a1为方位向的多项式常数项和一次项系数;b0、b1为距离向的多项式常数项和一次项系数;
步骤S5.3,根据步骤S5.2获取的方位向和距离向的偏移量拟合系数,更新坐标映射数据1得到新的坐标映射数据2。更新的公式如下:
对坐标映射数据1的每一个像元(i,j),得到其复数像元值(iValue,1,jValue,1),根据方位向和距离向的偏移量,有
iValue,2=iValue,1+a0+a1·iValue,1
jValue,2=jValue,1+b0+b1·jValue,1
其中,iValue,2、jValue,2为更新后的坐标映射数据2的方位向和距离向像元值;iValue,1、jValue,1为坐标映射数据1的方位向和距离向像元值;a0、a1为方位向的多项式常数项和一次项系数;b0、b1为距离向的多项式常数项和一次项系数。
9.根据权利要求1所述的影像正射纠正方法,其特征在于,步骤S6具体包括以下子步骤:
步骤S6.1,计算ZY3数据在外接矩形数字高程模型(RefDEM)中的坐标(iDEM,jDEM),计算公式如下所示:
iDEM=(Bi-Bref,1)/Bref,res
jDEM=(Li-Lref,1)/Lref,res
其中,(Bref,1,Lref,1)为外接矩形数字高程模型(RefDEM)的起始点地理坐标;(Bi,Li)为资源三号卫星影像控制点数据ZY3的经纬度地理坐标;Bref,res为纬度方向外接矩形数字高程模型(RefDEM)分辨率,Lref,res为经度方向外接矩形数字高程模型(RefDEM)分辨率;(iDEM,jDEM)为ZY3数据在外接矩形RefDEM中的坐标;
步骤S6.2,计算(iDEM,jDEM)在SAR模拟影像中的像素坐标(iMSAR,jMSAR);
步骤S6.3,将所有ZY3数据重复步骤S6.1和S6.2,获得ZY3在SAR模拟影像中的像点坐标(iMSAR,jMSAR),并基于所述像点坐标(iMSAR,jMSAR)和ZY3的地理坐标和高程构成ZY3控制点数据在多视SAR影像中的坐标集(Bi,Li,Hi,iMSAR,jMSAR);其中,(Bi,Li,Hi)为ZY3的地理坐标和高程;(iMSAR,jMSAR)为ZY3在多视SAR影像中的像点坐标。
10.根据权利要求9所述的影像正射纠正方法,其特征在于,步骤S9具体包括以下子步骤:
步骤S9.1,根据原始SAR影像的四个角点像素坐标、原始SAR成像参数和输入的DOM数据的纬度方向分辨率、经度方向分辨率,确定待校正的DOM数据的宽和高,以及DOM起始点地理坐标;
步骤S9.2,根据步骤S9.1的地理参数,利用步骤S3.4构建的RD模型,逐像素对待校正的DOM影像进行计算,得到待校正的DOM像素在原始SAR影像中的像素坐标(x',y');
步骤S9.3,对步骤S9.2得到的像素坐标(x',y')进行改正,得到改正后的原始SAR像素坐标(x,y),获取原始SAR像素坐标(x,y)对应的像元值PixDOM(x,y),将所述像元值PixDOM(x,y)保存到步骤S9.2所述的待校正的DOM像素坐标pix(xDOM,yDOM)中,遍历整个待校正的DOM像素,获得DOM正射纠正影像。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710554588.9A CN107341778B (zh) | 2017-07-10 | 2017-07-10 | 基于卫星控制点库和dem的sar影像正射纠正方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710554588.9A CN107341778B (zh) | 2017-07-10 | 2017-07-10 | 基于卫星控制点库和dem的sar影像正射纠正方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN107341778A true CN107341778A (zh) | 2017-11-10 |
CN107341778B CN107341778B (zh) | 2020-11-10 |
Family
ID=60218680
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201710554588.9A Active CN107341778B (zh) | 2017-07-10 | 2017-07-10 | 基于卫星控制点库和dem的sar影像正射纠正方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN107341778B (zh) |
Cited By (11)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108445457A (zh) * | 2018-02-12 | 2018-08-24 | 中国人民解放军61540部队 | 星载分布式干涉合成孔径雷达三维基线定标的方法 |
CN109063711A (zh) * | 2018-07-06 | 2018-12-21 | 航天星图科技(北京)有限公司 | 一种基于llts框架的卫星影像正射纠正算法 |
CN109886910A (zh) * | 2019-02-25 | 2019-06-14 | 榆林市国土资源局 | 外部数字高程模型dem修正方法及装置 |
CN109903352A (zh) * | 2018-12-24 | 2019-06-18 | 中国科学院遥感与数字地球研究所 | 一种卫星遥感影像大区域无缝正射影像制作方法 |
CN110133653A (zh) * | 2019-05-29 | 2019-08-16 | 中国空间技术研究院 | 一种基于dsm数据的星载sar图像快速间接定位方法 |
CN112433213A (zh) * | 2020-10-21 | 2021-03-02 | 中国电子科技集团公司第二十九研究所 | 一种sar干涉测量结果与光学影像位置偏移综合纠正方法 |
CN112862966A (zh) * | 2021-02-20 | 2021-05-28 | 中煤航测遥感集团有限公司 | 地表三维模型构建方法、装置、设备及存储介质 |
CN113902626A (zh) * | 2021-08-23 | 2022-01-07 | 桂林理工大学 | 超宽幅线阵影像附加约束条件的正射纠正方法 |
CN115423696A (zh) * | 2022-07-29 | 2022-12-02 | 上海海洋大学 | 一种自适应线程参数的遥感正射影像并行生成方法 |
CN115423702A (zh) * | 2022-08-23 | 2022-12-02 | 自然资源部国土卫星遥感应用中心 | 制作大区域星载光学和sar影像dom的方法和系统 |
CN117115621A (zh) * | 2023-10-24 | 2023-11-24 | 中国海洋大学 | 基于改进U-Net网络的卫星云图预测方法 |
Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2008133699A2 (en) * | 2006-07-20 | 2008-11-06 | Harris Corporation | Geospatial modeling system providing non-linear inpainting for voids in geospatial model frequency domain data and related methods |
US20090089018A1 (en) * | 2007-09-28 | 2009-04-02 | Harris Corporation | Geospatial modeling system providing building generation based upon user input on 3d model and related methods |
CN102968631A (zh) * | 2012-11-22 | 2013-03-13 | 中国科学院、水利部成都山地灾害与环境研究所 | 山区多光谱遥感卫星影像的自动几何纠正与正射校正方法 |
CN103218780A (zh) * | 2013-03-27 | 2013-07-24 | 中国科学院电子学研究所 | 基于逆rd定位模型的无控星载sar图像正射校正方法 |
US20140267250A1 (en) * | 2013-03-15 | 2014-09-18 | Intermap Technologies, Inc. | Method and apparatus for digital elevation model systematic error correction and fusion |
CN106097404A (zh) * | 2016-05-27 | 2016-11-09 | 山东科技大学 | 利用非线性矢量曲面构建InSAR相位图像模型的方法 |
-
2017
- 2017-07-10 CN CN201710554588.9A patent/CN107341778B/zh active Active
Patent Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2008133699A2 (en) * | 2006-07-20 | 2008-11-06 | Harris Corporation | Geospatial modeling system providing non-linear inpainting for voids in geospatial model frequency domain data and related methods |
US20090089018A1 (en) * | 2007-09-28 | 2009-04-02 | Harris Corporation | Geospatial modeling system providing building generation based upon user input on 3d model and related methods |
CN102968631A (zh) * | 2012-11-22 | 2013-03-13 | 中国科学院、水利部成都山地灾害与环境研究所 | 山区多光谱遥感卫星影像的自动几何纠正与正射校正方法 |
US20140267250A1 (en) * | 2013-03-15 | 2014-09-18 | Intermap Technologies, Inc. | Method and apparatus for digital elevation model systematic error correction and fusion |
CN103218780A (zh) * | 2013-03-27 | 2013-07-24 | 中国科学院电子学研究所 | 基于逆rd定位模型的无控星载sar图像正射校正方法 |
CN106097404A (zh) * | 2016-05-27 | 2016-11-09 | 山东科技大学 | 利用非线性矢量曲面构建InSAR相位图像模型的方法 |
Cited By (15)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108445457A (zh) * | 2018-02-12 | 2018-08-24 | 中国人民解放军61540部队 | 星载分布式干涉合成孔径雷达三维基线定标的方法 |
CN109063711A (zh) * | 2018-07-06 | 2018-12-21 | 航天星图科技(北京)有限公司 | 一种基于llts框架的卫星影像正射纠正算法 |
CN109063711B (zh) * | 2018-07-06 | 2021-10-29 | 中科星图股份有限公司 | 一种基于llts框架的卫星影像正射纠正算法 |
CN109903352A (zh) * | 2018-12-24 | 2019-06-18 | 中国科学院遥感与数字地球研究所 | 一种卫星遥感影像大区域无缝正射影像制作方法 |
CN109886910B (zh) * | 2019-02-25 | 2021-07-13 | 榆林市国土资源局 | 外部数字高程模型dem修正方法及装置 |
CN109886910A (zh) * | 2019-02-25 | 2019-06-14 | 榆林市国土资源局 | 外部数字高程模型dem修正方法及装置 |
CN110133653A (zh) * | 2019-05-29 | 2019-08-16 | 中国空间技术研究院 | 一种基于dsm数据的星载sar图像快速间接定位方法 |
CN110133653B (zh) * | 2019-05-29 | 2020-12-08 | 中国空间技术研究院 | 一种基于dsm数据的星载sar图像快速间接定位方法 |
CN112433213A (zh) * | 2020-10-21 | 2021-03-02 | 中国电子科技集团公司第二十九研究所 | 一种sar干涉测量结果与光学影像位置偏移综合纠正方法 |
CN112862966A (zh) * | 2021-02-20 | 2021-05-28 | 中煤航测遥感集团有限公司 | 地表三维模型构建方法、装置、设备及存储介质 |
CN112862966B (zh) * | 2021-02-20 | 2024-01-26 | 中煤航测遥感集团有限公司 | 地表三维模型构建方法、装置、设备及存储介质 |
CN113902626A (zh) * | 2021-08-23 | 2022-01-07 | 桂林理工大学 | 超宽幅线阵影像附加约束条件的正射纠正方法 |
CN115423696A (zh) * | 2022-07-29 | 2022-12-02 | 上海海洋大学 | 一种自适应线程参数的遥感正射影像并行生成方法 |
CN115423702A (zh) * | 2022-08-23 | 2022-12-02 | 自然资源部国土卫星遥感应用中心 | 制作大区域星载光学和sar影像dom的方法和系统 |
CN117115621A (zh) * | 2023-10-24 | 2023-11-24 | 中国海洋大学 | 基于改进U-Net网络的卫星云图预测方法 |
Also Published As
Publication number | Publication date |
---|---|
CN107341778B (zh) | 2020-11-10 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN107341778A (zh) | 基于卫星控制点库和dem的sar影像正射纠正方法 | |
Tang et al. | Triple linear-array image geometry model of ZiYuan-3 surveying satellite and its validation | |
CN102506824B (zh) | 一种城市低空无人机系统生成数字正射影像图的方法 | |
KR101165523B1 (ko) | 다중 소스의 지리 정보를 이용한 지리공간 모델링 시스템 및 관련 방법 | |
CN103129752B (zh) | 一种基于地面导航的光学遥感卫星姿态角误差动态补偿方法 | |
CN109977344B (zh) | 一种星载夜光遥感影像的区域网平差方法 | |
CN109801219A (zh) | 面向在线地图叠加的gis数据校正方法及装置 | |
Tang et al. | Verification of ZY-3 satellite imagery geometric accuracy without ground control points | |
CN105387846B (zh) | 卫星影像的正射校正方法和系统 | |
CN108226982B (zh) | 单线阵卫星激光联合高精度定位处理方法 | |
Gullu | Coordinate transformation by radial basis function neural network | |
US20040122633A1 (en) | Method for updating IKONOS RPC data by additional GCP | |
CN102508260A (zh) | 一种面向侧视中分辨率卫星的几何成像构建方法 | |
CN106600551A (zh) | 一种适用于大场景星载sar图像的高精度几何校正方法 | |
CN106443676A (zh) | 稀少控制点星载合成孔径雷达图像对地定位方法 | |
CN107330927A (zh) | 机载可见光图像定位方法 | |
Skaloud et al. | Theory and reality of direct georeferencing in national coordinates | |
CN110030968B (zh) | 一种基于星载立体光学影像的地面遮挡物仰角测量方法 | |
CN105571598B (zh) | 一种卫星激光高度计足印相机姿态的测定方法 | |
Xiao et al. | Prospects for mapping temporal height variations of the seasonal CO2 snow/ice caps at the Martian poles by co-registration of MOLA Profiles | |
CN104599285A (zh) | 一种基于遥感图像的水体信息提取方法及装置 | |
CN110189283B (zh) | 基于语义分割图的遥感图像dsm融合方法 | |
Crespi et al. | DSM generation from very high optical and radar sensors: Problems and potentialities along the road from the 3D geometric modeling to the Surface Model | |
Shaker | Satellite sensor modeling and 3D geo-positioning using empirical models | |
Kim et al. | Estimation and improvement in the geolocation accuracy of rational polynomial coefficients with minimum GCPs using KOMPSAT-3A |
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 | ||
CB02 | Change of applicant information |
Address after: 100048 No. 28 Lianhuachi West Road, Haidian District, Beijing Applicant after: Ministry of Natural Resources Land Satellite Remote Sensing Application Center Address before: 100048 No. 28 Lianhuachi West Road, Haidian District, Beijing Applicant before: Satellite Surveying and Mapping Application Center, NASG |
|
CB02 | Change of applicant information | ||
GR01 | Patent grant | ||
GR01 | Patent grant |