CN111239737A - 一种基于升降轨数据的区域sar影像叠掩补偿方法 - Google Patents
一种基于升降轨数据的区域sar影像叠掩补偿方法 Download PDFInfo
- Publication number
- CN111239737A CN111239737A CN202010200498.1A CN202010200498A CN111239737A CN 111239737 A CN111239737 A CN 111239737A CN 202010200498 A CN202010200498 A CN 202010200498A CN 111239737 A CN111239737 A CN 111239737A
- Authority
- CN
- China
- Prior art keywords
- image
- coordinates
- rpc
- overlap
- main
- 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
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
-
- 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
- G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
- G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
- G01S7/40—Means for monitoring or calibrating
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- Physics & Mathematics (AREA)
- Computer Networks & Wireless Communication (AREA)
- General Physics & Mathematics (AREA)
- Electromagnetism (AREA)
- Image Processing (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明公开了一种基于升降轨数据的区域SAR影像叠掩补偿方法,首先读取SAR区域升轨和降轨影像及辅助文件;根据文件内容构建RPC模型;确定主副影像集;主副影像集中主影像与副影像重叠度计算;主副影像叠掩补偿关系映射集合建立;获取控制点和连接点坐标;构建区域网平差误差方程,并进行逐点求解,获取仿射变换参数;更新RPC模型;主副影像正射纠正;依据RPC定位模型生成叠掩掩膜;纠正后的主影像集依据映射集合进行叠掩补偿。
Description
技术领域
本发明属于微波遥感技术领域,涉及一种SAR影像叠掩补偿方法,尤其涉及一种基于升降轨数据的区域SAR影像叠掩补偿方法。
背景技术
合成孔径雷达(SAR),作为主动式微波遥感技术,具有全天时、全天候成像能力,对长期云、雨覆盖的地方提供有效的观测手段,能为光学遥感提供补充。目前已经应用在测绘,海洋,等方面。然而,SAR由于侧视成像的特点,当地形起伏较大时发生严重几何畸变,如叠掩、阴影等。相比阴影,在山区,叠掩是制约着SAR进一步应用的主要因素。针对叠掩现象,单景影像无法解决,通常做法采用升降轨(或多角度)多张SAR影像,利用信息冗余,将主影像异常像元值用对应地理位置副影像正常像元值进行替代,从而对丢失信息进行补偿。目前关于叠掩补偿研究所采用的影像数量一般不超过三景,而SAR影像的应用已不仅仅局限于单景影像,而更偏向于区域SAR影像。针对区域SAR影像,现有方法可能会导致影像集无法完整覆盖研究区,且一景主影像可能会与多张副影像存在重叠,主副影像补偿关系缺乏定量化指标。因此迫切需要一种以区域SAR影像为对象的叠掩补偿方法。
发明内容
本发明的目的是提供一种基于升降轨数据的区域SAR影像叠掩补偿方法。
本发明所采用的技术方案是:一种基于升降轨数据的区域SAR影像叠掩补偿,包括如下步骤:
步骤1,读取升降轨道星载SAR影像数据,整理RPC参数文件格式符合国际规范;
步骤2,读取RPC定位参数文件,利用影像的RPC参数文件里的模型参数,构建影像的RPC模型;
步骤3,根据影像覆盖关系,进行主副影像集判断,通过判读升轨降轨影像集对研究区的覆盖情况,将覆盖程度较高视向的影像作为主影像集,另一视向的影像作为副影像集;
步骤4,逐个计算主影像与副影像的重叠度,将重叠度作为叠掩补偿的阈值;
步骤5,依据重叠度,建立主副影像叠掩补偿关系映射集合;
步骤6,获取区域影像连接点、控制点坐标,并将各点坐标保存为对应的点文件中;
步骤7,根据步骤6中获取的控制点、连接点坐标,建立区域网平差误差方程,并对误差方程利用最小二乘原理进行求解,最终获取仿射变换参数;
步骤8,利用仿射变换参数对原始RPC模型进行更新;
步骤9,结合对应DEM数据,采用间接校正方法分别对主副影像集进行正射纠正;
步骤10,结合DEM数据、SAR影像、更新后的RPC模型进行叠掩掩膜生成;
步骤11,将纠正后的主影像集,结合叠掩掩膜与主副影像映射集合分别进行叠掩补偿。
进一步的,根据权利要求1所述的基于升降轨数据的区域SAR影像叠掩补偿方法,其特征在于:步骤2中,RPC模型定义为:
为方便表达,将式(1)也写成如下形式:
其中:(Xn,Yn,Zn)为正则化的地面点大地坐标,(rn,cn)为正则化的影像坐标,Fr(Xn,Yn,Zn)、Fc(Xn,Yn,Zn)分别与rn、cn含义相同,P1(Xn,Yn,Zn)、P2(Xn,Yn,Zn)、P3(Xn,Yn,Zn)、P4(Xn,Yn,Zn)记为P1、P2、P3、P4,其中P1、P2、P3、P4是关于Xn、Yn、Zn的三次多项式,没有实际的物理意义;
P1=a0+a1Zn+a2Yn+a3Xn+a4ZnYn+a5ZnXn+a6YnXn+a7Zn2+a8Yn 2+a9Xn 2+a10ZnYnXn+a11Zn 2Yn+a12Zn 2Xn+a13Yn 2Zn+a14Yn 2Xn+a15ZnXn 2+a16YnXn 2+a17Zn 3+a18Yn 3+a19Xn 3
P2=b0+b1Zn+b2Yn+b3Xn+b4ZnYn+b5ZnXn+b6YnXn+b7Zn 2+b8Yn 2+b9Xn 2+b10ZnYnXn+b11Zn 2Yn+b12Zn 2Xn+b13Yn 2Zn+b14Yn 2Xn+b15ZnXn 2+b16YnXn 2+b17Zn 3+b18Yn 3+b19Xn 3
P3=c0+c1Zn+c2Yn+c3Xn+c4ZnYn+c5ZnXn+c6YnXn+c7Zn 2+c8Yn 2+c9Xn 2+c10ZnYnXn+c11Zn 2Yn+c12Zn 2Xn+c13Yn 2Zn+c14Yn 2Xn+c15ZnXn 2+c16YnXn 2+c17Zn 3+c18Yn 3+c19Xn 3
P4=d0+d1Zn+d2Yn+d3Xn+d4ZnYn+d5ZnXn+d6YnXn+d7Zn 2+d8Yn 2+d9Xn 2+d10ZnYnXn+d11Zn 2Yn+d12Zn 2Xn+d13Yn 2Zn+d14Yn 2Xn+d15ZnXn 2+d16YnXn 2+d17Zn 3+d18Yn 3+d19Xn 3
其中,a0,…,a19,b0,…,b19,c0,…,c19,d0,…,d19是RPC参数文件自带模型参数;
正则化的地面点大地坐标(Xn,Yn,Zn)、正则化的影像坐标(rn,cn)与地面点大地坐标(X,Y,Z)和影像坐标(R,C)对应关系如式(3)、(4)所示:
式(3)、(4)中,X为经度、Y为纬度,Z为高程,Xo、Yo、Zo、Ro、Co为正则化的平移参数;Xs、Ys、Zs、Rs、Cs为正则化比例系数,均为RPC文件自带模型参数。
进一步的,步骤4中重叠度计算方法如下,
步骤4.1,获取副影像四角点影像坐标,记为(xm,ym),m=1,2,3,4;
步骤4.2,结合副影像RPC定位文件与给定参考高程面H,采用RPC模型正变换计算(xm,ym),m=1,2,3,4,对应的地面点大地坐标(Bm,Lm,H),m=1,2,3,4;其中,Bm、Lm和H分别代表对应点的纬度、经度和高程;
步骤4.3,将地面点大地坐标(Bm,Lm,H),m=1,2,3,4,利用主影像RPC定位文件,采用RPC模型反变换,计算其在主影像上影像坐标(x′m,y′m),m=1,2,3,4;
步骤4.4,结合主影像四角点影像坐标、步骤4.3计算的副影像四角点投影到主影像上影像坐标(x′m,y′m),m=1,2,3,4,通过建立多边形并判断是否相交从而统计重叠像元个数;
步骤4.5,重叠度计算,公式如下:
其中,Rateoverlap为重叠度,Soverlap为重叠像元个数,Smaster表示主影像像元个数。
进一步的,重叠度Rateoverlap阈值设置20%,即大于20%的情况下副影像才参与主影像进行叠掩补偿。
进一步的,步骤5中主副影像映射集合建立的具体实现方式如下,
步骤5.1,映射关系初步表达:首先将主影像集{M1,m2...,Mp}第k景影像与副影像集{S1,S2,...,Sq}中影像对应关系表示为{S1,S2,...,Sq}Mk,其中,p表示主影像为p景,q表示副影像为q景,且k≤p;引入重叠度后,对应关系重新表示为其中中rate_kq代表主影像k与副影像q的重叠度;
进一步的,骤6中,利用自动匹配技术或者人工判读方式,获取主副影像集中相邻影像重叠区域的同名点,并将获得的同名点作为影像连接点;根据控制点坐标信息,在影像上转刺出对应的像方坐标;其中,控制点布设在研究区四角,以便于控制整个研究区;整理控制点物方、像方坐标和连接点像方信息,并保存成相对应格式的点文件。
进一步的,步骤7中的具体实现方式如下,
步骤7.1,确定仿射变换模型;
使用像方补偿方案对RFM中的系统误差进行补偿时,RPC模型描述的影像坐标(R,C)和地面点大地坐标(X,Y,Z)之间的关系被修正为:
式中,(ΔR,ΔC)为影像坐标(R,C)的系统误差补偿值,采用仿射变换模型:
(e0,e1,e2,f0,f1,f2)为仿射变换参数,即影像定向参数;
步骤7.2,误差方程建立;
将仿射变换系数和目标点的地面点大地坐标作为未知数,得到基于RPC模型的光束法平差误差方程矩阵形式为:
上式进一步表示为:
V=At+Bx-l (8)
式中,V=[vR vC]T为影像行、列坐标观测值残差向量;t=[Δe0 Δe1 Δe2 Δf0Δf1 Δf2]T为影像坐标系统误差补偿参数增量向量;x=[ΔX ΔY ΔZ]T为地面点大地坐标增量向量;
具体形式为:
对于控制点而言,对应的x=[ΔX ΔY ΔZ]T为零向量;
步骤7.3,待定参数求解;
每一个像点均可建立两个误差方程,当量测了足够多的像点时,根据最小二乘平差原理形成法方程:
按照光束法区域网平差中大规模法方程解算策略进行求解,当参数求解完成后,更新地面点大地坐标和仿射变换参数;然后将更新的值参作为下一次迭代计算的初值,直至影像仿射变换参数中的平移参数e0、f0最后两次差值小于阈值时,平差迭代结束;当不满足预设条件时,继续迭代计算,直到满足迭代收敛条件;如果迭代次数达到预设迭代次数,仍然不能收敛,那么平差失败退出,平差完成后输出最终的仿射变换参数。
进一步的,步骤8中的具体实现方式如下,
步骤8.1,建立虚拟格网获取虚拟控制点空间坐标;
(1)根据正则化偏移参数,在高程方向以一定的间隔分层,在平面上,以一定的格网大小建立影像规则格网,生成控制点,即控制点大地坐标(Latitude,Longitude,Height),其中Latitude表示纬度、Longitude表示经度、Height表示高程;
(2)参数正则化
利用公式(3),对获取的控制点大地坐标(Latitude,Longitude,Height)进行正则化处理,获得正则化的控制点大地坐标(Xn,Yn,Zn);
步骤8.2,控制点像方坐标求解;
(1)求解控制点初始影像坐标:将正则化的控制点大地坐标(Xn,Yn,Zn)结合初始RPC参数,利用公式(1)求解对应的初始正则化影像坐标(rn,cn);
(2)初始影像坐标获取:将公式(4)改写成下式
将正则化影像坐标(rn,cn)代入式(10)中,获取初始影像坐标(R,C);
(3)获取仿射变换后的坐标:利用求解的仿射变换参数,对像方误差进行补偿,获得真实的影像坐标(R′,C′),
最终获取仿射变化后的真实影像坐标;
步骤8.3,RPC参数重新求解,
将公式(1)变形为
进而得到误差方程如下:
V=Kr-m,P (13)
式中,r=[au bv cu dv]T,P为权矩阵,一般为单位阵;
根据最小而成原理求解待求变量:
r=(KTK)-1(KTm) (14)
平差的误差方程为线性模型,在求解过程中无需初值;
步骤8.4,将更新的RPC参数重新写入新的定位文件,其中公式(3)和(4)中正则化的平移参数Xo、Yo、Zo、Ro、Co;正则化比例系数Xs、Ys、Zs、Rs、Cs均不变。
进一步的,步骤10中的具体实现方式如下,
步骤10.1,根据原始SAR影像四角点及RPC参数确定所需DEM范围,并对该区域DEM重采样至纠正后的影像分辨率相同大小;
步骤10.2,根据重采样后DEM行列号和像素值获取对应的地面点大地坐标(B,L,H),其中,B、L和H分别代表纬度、经度和高程;随后利用RPC定位模型计算对应SAR影像坐标(x,y),将影像坐标(x,y)对应的地面点大地坐标用(Bg,Lg,Hg)(x,y),g=1,2...,w表示,其中g的取值为1-w,表示同一个影像坐标(x,y)对应地面点大地坐标不超过w个;
步骤10.3,遍历SAR影像坐标,判断(x,y)对应的地面点大地坐标(Bg,Lg,Hg)(x,y)的个数,当g>2,则为叠掩区;此时将地面点大地坐标(Bg,Lg,Hg)(x,y)对应位置标记为flag=2,否则设置为正常区域记为flag=1,生成原始掩膜图像;
步骤10.4,对步骤10.3中原始掩膜图进行编辑,采用形态学闭运算减少掩膜图中的“孔洞”效应。
进一步的,步骤11中的具体实现方式如下,
其中,∈表示像元属于图像范畴;Λ和→分别表示“且”和“蕴含”运算;M(i,j),AL和AN分别代表主影像上任意像元、叠掩区域范围和正常区域范围;S(i′,j′)、BL和BN分别代表副影像上对应主影像M(i,j)位置像元、副影像上叠掩区域范围和副影像上正常区域范围;R(i,j)代表补偿后影像上像元值;
步骤11.2,当主影像上叠掩像元被补偿时,该像元对应的掩膜位置标记为正常,补偿前叠掩像元标记为flag=2,补偿后该位置对应标记为flag=1,遍历整张影像更新主影像掩膜;
步骤11.4,对研究区主影像集{M1,M2...,Mp}每景影像执行步骤11.1到步骤11.3,直到所有主影像都完成叠掩补偿为止。
本发明产生的有益效果是:
(1)针对研究区多景升降轨影像数据,采用升轨影像或降轨影像对研究区的覆盖程度高优先选为主影像集,从而确保主影像完整覆盖研究区;
(2)采用影像间重叠度构建主副影像映射集合,为叠掩补偿提供定量化指标,提出的整个流程,易于自动化处理;
(3)采用RPC模型替代RD模型,避免了针对不同卫星平台的分别参数设置,实现了各类影像类型的统一处理;
(4)采用结合控制点的区域网平差方法,保证影像间的相对和绝对定位精度,同时最终的高精度RPC文件用于生成叠掩掩膜,有效利用辅助文件,提升文件利用率。
附图说明
图1是本发明实施例的方法流程图。
具体实施方式
为了便于本领域普通技术人员理解和实施本发明,附图1为方法流程图,下面对本发明作进一步的详细描述,应当理解,此处所描述的实施示例仅用于说明和解释本发明,并不用于限定本发明。
请见图1,本发明提供的一种基于升降轨数据的区域SAR影像叠掩补偿方法,包括以下步骤:
步骤1:读取升降轨道星载SAR影像,辅助RPC定位文件,并检查RPC文件格式符合国际规范,若不符合规范要求,整理成国际规范格式。
步骤2:读取步骤1中的RPC定位文件,利用影像的RPC参数文件里的模型参数,构建影像的RPC模型;
RPC模型定义为:
为方便表达,将式(1)也可写成如下形式:
其中:(Xn,Yn,Zn)为正则化的地面点大地坐标,(rn,cn)为正则化的影像坐标,Fr(Xn,Yn,Zn)、Fc(Xn,Yn,Zn)分别与rn、cn含义相同,为方便后文求偏导,用函数方式进行表达;P1(Xn,Yn,Zn)、P2(Xn,Yn,Zn)、P3(Xn,Yn,Zn)、P4(Xn,Yn,Zn)记为P1、P2、P3、P4,P1、P2、P3、P4是关于Xn、Yn、Zn的三次多项式,没有实际的物理意义。
P1=a0+a1Zn+a2Yn+a3Xn+a4ZnYn+a5ZnXn+a6YnXn+a7Zn 2+a8Yn 2+a9Xn 2+a10ZnYnXn+a11Zn 2Yn+a12Zn 2Xn+a13Yn 2Zn+a14Yn 2Xn+a15ZnXn 2+a16YnXn 2+a17Zn 3+a18Yn 3+a19Xn 3
P2=b0+b1Zn+b2Yn+b3Xn+b4ZnYn+b5ZnXn+b6YnXn+b7Zn 2+b8Yn 2+b9Xn 2+b10ZnYnXn+b11Zn 2Yn+b12Zn 2Xn+b13Yn 2Zn+b14Yn 2Xn+b15ZnXn 2+b16YnXn 2+b17Zn 3+b18Yn 3+b19Xn 3
P3=c0+c1Zn+c2Yn+c3Xn+c4ZnYn+c5ZnXn+c6YnXn+c7Zn 2+c8Yn 2+c9Xn 2+c10ZnYnXn+c11Zn 2Yn+c12Zn 2Xn+c13Yn 2Zn+c14Yn 2Xn+c15ZnXn 2+c16YnXn 2+c17Zn 3+c18Yn 3+c19Xn 3
P4=d0+d1Zn+d2Yn+d3Xn+d4ZnYn+d5ZnXn+d6YnXn+d7Zn 2+d8Yn 2+d9Xn 2+d10ZnYnXn+d11Zn 2Yn+d12Zn 2Xn+d13Yn 2Zn+d14Yn 2Xn+d15ZnXn 2+d16YnXn 2+d17Zn 3+d18Yn 3+d19Xn 3
其中,a0,…,a19,b0,…,b19,c0,…,c19,d0,…,d19是RPC参数文件自带模型参数;b0与d0通常为1。
正则化的地面点大地坐标(Xn,Yn,Zn)、正则化的影像坐标(rn,cn)与地面点大地坐标(X,Y,Z)和影像坐标(R,C)对应关系如式(3)、(4)所示:
式(3)、(4)中,X为经度、Y为纬度,Z为高程,Xo、Yo、Zo、Ro、Co为正则化的平移参数;Xs、Ys、Zs、Rs、Cs为正则化比例系数,均为RPC文件自带模型参数。
步骤3:根据影像覆盖关系,进行主副影像集判断。首先加载研究区的矢量边界,并利用影像自带RPC定位文件,将升降轨影像可视化展现在专业软件(如ArcGIS,ENVI等);人工判读升轨降轨影像集对研究区的覆盖情况;将覆盖程度较高视向的影像作为主影像集,另一视向的影像作为副影像集。
步骤4:逐个计算主影像与副影像的重叠度,将重叠度作为叠掩补偿的阈值;具体实现包括如下步骤:
步骤4.1:获取副影像四角点影像坐标,记为(xm,ym),m=1,2,3,4.
步骤4.2:结合副影像RPC定位文件与给定参考高程面H,采用RPC模型正变换计算(xm,ym),m=1,2,3,4.对应的地面点大地坐标(Bm,Lm,H),m=1,2,3,4.。其中,Bm、Lm和H分别代表对应点的纬度、经度和高程。
步骤4.3:将地面点大地坐标(Bm,Lm,H),m=1,2,3,4.,利用主影像RPC定位文件,采用RPC模型反变换,计算其在主影像上影像坐标(x′m,y′m),m=1,2,3,4.
步骤4.4:结合主影像四角点影像坐标、步骤4.3计算的副影像四角点投影到主影像上影像坐标(x′m,y′m),m=1,2,3,4.,通过建立多边形并判断是否相交从而统计重叠像元个数。
步骤4.5,重叠度计算,公式如下:
其中,Rateoverlap为重叠度,Soverlap为重叠像元个数,Smaster表示主影像像元个数,可由主影像宽和高相乘计算得出。
在本实例中,将重叠度Rateoverlap阈值设置20%,即大于20%的情况下副影像才参与主影像进行叠掩补偿,低于20%的情况由于对主影像的补偿部分十分有限,这里不予考虑,但不限于此。
步骤5:依据重叠度,建立主副影像叠掩补偿关系映射集合;
步骤5.1:映射关系初步表达:首先将主影像集{M1,M2...,Mp}第k景影像与副影像集{S1,S2,...,Sq}中影像对应关系表示为其中,p表示主影像为p景,q表示副影像为q景,且k≤p;引入重叠度后,对应关系重新表示为其中中rate_kq代表主影像k与副影像q的重叠度。
步骤5.2:重叠度阈值筛选:将中满足rate_ki<20%,的影像删除,其中,i表示为副影像中任意一景影像且i≤q(这里20%为步骤4设置的阈值)。记为其中t≤q,表示参与主影像k对应的副影像有t景。
步骤6:利用自动匹配技术或者人工判读方式,获取主副影像集中相邻影像重叠区域的同名点,并将获得的同名点作为影像连接点;根据控制点坐标信息,在影像上转刺出对应的像方坐标。其中,控制点尽量布设在研究区四角,以便于控制整个研究区;整理控制点物方、像方坐标和连接点像方信息,并保存成相对应格式的点文件。
步骤7:根据步骤6中获取的控制点、连接点坐标,建立区域网平差误差方程,并对误差方程利用最小二乘原理进行求解。
步骤7.1:确定仿射变换模型
使用像方补偿方案对PRC模型中的系统误差进行补偿时,RPC模型描述的影像坐标(R,C)和地面点大地坐标(X,Y,Z)之间的关系被修正为:
式中,(ΔR,ΔC)为影像坐标(R,C)的系统误差补偿值,一般采用仿射变换模型:
(e0,e1,e2,f0,f1,f2)为仿射变换参数,即影像定向参数。
步骤7.2:误差方程建立
将仿射变换系数和地面点的大地坐标作为未知数,则可得到基于RPC模型平差误差方程矩阵形式为:
上式可以表示为:
V=At+Bx-l (8)
式中,V=[vR vC]T为影像行、列坐标观测值残差向量;t=[Δe0 Δe1 Δe2 Δf0Δf1 Δf2]T为影像坐标系统误差补偿参数增量向量;x=[ΔX ΔY ΔZ]T为地面点大地坐标增量向量;
具体形式为:
对于控制点而言,对应的x=[ΔX ΔY ΔZ]T为零向量。
步骤7.3,待定参数求解
每一个像点均可建立两个误差方程,当量测了足够多的像点时,根据最小二乘平差原理形成法方程:
按照光束法区域网平差中大规模法方程解算策略进行求解,当参数求解完成后,更新地面点大地坐标和仿射变换参数;然后将更新的值作为下一次迭代计算的初值,直至影像仿射变换参数中的平移参数e0、f0最后两次差值小于阈值时(本实施例为0.1个像元pixel,但是不限于此)时,平差迭代结束;当不满足预设条件时,继续迭代计算,直到满足迭代收敛条件;如果迭代次数达到预设迭代次数(本实施例设定的迭代次数为30次,但是不限于此),仍然不能收敛,那么平差失败退出。平差完成后输出最终的仿射变换参数。
步骤8:利用定向参数(仿射变换参数)对原始RPC模型进行更新;
步骤8.1:建立虚拟格网获取虚拟控制点空间坐标
(1)根据正则化偏移参数,在高程方向以一定的间隔分层,在平面上,以一定的格网大小建立影像规则格网,生成控制点,即控制点大地坐标(Latitude,Longitude,Height),其中Latitude表示纬度、Longitude表示经度、Height表示高程;
(2)参数正则化
利用公式(3),对获取的控制点大地坐标(Latitude,Longitude,Height)进行正则化处理,获得正则化的控制点大地坐标(Xn,Yn,Zn)。
步骤8.2控制点像方坐标求解
(1)求解控制点初始影像坐标:将正则化的控制点大地坐标(Xn,Yn,Zn结合初始RPC参数,利用公式(1)求解对应的初始正则化影像坐标(rn,cn)。
(2)初始影像坐标获取:公式(4)可以改写成下式
将正则化影像坐标(rn,cn)代入式(10)中,获取初始影像坐标(R,C)。
(3)获取仿射变换后的坐标
利用求解的仿射变换参数,对像方误差进行补偿,获得真实的影像坐标(R′,C′)
最终获取仿射变化后的真实影像坐标
步骤8.3,RPC参数重新求解
将公式(1)变形为
进而得到误差方程如下:
V=Kr-m,P (13)
式中,r=[au bv cu dv]T,P为权矩阵,一般为单位阵;
根据最小而成原理求解待求变量:
r=(KTK)-1(KTm) (14)
平差的误差方程为线性模型,在求解过程中无需初值。
步骤8.4将更新的RPC参数重新写入新的定位文件,其中公式(3)和(4)中正则化的平移参数Xo、Yo、Zo、Ro、Co;正则化比例系数Xs、Ys、Zs、Rs、Cs均不变。
步骤9:结合对应DEM数据,采用间接校正方法分别对主副影像集进行正射纠正。首先在DEM上根据空间坐标计算影像上对应的像点坐标;然后根据像点坐标在影像上对灰度值进行重采样;将重采样的灰度值赋予DEM对应位置;遍历影像对应的DEM,执行相同操作,完成正射纠正。
步骤10:结合DEM数据、SAR影像、更新后的RPC模型进行叠掩掩膜生成;
采用RPC定位模型,本实例叠掩区域判定方法步骤如下:
步骤10.1:根据原始SAR影像四角点及RPC参数确定所需DEM范围,并对该区域DEM重采样至纠正后的影像分辨率相同大小;
步骤10.2:根据重采样后DEM行列号和像素值获取对应的地面点大地坐标(B,L,H),其中,B、L和H分别代表纬度、经度和高程;随后利用RPC定位模型计算对应SAR影像坐标(x,y),将影像坐标(x,y)对应的地面点大地坐标用(Bg,Lg,Hg)(x,y),g=1,2...,w.表示,其中g的取值为1-w,表示同一个影像坐标(x,y)对应地面点大地坐标不超过w个。
步骤10.3:遍历SAR影像坐标,判断(x,y)对应的地面点大地坐标(Bg,Lg,Hg)(x,y)的个数,当g>2,则为叠掩区。此时将地面点大地坐标(Bg,Lg,Hg)(x,y)对应位置标记为flag=2,否则设置为正常区域记为flag=1,生成原始掩膜图像;
步骤10.4:对步骤10.3中原始掩膜图进行编辑,采用形态学闭运算减少掩膜图中的“孔洞”效应。
步骤11:将纠正后的主影像集,结合叠掩掩膜与主副影像映射集合分别进行叠掩补偿;
其中,∈表示像元属于图像范畴;Λ和→分别表示“且”和“蕴含”运算;M(i,j)、AL和AN分别代表主影像上任意像元、叠掩区域范围和正常区域范围;S(i′,j′)、BL和BN分别代表副影像上对应主影像M(i,j)位置像元、副影像上叠掩区域范围和副影像上正常区域范围;R(i,j)代表补偿后影像上像元值。这里的正常区域范围表示通过叠掩阴影掩膜对图像中的叠掩进行掩膜后剩下的区域。
步骤11.2:当主影像上叠掩像元被补偿时,该像元对应的掩膜位置标记为正常。在步骤10.3中提到,补偿前叠掩像元标记为flag=2,补偿后该位置对应标记为flag=1,遍历整张影像更新主影像掩膜。
步骤11.4:对研究区主影像集{M1,M2…,Mp}每景影像执行步骤11.1到步骤11.3,直到所有主影像都完成叠掩补偿为止。
应当理解的是,本说明书未详细阐述的部分均属于现有技术。本领域的普通技术人员在本发明的启示下,在不脱离本发明权利要求所保护的范围情况下,还可以做出替换或变形,均落入本发明的保护范围之内,本发明的请求保护范围应以所附权利要求为准。
Claims (10)
1.一种基于升降轨数据的区域SAR影像叠掩补偿方法,其特征在于,包括以下步骤:
步骤1,读取升降轨道星载SAR影像数据,整理RPC参数文件格式符合国际规范;
步骤2,读取RPC定位参数文件,利用影像的RPC参数文件里的模型参数,构建影像的RPC模型;
步骤3,根据影像覆盖关系,进行主副影像集判断,通过判读升轨降轨影像集对研究区的覆盖情况,将覆盖程度较高视向的影像作为主影像集,另一视向的影像作为副影像集;
步骤4,逐个计算主影像与副影像的重叠度,将重叠度作为叠掩补偿的阈值;
步骤5,依据重叠度,建立主副影像叠掩补偿关系映射集合;
步骤6,获取区域影像连接点、控制点坐标,并将各点坐标保存为对应的点文件中;
步骤7,根据步骤6中获取的控制点、连接点坐标,建立区域网平差误差方程,并对误差方程利用最小二乘原理进行求解,最终获取仿射变换参数;
步骤8,利用仿射变换参数对原始RPC模型进行更新;
步骤9,结合对应DEM数据,采用间接校正方法分别对主副影像集进行正射纠正;
步骤10,结合DEM数据、SAR影像、更新后的RPC模型进行叠掩掩膜生成;
步骤11,将纠正后的主影像集,结合叠掩掩膜与主副影像映射集合分别进行叠掩补偿。
2.根据权利要求1所述的基于升降轨数据的区域SAR影像叠掩补偿方法,其特征在于:步骤2中,RPC模型定义为:
为方便表达,将式(1)也写成如下形式:
其中:(Xn,Yn,Zn)为正则化的地面点大地坐标,(rn,cn)为正则化的影像坐标,Fr(Xn,Yn,Zn)、Fc(Xn,Yn,Zn)分别与rn、cn含义相同,P1(Xn,Yn,Zn)、P2(Xn,Yn,Zn)、P3(Xn,Yn,Zn)、P4(Xn,Yn,Zn)记为P1、P2、P3、P4,其中P1、P2、P3、P4是关于Xn、Yn、Zn的三次多项式,没有实际的物理意义;
P1=a0+a1Zn+a2Yn+a3Xn+a4ZnYn+a5ZnXn+a6YnXn+a7Zn 2+a8Yn 2
+a9Xn 2+a10ZnYnXn+a11Zn 2Yn+a12Zn 2Xn+a13Yn 2Zn
+a14Yn 2Xn+a15ZnXn 2+a16YnXn 2+a17Zn 3+a18Yn 3+a19Xn 3
P2=b0+b1Zn+b2Yn+b3Xn+b4ZnYn+b5ZnXn+b6YnXn+b7Zn 2+b8Yn 2
+b9Xn 2+b10ZnYnXn+b11Zn 2Yn+b12Zn 2Xn+b13Yn 2Zn
+b14Yn 2Xn+b15ZnXn 2+b16YnXn 2+b17Zn 3+b18Yn 3+b19Xn 3
P3=c0+c1Zn+c2Yn+c3Xn+c4ZnYn+c5ZnXn+c6YnXn+c7Zn 2+c8Yn 2
+c9Xn 2+c10ZnYnXn+c11Zn 2Yn+c12Zn 2Xn+c13Yn 2Zn
+c14Yn 2Xn+c15ZnXn 2+c16YnXn 2+c17Zn 3+c18Yn 3+c19Xn 3
P4=d0+d1Zn+d2Yn+d3Xn+d4ZnYn+d5ZnXn+d6YnXn+d7Zn 2+d8Yn 2
+d9Xn 2+d10ZnYnXn+d11Zn 2Yn+d12Zn 2Xn+d13Yn 2Zn
+d14Yn 2Xn+d15ZnXn 2+d16YnXn 2+d17Zn 3+d18Yn 3+d19Xn 3
其中,a0,…,a19,b0,…,b19,c0,…,c19,d0,…,d19是RPC参数文件自带模型参数;
正则化的地面点大地坐标(Xn,Yn,Zn)、正则化的影像坐标(rn,cn)与地面点大地坐标(X,Y,Z)和影像坐标(R,C)对应关系如式(3)、(4)所示:
式(3)、(4)中,X为经度、Y为纬度,Z为高程,Xo、Yo、Zo、Ro、Co为正则化的平移参数;Xs、Ys、Zs、Rs、Cs为正则化比例系数,均为RPC文件自带模型参数。
3.根据权利要求1所述的基于升降轨数据的区域SAR影像叠掩补偿方法,其特征在于:步骤4中重叠度计算方法如下,
步骤4.1,获取副影像四角点影像坐标,记为(xm,ym),m=1,2,3,4;
步骤4.2,结合副影像RPC定位文件与给定参考高程面H,采用RPC模型正变换计算(xm,ym),m=1,2,3,4,对应的地面点大地坐标(Bm,Lm,H),m=1,2,3,4;其中,Bm、Lm和H分别代表对应点的纬度、经度和高程;
步骤4.3,将地面点大地坐标(Bm,Lm,H),m=1,2,3,4,利用主影像RPC定位文件,采用RPC模型反变换,计算其在主影像上影像坐标(x′m,y′m),m=1,2,3,4;
步骤4.4,结合主影像四角点影像坐标、步骤4.3计算的副影像四角点投影到主影像上影像坐标(x′m,y′m),m=1,2,3,4,通过建立多边形并判断是否相交从而统计重叠像元个数;
步骤4.5,重叠度计算,公式如下:
其中,Rateoverlap为重叠度,Soverlap为重叠像元个数,Smaster表示主影像像元个数。
4.根据权利要求3所述的基于升降轨数据的区域SAR影像叠掩补偿方法,其特征在于:重叠度Rateoverlap阈值设置20%,即大于20%的情况下副影像才参与主影像进行叠掩补偿。
5.根据权利要求1所述的基于升降轨数据的区域SAR影像叠掩补偿方法,其特征在于:步骤5中主副影像映射集合建立的具体实现方式如下,
步骤5.1,映射关系初步表达:首先将主影像集{M1,M2…,Mv}第k景影像与副影像集{S1,S2,...,Sq}中影像对应关系表示为其中,p表示主影像为p景,q表示副影像为q景,且k≤p;引入重叠度后,对应关系重新表示为其中中rate_kq代表主影像k与副影像q的重叠度;
6.根据权利要求1所述的基于升降轨数据的区域SAR影像叠掩补偿方法,其特征在于:步骤6中,利用自动匹配技术或者人工判读方式,获取主副影像集中相邻影像重叠区域的同名点,并将获得的同名点作为影像连接点;根据控制点坐标信息,在影像上转刺出对应的像方坐标;其中,控制点布设在研究区四角,以便于控制整个研究区;整理控制点物方、像方坐标和连接点像方信息,并保存成相对应格式的点文件。
7.根据权利要求2所述的基于升降轨数据的区域SAR影像叠掩补偿方法,其特征在于:步骤7中的具体实现方式如下,
步骤7.1,确定仿射变换模型;
使用像方补偿方案对RFM中的系统误差进行补偿时,RPC模型描述的影像坐标(R,C)和地面点大地坐标(X,Y,Z)之间的关系被修正为:
式中,(ΔR,ΔC)为影像坐标(R,C)的系统误差补偿值,采用仿射变换模型:
(e0,e1,e2,f0,f1,f2)为仿射变换参数,即影像定向参数;
步骤7.2,误差方程建立;
将仿射变换系数和目标点的地面点大地坐标作为未知数,得到基于RPC模型的光束法平差误差方程矩阵形式为:
上式进一步表示为:
V=At+Bx-l (8)
式中,V=[vR vC]T为影像行、列坐标观测值残差向量;t=[Δe0 Δe1 Δe2 Δf0 Δf1Δf2]T为影像坐标系统误差补偿参数增量向量;x=[ΔX ΔY ΔZ]T为地面点大地坐标增量向量;
具体形式为:
对于控制点而言,对应的x=[ΔX ΔY ΔZ]T为零向量;
步骤7.3,待定参数求解;
每一个像点均可建立两个误差方程,当量测了足够多的像点时,根据最小二乘平差原理形成法方程:
按照光束法区域网平差中大规模法方程解算策略进行求解,当参数求解完成后,更新地面点大地坐标和仿射变换参数;然后将更新的值参作为下一次迭代计算的初值,直至影像仿射变换参数中的平移参数e0、f0最后两次差值小于阈值时,平差迭代结束;当不满足预设条件时,继续迭代计算,直到满足迭代收敛条件;如果迭代次数达到预设迭代次数,仍然不能收敛,那么平差失败退出,平差完成后输出最终的仿射变换参数。
8.根据权利要求2所述的基于升降轨数据的区域SAR影像叠掩补偿方法,其特征在于:步骤8中的具体实现方式如下,
步骤8.1,建立虚拟格网获取虚拟控制点空间坐标;
(1)根据正则化偏移参数,在高程方向以一定的间隔分层,在平面上,以一定的格网大小建立影像规则格网,生成控制点,即控制点大地坐标(Latitude,Longitude,Height),其中Latitude表示纬度、Longitude表示经度、Height表示高程;
(2)参数正则化
利用公式(3),对获取的控制点大地坐标(Latitude,Longitude,Height)进行正则化处理,获得正则化的控制点大地坐标(Xn,Yn,Zn);
步骤8.2,控制点像方坐标求解;
(1)求解控制点初始影像坐标:将正则化的控制点大地坐标(Xn,Yn,Zn)结合初始RPC参数,利用公式(1)求解对应的初始正则化影像坐标(rn,cn);
(2)初始影像坐标获取:将公式(4)改写成下式
将正则化影像坐标(rn,cn)代入式(10)中,获取初始影像坐标(R,C);
(3)获取仿射变换后的坐标:利用求解的仿射变换参数,对像方误差进行补偿,获得真实的影像坐标(R′,C′),
最终获取仿射变化后的真实影像坐标;
步骤8.3,RPC参数重新求解,
将公式(1)变形为
进而得到误差方程如下:
V=Kr-m,P (13)
根据最小而成原理求解待求变量:
r=(KTK)-1(KTm) (14)
平差的误差方程为线性模型,在求解过程中无需初值;
步骤8.4,将更新的RPC参数重新写入新的定位文件,其中公式(3)和(4)中正则化的平移参数Xo、Yo、Zo、Ro、Co;正则化比例系数Xs、Ys、Zs、Rs、Cs均不变。
9.根据权利要求1所述的基于升降轨数据的区域SAR影像叠掩补偿方法,其特征在于:步骤10中的具体实现方式如下,
步骤10.1,根据原始SAR影像四角点及RPC参数确定所需DEM范围,并对该区域DEM重采样至纠正后的影像分辨率相同大小;
步骤10.2,根据重采样后DEM行列号和像素值获取对应的地面点大地坐标(B,L,H),其中,B、L和H分别代表纬度、经度和高程;随后利用RPC定位模型计算对应SAR影像坐标(x,y),将影像坐标(x,y)对应的地面点大地坐标用(Bg,Lg,Hg)(x,y),g=1,2...,w表示,其中g的取值为1至w,表示同一个影像坐标(x,y)对应地面点大地坐标不超过w个;
步骤10.3,遍历SAR影像坐标,判断(x,y)对应的地面点大地坐标(Bg,Lg,Hg)(x,y)的个数,当g>2,则为叠掩区;此时将地面点大地坐标(Bg,Lg,Hg)(x,y)对应位置标记为flag=2,否则设置为正常区域记为flag=1,生成原始掩膜图像;
步骤10.4,对步骤10.3中原始掩膜图进行编辑,采用形态学闭运算减少掩膜图中的“孔洞”效应。
10.根据权利要求5所述的基于升降轨数据的区域SAR影像叠掩补偿方法,其特征在于:步骤11中的具体实现方式如下,
其中,∈表示像元属于图像范畴;∧和→分别表示“且”和“蕴含”运算;M(i,j),AL和AN分别代表主影像上任意像元、叠掩区域范围和正常区域范围;S(i′,j′)、BL和BN分别代表副影像上对应主影像M(i,j)位置像元、副影像上叠掩区域范围和副影像上正常区域范围;R(i,j)代表补偿后影像上像元值;
步骤11.2,当主影像上叠掩像元被补偿时,该像元对应的掩膜位置标记为正常,补偿前叠掩像元标记为flag=2,补偿后该位置对应标记为flag=1,遍历整张影像更新主影像掩膜;
步骤11.4,对研究区主影像集{M1,M2...,Mp}中每景影像执行步骤11.1到步骤11.3,直到所有主影像都完成叠掩补偿为止。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010200498.1A CN111239737B (zh) | 2020-03-20 | 2020-03-20 | 一种基于升降轨数据的区域sar影像叠掩补偿方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010200498.1A CN111239737B (zh) | 2020-03-20 | 2020-03-20 | 一种基于升降轨数据的区域sar影像叠掩补偿方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN111239737A true CN111239737A (zh) | 2020-06-05 |
CN111239737B CN111239737B (zh) | 2021-12-03 |
Family
ID=70864515
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202010200498.1A Active CN111239737B (zh) | 2020-03-20 | 2020-03-20 | 一种基于升降轨数据的区域sar影像叠掩补偿方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN111239737B (zh) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112082576A (zh) * | 2020-09-09 | 2020-12-15 | 桂林理工大学 | 一种“三无情况下”的卫星影像正射纠正方法 |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101739677A (zh) * | 2009-12-17 | 2010-06-16 | 中国测绘科学研究院 | Sar正射影像图升降轨融合方法 |
CN102628942A (zh) * | 2012-04-24 | 2012-08-08 | 中国科学院遥感应用研究所 | 一种雷达影像双视向信息补偿方法 |
CN106949880A (zh) * | 2017-03-10 | 2017-07-14 | 中国电建集团昆明勘测设计研究院有限公司 | 海拔起伏较大测区无人机影像局部重叠度过高处理的方法 |
CN109709551A (zh) * | 2019-01-18 | 2019-05-03 | 武汉大学 | 一种星载合成孔径雷达影像的区域网平面平差方法 |
CN109977344A (zh) * | 2019-03-20 | 2019-07-05 | 武汉大学 | 一种星载夜光遥感影像的区域网平差方法 |
-
2020
- 2020-03-20 CN CN202010200498.1A patent/CN111239737B/zh active Active
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101739677A (zh) * | 2009-12-17 | 2010-06-16 | 中国测绘科学研究院 | Sar正射影像图升降轨融合方法 |
CN102628942A (zh) * | 2012-04-24 | 2012-08-08 | 中国科学院遥感应用研究所 | 一种雷达影像双视向信息补偿方法 |
CN106949880A (zh) * | 2017-03-10 | 2017-07-14 | 中国电建集团昆明勘测设计研究院有限公司 | 海拔起伏较大测区无人机影像局部重叠度过高处理的方法 |
CN109709551A (zh) * | 2019-01-18 | 2019-05-03 | 武汉大学 | 一种星载合成孔径雷达影像的区域网平面平差方法 |
CN109977344A (zh) * | 2019-03-20 | 2019-07-05 | 武汉大学 | 一种星载夜光遥感影像的区域网平差方法 |
Non-Patent Citations (1)
Title |
---|
程前等: "基于RFM模型的叠掩区域定位方法", 《航天返回与遥感》 * |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112082576A (zh) * | 2020-09-09 | 2020-12-15 | 桂林理工大学 | 一种“三无情况下”的卫星影像正射纠正方法 |
Also Published As
Publication number | Publication date |
---|---|
CN111239737B (zh) | 2021-12-03 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN101604018B (zh) | 高分辨率遥感影像数据处理方法及其系统 | |
Li et al. | Rigorous photogrammetric processing of HiRISE stereo imagery for Mars topographic mapping | |
CN108305237A (zh) | 考虑不同光照成像条件的多立体影像融合制图方法 | |
CN107144293A (zh) | 一种视频卫星面阵相机的几何定标方法 | |
CN109709551B (zh) | 一种星载合成孔径雷达影像的区域网平面平差方法 | |
CN109696182A (zh) | 一种星载推扫式光学传感器内方位元素定标方法 | |
Kornus et al. | DEM generation from SPOT-5 3-fold along track stereoscopic imagery using autocalibration | |
Albertz et al. | HRSC on Mars Express–Photogrammetric and cartographic research | |
CN114241125B (zh) | 一种基于多视角卫星影像精细三维建模方法及系统 | |
CN110853140A (zh) | 一种dem辅助的光学视频卫星影像稳像方法 | |
CN112465849B (zh) | 一种无人机激光点云与序列影像的配准方法 | |
CN113514829A (zh) | 面向InSAR的初始DSM的区域网平差方法 | |
CN111239737B (zh) | 一种基于升降轨数据的区域sar影像叠掩补偿方法 | |
CN115471619A (zh) | 基于立体成像高分辨率卫星影像的城市三维模型构建方法 | |
CN117576343B (zh) | 基于高分辨率卫星立体影像的三维mesh模型制作方法 | |
CN109029379B (zh) | 一种高精度小基高比立体测绘方法 | |
CN113324527B (zh) | 一种同轨激光测高点与三线阵立体影像联合测绘处理方法 | |
CN114913297A (zh) | 一种基于mvs稠密点云的场景正射影像生成方法 | |
Jacobsen | Calibration and validation of Corona KH-4B to generate height models and orthoimages | |
Zhang et al. | Tests and performance evaluation of DMC images and new methods for their processing | |
CN114777745A (zh) | 一种基于无迹卡尔曼滤波的倾斜取证建模方法 | |
CN113593026A (zh) | 车道线标注辅助地图生成方法、装置和计算机设备 | |
Cao et al. | Precise sensor orientation of high-resolution satellite imagery with the strip constraint | |
Guo et al. | Verification and Analysis of Overseas 1: 50000 Scale Surveying and Mapping 4D Products Quality Generated by the ZY3 and TH1 Satellite Imagery | |
Elshehaby et al. | Assessment of cartographic potential of EgyptSat-1 satellite image (case study in flat areas) |
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 |