CN109116321A - 一种星载干涉成像高度计的相位滤波方法及高度测量方法 - Google Patents
一种星载干涉成像高度计的相位滤波方法及高度测量方法 Download PDFInfo
- Publication number
- CN109116321A CN109116321A CN201810777772.4A CN201810777772A CN109116321A CN 109116321 A CN109116321 A CN 109116321A CN 201810777772 A CN201810777772 A CN 201810777772A CN 109116321 A CN109116321 A CN 109116321A
- Authority
- CN
- China
- Prior art keywords
- pixel
- filtering
- phase
- azimuth
- image
- 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
- 238000001914 filtration Methods 0.000 title claims abstract description 91
- 238000000034 method Methods 0.000 title claims abstract description 64
- 238000003384 imaging method Methods 0.000 title claims abstract description 60
- 238000000691 measurement method Methods 0.000 title abstract description 4
- 238000012545 processing Methods 0.000 claims abstract description 13
- 238000004364 calculation method Methods 0.000 claims description 15
- 238000005259 measurement Methods 0.000 claims description 13
- BTCSSZJGUNDROE-UHFFFAOYSA-N gamma-aminobutyric acid Chemical compound NCCCC(O)=O BTCSSZJGUNDROE-UHFFFAOYSA-N 0.000 claims description 6
- 238000012935 Averaging Methods 0.000 claims description 5
- 230000008569 process Effects 0.000 claims description 3
- 238000009826 distribution Methods 0.000 description 8
- 238000010586 diagram Methods 0.000 description 7
- 230000008859 change Effects 0.000 description 5
- 230000006872 improvement Effects 0.000 description 5
- 238000013507 mapping Methods 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 3
- 238000010587 phase diagram Methods 0.000 description 3
- 230000009897 systematic effect Effects 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 230000000875 corresponding effect Effects 0.000 description 2
- 238000013461 design Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 2
- 101100420946 Caenorhabditis elegans sea-2 gene Proteins 0.000 description 1
- 241000321538 Candidia Species 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 201000003984 candidiasis Diseases 0.000 description 1
- 230000002596 correlated effect Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 230000009977 dual effect Effects 0.000 description 1
- 238000002592 echocardiography Methods 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- PCHJSUWPFVWCPO-UHFFFAOYSA-N gold Chemical compound [Au] PCHJSUWPFVWCPO-UHFFFAOYSA-N 0.000 description 1
- 239000010931 gold Substances 0.000 description 1
- 229910052737 gold Inorganic materials 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 230000002452 interceptive effect Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000009827 uniform distribution Methods 0.000 description 1
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
- 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/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
-
- 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/882—Radar or analogous systems specially adapted for specific applications for altimeters
-
- 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/9021—SAR image post-processing techniques
- G01S13/9023—SAR image post-processing techniques combined with interferometric techniques
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)
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明公开了一种星载干涉成像高度计的相位滤波方法及高度测量方法，所述相位滤波方法包括：对星载干涉成像高度计两个接收通道的数据进行合成孔径成像处理，得到两幅复图像，以其中一幅图像作为基准对两幅复图像进行配准，再计算配准后两幅图像的相干系数；依次计算图像中每个像元的几何位置、观测角和多视视数；由此确定每个像元的方位向滤波窗口大小和距离向滤波窗口大小，并确定每个像元的方位向滤波区间和距离向滤波区间；将两个复图像进行共轭相乘，同时去除地平相位，然后在方位向滤波区间和距离向滤波区间对每个像元分别进行滤波，最终得到滤波后的干涉相位。基于上述相位滤波方法，本发明还提出了一种星载干涉成像高度计的高度测量方法。
Description
技术领域
本发明涉及雷达信号处理技术领域，特别涉及一种星载干涉成像高度计的相位滤波方法及高度测量方法。
背景技术
星载干涉成像高度计是一种利用短基线、小角度干涉技术测量海面和陆地高度的设备，星载干涉成像高度计通过偏离天顶点观测地球表面，通过两个接收天线同时接收地球表面反射的雷达回波，经过处理得到两天线观测图像的相位差，又称为干涉相位，结合雷达高度计天线相位中心的位置得到地球表面的高度(文献[1]：Y.Zhang,J.Jiang,H.Zhang,and D.Zhang,“Spaceborne imaging altimeter for topographic mapping,”2000IEEEInternational Geoscience and Remote Sensing Symposium,2000,pp.2349-2351；文献[2]：Y.Zhang,X.Zhang,X.Meng,W.Luo,Z.Zhou,and J.Jiang,“An interferometricimaging altimeter applied for both ocean and land observation,”2007IEEEInternational Geoscience and Remote Sensing Symposium,2007,pp.3821-3824)。
同传统雷达高度计相比，星载干涉成像高度计采用偏离天顶点干涉合成孔径技术得到目标高度，因此具备宽刈幅、高分辨率的优点；同传统合成孔径雷达相比，星载干涉成像高度计入射角更小，水体表面在小入射角下呈准镜面散射，散射系数和信噪比均高于大角度入射情况，因此更适合观测海洋、湖泊、河流等水体表面。2016年发射的天宫二号空间实验室上装载了国际上第一个星载干涉成像高度计，该高度计工作于Ku波段，具备海洋和陆地的三维测绘能力(文献[3]：Y.Zhang,J.Jiang,X.Zhang,K.Xu,J.Yan,C.Jiang,et al.,“Design and preliminary experiment of china imaging altimeter,”MicrowaveRemote Sensing of the Atmosphere and Environment III,2003,pp.190-199；文献[4]：X.Dong,Y.Zhang,and W.Zhai,“Design and algorithms of the Tiangong-2interferometric imaging radar altimeter processor,”2017Progress InElectromagnetics Research Symposium-Spring(PIERS 2017),2017,pp.3802-3803)，对海洋观测和研究具有重要的意义。
干涉成像高度计的高度测量误差包括系统性误差和随机性误差，系统性误差包括仪器测量漂移、基线测量误差等因素，可通过定标进行校正；随机性误差包是指干涉相位噪声，是由系统热噪声、几何去相关和处理误差导致的。干涉相位噪声无法通过定标方法消除，只能通过干涉滤波方法降低相位噪声均方差，是影响干涉雷达高度计测量精度的主要因素。
传统的雷达干涉相位滤波方法包括多视处理、中值滤波、高斯滤波、Lee滤波和Goldstein滤波等(文献[5]：A.L.B.Candeias,J.C.Mura,L.V.Dutra,and J.R.Moreira,“Interferogram phase noise reduction using morphological and modified medianfilters,”1995IEEE Geoscience and Remote Sensing Symposium,1995,pp.166-168；Jong-Sen,L.,文献[6]：Papathanassiou,K.P.,Ainsworth,T.L.,Grunes,M.R.,andReigber,A.:“A new technique for noise filtering of SAR interferometric phaseimages”,IEEE Transactions on Geoscience and Remote Sensing,1998,36,(5),pp.1456-1465；文献[7]：Goldstein,R.M.,and Werner,C.L.:“Radar interferogramfiltering for geophysical applications”,Geophysical Research Letters,1998,25,(21),pp.4035-4038)，其中空间多视滤波法通过相邻/相关像素的平均降低噪声方差，是雷达干涉相位滤波中最常用的方法。成像雷达高度计采用小入射角干涉技术，通常在整个观测刈幅范围内入射角变化大，最大入射角可达最小入射角的3～5倍，导致交轨向的空间分辨率和高度误差变化剧烈。传统雷达干涉相位滤波方法是针对大角度干涉合成孔径雷达设计的，在小角度干涉情况下，会出现高度误差沿交轨向不均匀的问题，降低高度测量精度。
发明内容
本发明的目的在于克服小角度干涉下高度误差非均匀的问题，从而提供一种星载干涉成像高度计的相位滤波方法及高度测量方法，在小角度观测的情况下，该方法能够降低干涉相位噪声，实现星载干涉成像高度计的高精度的高度测量。
为了实现上述目的，本发明提出了一种星载干涉成像高度计的相位滤波方法，所述方法包括：
对星载干涉成像高度计两个接收通道的数据进行合成孔径成像处理，得到两幅复图像，以其中一幅图像作为基准对两幅复图像进行配准，再计算配准后两幅图像的相干系数；
依次计算图像中每个像元的几何位置、观测角和多视视数；由此确定每个像元的方位向滤波窗口大小和距离向滤波窗口大小，并确定每个像元的方位向滤波区间和距离向滤波区间；
将两个复图像进行共轭相乘，同时去除地平相位，然后在方位向滤波区间和距离向滤波区间对每个像元分别进行滤波，最终得到滤波后的干涉相位。
作为上述方法的一种改进，所述方法具体包括：
步骤1)对星载干涉成像高度计两个接收通道的数据进行合成孔径成像处理，得到两幅复数图像s1(na,nr)和s2(na,nr)，其中nr和na两个图像像素横纵坐标的序号，以s1(na,nr)作为基准对两幅复图像进行配准，配准后再计算相干系数γ(na,nr)：
其中，Πa(na)为方位向慢时间的有效区间，Πr(nr)为距离向快时间的有效区间，
Πa(na)＝[na-Ma/2,na+Ma/2]，Πr(nr)＝[nr-Mr/2,nr+Mr/2]
其中，Ma和Mr分别是方位向和距离向的窗口；
步骤2)利用距离-多普勒方程分别得到s1(na,nr)和s2(na,nr)中每个像元的几何位置Pt(na,nr)；计算s1(na,nr)和s2(na,nr)中每个像元的观测角θ(na,nr)；
以s1(na,nr)为基准将两幅图像进行配准后，s1(na,nr)和s2(na,nr)图像每个像元的几何位置和观测角相同；
步骤3)对于给定的高度误差门限值σh，计算s1(na,nr)和s2(na,nr)中每个像元的多视视数L(na,nr)为：
其中，λ为波长，r1(na,nr)为第一天线相位中心到目标像元(na,nr)的距离，α是基线倾角，B为两个天线的基线长度；
步骤4)确定s1(na,nr)和s2(na,nr)中每个像元的方位向滤波窗口大小Wa(na,nr)和距离向滤波窗口大小Wr(na,nr)，由此确定每个像元的方位向滤波区间Γa(na,nr)和距离向滤波区间Γr(na,nr)；
步骤5)对两个通道复图像s1(na,nr)和s2(na,nr)进行共轭相乘，同时去除地平相位，然后对每个像元分别进行滤波，最终得到滤波后的干涉相位为：
其中，为(i,j)处的地平相位。
作为上述方法的一种改进，所述步骤2)具体为：
步骤2-1)建立以下方程组计算s1(na,nr)和s2(na,nr)图像像元的几何位置Pt(na,nr)，
其中，Pt(na,nr)为(na,nr)处目标像元的坐标，三维坐标为xt(na,nr)、yt(na,nr)和zt(na,nr)；S1(na)为第一天线在方位序号na处的坐标；V1(na)为第一天线在方位序号na处的速度；Vt(na,nr)为(na,nr)处目标的速度；fdc为多普勒中心；Re为地球赤道半径；h(na,nr)为(na,nr)处目标相对于地球模型的高度，Rp为地球极地半径，Rp＝(1-f)(Re+h)f为地球扁率因子；
采用牛顿法迭代求解上述方程组得到s1(na,nr)图像像元的几何位置Pt(na,nr)；
步骤2-2)分别计算s1(na,nr)和s2(na,nr)中每个像元的观测角θ(na,nr)：
作为上述方法的一种改进，所述步骤3)中的高度误差门限值σh为：
作为上述方法的一种改进，所述步骤4)具体包括：
步骤4-1)确定距离向和方位向滤波窗口的过程为：
其中，ρa(na,nr)为每个像元的原始方位分辨率，ρr(na,nr)为每个像元的原始地距分辨率：
其中，D为天线方位向有效口径，Bw为发射信号带宽，θinc(na,nr)为像元的本地入射角：
步骤4-2)每个像元的方位向滤波区间Γa(na,nr)和距离向滤波区间Γr(na,nr)为：
作为上述方法的一种改进，所述的步骤5)的地平相位的计算公式为
其中，S2(na)为第二天线在方位序号na处的坐标。
基于上述相位滤波方法，本发明还提出了一种星载干涉成像高度计的高度测量方法，所述方法包括：
步骤S1)对滤波后的干涉相位进行解缠绕，对于图像中的每个像元，利用该像元对应的距离测量值和解缠后相位φ(na,nr)，根据第一天线和第二天线的位置和速度，求出该像元的绝对位置Pt′(na,nr)；
步骤S2)将像元的绝对位置Pt′(na,nr)转换成经度、纬度和高度坐标，得到每个像元的高度h(na,nr)；
步骤S3)建立经纬度均匀分布的网格，将图像高度坐标二维插值到均匀经纬度坐标下，针对目标经纬度所在区域和分辨率要求，选用滤波窗口对高度进行滤波。
作为上述方法的一种改进，所述步骤S1)中通过求解下面的方程组计算该像元的绝对位置Pt′(na,nr)：
其中，λ为波长，S1(na)为第一天线在方位序号na处的坐标，S2(na)为第二天线在方位序号na处的坐标；Vt(na,nr)为(na,nr)处目标的速度；r1(na,nr)为(na,nr)处目标的距离测量值；fdc为多普勒中心。
作为上述方法的一种改进，所述的步骤S3)分成海洋区域和陆地区域分别进行高度滤波：对于海洋区域，所述滤波窗口选用公里级别高斯窗口进行多视平均；对于陆地区域，所述滤波窗口选用百米级别高斯窗口进行平均。
本发明的优点在于：
1、本发明为实现星载干涉成像高度计的高精度高度测量提供了新的技术途径：它根据不同目标的相干系数与高度误差的关系，利用真实测量的相干系数计算所需的滤波窗口，每个像元的滤波窗口是独立计算得到，一方面可以使所有像元的高度误差保持一致，另一方面可以通过并行计算快速实现；
2、本发明方法针对星载干涉成像高度计小角度测量的条件下，实现了刈幅内高度误差的均匀分布，计算效率高，易于实现，对于实际系统有重要的现实意义。
附图说明
图1为本发明提出的星载干涉成像高度计的高度测量方法流程图；
图2为星载干涉成像高度计的观测几何关系示意图；
图3为天宫二号成像高度计未滤波的干涉相位条纹图；
图4为天宫二号成像高度计相干系数图；
图5为天宫二号成像高度计去地平并解缠后相位图；
图6为采用传统多视滤波方法得到的海面高度误差分布图；
图7为采用本发明提出的一种星载干涉成像高度计的高度测量方法得到的海面高度误差分布图；
图8为本发明提出的星载干涉成像高度计的高度测量方法得到的天宫二号成像高度计三维海面高度图。
具体实施方式
为了更好的理解本发明的技术方案，以下结合附图对本发明的实施方式作进一步的描述。
如图1所示，本发明提出的一种干涉雷达高度计的高度测量方法，该方法具体步骤包括：
步骤1)对星载干涉成像高度计两个接收通道的数据进行合成孔径成像处理，得到两幅复数图像s1(na,nr)和s2(na,nr)，其中nr和na两个图像像素横纵坐标的序号，以s1(na,nr)作为基准对两幅复图像进行配准，配准后再计算相干系数γ(na,nr)：
其中，Πa(na)为方位向慢时间的有效区间，Πr(nr)为距离向快时间的有效区间，
Πa(na)＝[na-Ma/2,na+Ma/2]，Πr(nr)＝[nr-Mr/2,nr+Mr/2] (2)
其中，Ma和Mr分别是方位向和距离向的窗口；
步骤2)利用距离-多普勒方程分别得到s1(na,nr)和s2(na,nr)中每个像元的几何位置Pt(na,nr)；计算s1(na,nr)和s2(na,nr)中每个像元的观测角θ(na,nr)；
由于步骤1)以第一幅图像为基准完成两幅图像的配准，故两幅图像的像元位置和观测角度相同。
其中各变量的定义如下：
Pt(na,nr)为(na,nr)处目标像元的坐标，三维坐标为xt(na,nr)、yt(na,nr)和zt(na,nr)；
S1(na)为第一天线在方位序号na处的坐标；
r1(na,nr)为(na,nr)处目标的距离测量值；
Vt(na,nr)为(na,nr)处目标的速度；
V1(na)为第一天线在方位序号na处的速度；
fdc为多普勒中心；
Re为地球赤道半径；
h(na,nr)为(na,nr)处目标相对于地球模型的高度；
Rp为地球极地半径，Rp＝(1-f)(Re+h)，f为地球扁率因子；
对每个像元均按(3)中方程组进行求解得到三维位置Pt(na,nr)，采用牛顿法迭代求解，每个像元的求解相互独立，利用并行化处理提高计算效率；
根据星载干涉成像高度计的观测几何，如图2所示，计算每个像元的观测角θ(na,nr)，计算公式为，
步骤3)对于给定的高度误差门限值σh，计算s1(na,nr)和s2(na,nr)中每个像元的多视视数L(na,nr)为：
其中，λ为波长，r1(na,nr)为第一天线相位中心到(na,nr)处目标像元的距离，α是基线倾角，B为两个天线的基线长度；
根据相位噪声方差σφ与相干系数的关系：
相位噪声φ与高度误差δh的关系是：
其中，B⊥为有效基线长度，如图2所示，B⊥＝Bcos(θ(na,nr)-α)，则高度误差门限值σh与相干系数的关系是：
步骤4)确定s1(na,nr)和s2(na,nr)中每个像元的方位向滤波窗口大小Wa(na,nr)和距离向滤波窗口大小Wr(na,nr)，由此确定每个像元的方位向滤波区间Γa(na,nr)和距离向滤波区间Γr(na,nr)；
首先，二维滤波窗口的大小Wa(na,nr)×Wr(na,nr)为：
Wa(na,nr)×Wr(na,nr)＝L(na,nr)
其次，根据每个像元的原始地距分辨率ρr(na,nr)和原始方位分辨率ρa(na,nr)，使得滤波后的分辨率在两个方向上相同
Wa(na,nr)ρa(na,nr)＝Wr(na,nr)ρr(na,nr)
其中ρa(na,nr)为天线有效孔径D的一半，ρr(na,nr)取决于发射信号带宽Bw与像元的本地入射角θinc(na,nr)：
得到
则每个像元的方位向滤波区间Γa(na,nr)和距离向滤波区间Γr(na,nr)为：
步骤5)对步骤2)配准后的两个通道复图像s1(na,nr)和s2(na,nr)进行共轭相乘，同时去除地平相位，地平相位的计算公式为
然后对每个像元分别进行滤波，每个像元的滤波窗口由步骤4)计算得到，最终得到干涉相位
步骤6)对滤波后的干涉相位进行解缠绕，对于图像中的每个像元，利用该像元对应的距离测量值和解缠后相位φ(na,nr)，根据第一天线和第二天线的位置和速度，求出该像元的绝对位置Pt′(na,nr)；
方程组的求解可以有两种处理方式：1数值解法，利用牛顿法迭代求解；2.解析法，直接根据解析解计算得到结果；为了获得高精度的三维位置信息，需要对每个像元都进行求解上面的方程组，在数据量大的情况下计算复杂度高；提高计算速度的手段是并行计算，由于每个像素的计算方法相同，只是数据不同，采用单指令多数据的并行计算方法可以快速的求解出每个像元的三维位置。
步骤7)将像元的绝对位置Pt′(na,nr)转换成经度、纬度和高度坐标，得到每个像元的高度h(na,nr)，由于入射角在整个观测刈幅内变化大，远地点入射角可达近地点入射角的3～5倍，此时高度h(na,nr)的地理分布非均匀；
步骤8)建立经纬度均匀分布的网格，将图像高度坐标二维插值到均匀经纬度坐标下，针对目标经纬度所在区域，分成海洋区域和陆地区域分别进行高度滤波：对于海洋区域，由于目标高度变化缓慢且均匀，滤波窗口选用公里级别高斯窗口进行多视平均；对于陆地区域，由于目标高度起伏剧烈，滤波窗口选用百米级别高斯窗口进行平均；经过高度滤波后，得到观测区域的三维高度图像。
如图3、图4、图5、图6、图7和图8所示，图3为天宫二号成像高度计未滤波的干涉相位条纹图；图4为天宫二号成像高度计相干系数图；图5为天宫二号成像高度计去地平并解缠后相位图；图6为采用传统多视滤波方法得到的海面高度误差分布图，该图未做高度滤波，因此波动较大；图7为本发明提出的一种星载干涉成像高度计的相位滤波方法得到的海面高度误差分布图，该图未做高度滤波，分辨率高，但波动较大，后续的高度滤波会降低分辨率减小波动；图8为本发明提出的一种星载干涉成像高度计的相位滤波方法得到的天宫二号成像高度计三维海面高度图，该图经过高度滤波，得到平滑的三维海面高度。所用的数据为天宫二号三维成像微波高度计对太平洋海域的观测结果。首先利用合成孔径方法进行成像，经双通道干涉得到图3的干涉相位图和图4的相干系数图，去除地平相位后经过解缠得到图5；然后利用本发明的星载干涉成像高度计的相位滤波方法进行相位滤波消除相位噪声，高度误差沿距离向分布如图7所示；最后利用本发明的星载干涉成像高度计的相位滤波方法进行二次高程滤波得到图8中的三维高度图，分辨率为5公里。
从图3、图4和图5中可以看出，干涉相位在近地点变化速率快，远地点变换速率慢，因为近地点入射角约为1～3度，而远地点入射角约为6～8度，二者差距大，这导致相位噪声在刈幅内分布不均匀；相干系数受到入射角、天线方向图、双天线几何关系等因素的影响，在整个刈幅范围内非均匀变化。从图6、图7和图8可以看出，采用传统的固定窗口多视滤波方法对干涉成像高度计进行处理，得到的高度误差变化在刈幅内不均匀，这是因为近地点和远地点的相位噪声分布不同；采用本发明的相位滤波方法，可以得到均匀的高度误差分布，经过二次高度滤波后可以得到平滑的三维海面拓扑结构。
本发明针对星载干涉成像高度计，在小入射角干涉反演高度的情况下，保证了三维高度的均匀误差分布，实现了高精度的三维高度处理。理论上，本发明所提出的方法可以实现任意高度误差需求下的滤波处理，同时像元间计算可以通过并行加速，对于星载干涉成像高度计准确测量海面高度有重要的现实意义。
Claims (9)
1.一种星载干涉成像高度计的相位滤波方法，所述方法包括：
对星载干涉成像高度计两个接收通道的数据进行合成孔径成像处理，得到两幅复图像，以其中一幅图像作为基准对两幅复图像进行配准，再计算配准后两幅图像的相干系数；
依次计算图像中每个像元的几何位置、观测角和多视视数；由此确定每个像元的方位向滤波窗口大小和距离向滤波窗口大小，并确定每个像元的方位向滤波区间和距离向滤波区间；
将两个复图像进行共轭相乘，同时去除地平相位，然后在方位向滤波区间和距离向滤波区间对每个像元分别进行滤波，最终得到滤波后的干涉相位。
2.根据权利要求1所述的星载干涉成像高度计的相位滤波方法，其特征在于，所述方法具体包括：
步骤1)对星载干涉成像高度计两个接收通道的数据进行合成孔径成像处理，得到两幅复数图像s1(na,nr)和s2(na,nr)，其中nr和na两个图像像素横纵坐标的序号，以s1(na,nr)作为基准对两幅复图像进行配准，配准后再计算相干系数γ(na,nr)：
其中，Πa(na)为方位向慢时间的有效区间，Πr(nr)为距离向快时间的有效区间，
Πa(na)＝[na-Ma/2,na+Ma/2]，Πr(nr)＝[nr-Mr/2,nr+Mr/2]
其中，Ma和Mr分别是方位向和距离向的窗口；
步骤2)利用距离-多普勒方程分别得到s1(na,nr)和s2(na,nr)中每个像元的几何位置Pt(na,nr)；计算s1(na,nr)和s2(na,nr)中每个像元的观测角θ(na,nr)；
以s1(na,nr)为基准将两幅图像进行配准后，s1(na,nr)和s2(na,nr)图像每个像元的几何位置和观测角相同；
步骤3)对于给定的高度误差门限值σh，计算s1(na,nr)和s2(na,nr)中每个像元的多视视数L(na,nr)为：
其中，λ为波长，r1(na,nr)为第一天线相位中心到目标像元(na,nr)的距离，α是基线倾角，B为两个天线的基线长度；
步骤4)确定s1(na,nr)和s2(na,nr)中每个像元的方位向滤波窗口大小Wa(na,nr)和距离向滤波窗口大小Wr(na,nr)，由此确定每个像元的方位向滤波区间Γa(na,nr)和距离向滤波区间Γr(na,nr)；
步骤5)对两个通道复图像s1(na,nr)和s2(na,nr)进行共轭相乘，同时去除地平相位，然后对每个像元分别进行滤波，最终得到滤波后的干涉相位为：
其中，为(i,j)处的地平相位。
3.根据权利要求2所述的星载干涉成像高度计的相位滤波方法，其特征在于，所述步骤2)具体为：
步骤2-1)建立以下方程组计算s1(na,nr)和s2(na,nr)图像像元的几何位置Pt(na,nr)，
其中，Pt(na,nr)为(na,nr)处目标像元的坐标，三维坐标为xt(na,nr)、yt(na,nr)和zt(na,nr)；S1(na)为第一天线在方位序号na处的坐标；V1(na)为第一天线在方位序号na处的速度；Vt(na,nr)为(na,nr)处目标的速度；fdc为多普勒中心；Re为地球赤道半径；h(na,nr)为(na,nr)处目标相对于地球模型的高度，Rp为地球极地半径，Rp＝(1-f)(Re+h)f为地球扁率因子；
采用牛顿法迭代求解上述方程组得到s1(na,nr)图像像元的几何位置Pt(na,nr)；
步骤2-2)分别计算s1(na,nr)和s2(na,nr)中每个像元的观测角θ(na,nr)：
4.根据权利要求3所述的星载干涉成像高度计的相位滤波方法，其特征在于，所述步骤3)中的高度误差门限值σh为：
5.根据权利要求2所述的星载干涉成像高度计的相位滤波方法，其特征在于，所述步骤4)具体包括：
步骤4-1)确定距离向和方位向滤波窗口的过程为：
其中，ρa(na,nr)为每个像元的原始方位分辨率，ρr(na,nr)为每个像元的原始地距分辨率：
其中，D为天线方位向有效口径，Bw为发射信号带宽，θinc(na,nr)为像元的本地入射角：
步骤4-2)每个像元的方位向滤波区间Γa(na,nr)和距离向滤波区间Γr(na,nr)为：
6.根据权利要求3所述的星载干涉成像高度计的相位滤波方法，其特征在于，所述的步骤5)的地平相位的计算公式为
其中，S2(na)为第二天线在方位序号na处的坐标。
7.一种基于权利要求1-6之一所述的相位滤波方法实现的星载干涉成像高度计的高度测量方法，所述方法包括：
步骤S1)对滤波后的干涉相位进行解缠绕，对于图像中的每个像元，利用该像元对应的距离测量值和解缠后相位φ(na,nr)，根据第一天线和第二天线的位置和速度，求出该像元的绝对位置Pt′(na,nr)；
步骤S2)将像元的绝对位置Pt′(na,nr)转换成经度、纬度和高度坐标，得到每个像元的高度h(na,nr)；
步骤S3)建立经纬度均匀分布的网格，将图像高度坐标二维插值到均匀经纬度坐标下，针对目标经纬度所在区域和分辨率要求，选用滤波窗口对高度进行滤波。
8.根据权利要求7所述的星载干涉成像高度计的高度测量方法，其特征在于，所述步骤S1)中通过求解下面的方程组计算该像元的绝对位置Pt′(na,nr)：
其中，λ为波长，S1(na)为第一天线在方位序号na处的坐标，S2(na)为第二天线在方位序号na处的坐标；Vt(na,nr)为(na,nr)处目标的速度；r1(na,nr)为(na,nr)处目标的距离测量值；fdc为多普勒中心。
9.根据权利要求7或8所述的星载干涉成像高度计的高度测量方法，其特征在于，所述的步骤S3)分成海洋区域和陆地区域分别进行高度滤波：对于海洋区域，所述滤波窗口选用公里级别高斯窗口进行多视平均；对于陆地区域，所述滤波窗口选用百米级别高斯窗口进行平均。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810777772.4A CN109116321B (zh) | 2018-07-16 | 2018-07-16 | 一种星载干涉成像高度计的相位滤波方法及高度测量方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810777772.4A CN109116321B (zh) | 2018-07-16 | 2018-07-16 | 一种星载干涉成像高度计的相位滤波方法及高度测量方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109116321A true CN109116321A (zh) | 2019-01-01 |
| CN109116321B CN109116321B (zh) | 2019-09-24 |
Family
ID=64862164
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810777772.4A Active CN109116321B (zh) | 2018-07-16 | 2018-07-16 | 一种星载干涉成像高度计的相位滤波方法及高度测量方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109116321B (zh) |
Cited By (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110764087A (zh) * | 2019-10-15 | 2020-02-07 | 中国科学院国家空间科学中心 | 一种基于干涉成像高度计的海面风向反加权反演方法 |
| CN110907132A (zh) * | 2019-12-13 | 2020-03-24 | 中国人民解放军军事科学院国防科技创新研究院 | 一种波向探测方法、系统、设备、介质 |
| CN111505592A (zh) * | 2020-04-26 | 2020-08-07 | 中国科学院国家空间科学中心 | 一种利用底视差分干涉实现星载雷达横滚角测定的方法 |
| CN112308889A (zh) * | 2020-10-23 | 2021-02-02 | 香港理工大学深圳研究院 | 一种利用矩形及扁率信息的点云配准方法及存储介质 |
| CN113406628A (zh) * | 2021-05-12 | 2021-09-17 | 中国科学院国家空间科学中心 | 一种用于干涉成像高度计的数据重采样方法 |
| CN113589282A (zh) * | 2021-07-12 | 2021-11-02 | 中国科学院国家空间科学中心 | 一种基于图像域变换的星载干涉成像高度计去平地效应方法 |
| CN114359876A (zh) * | 2022-03-21 | 2022-04-15 | 成都奥伦达科技有限公司 | 一种车辆目标识别方法及存储介质 |
| CN118191741A (zh) * | 2024-03-25 | 2024-06-14 | 中国科学院国家空间科学中心 | 一种开环模式下合成孔径雷达高度计信号处理方法及系统 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105487065A (zh) * | 2016-01-08 | 2016-04-13 | 香港理工大学深圳研究院 | 一种时序星载雷达数据处理方法和装置 |
| CN106249236A (zh) * | 2016-07-12 | 2016-12-21 | 北京航空航天大学 | 一种星载InSAR长短基线图像联合配准方法 |
-
2018
- 2018-07-16 CN CN201810777772.4A patent/CN109116321B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105487065A (zh) * | 2016-01-08 | 2016-04-13 | 香港理工大学深圳研究院 | 一种时序星载雷达数据处理方法和装置 |
| CN106249236A (zh) * | 2016-07-12 | 2016-12-21 | 北京航空航天大学 | 一种星载InSAR长短基线图像联合配准方法 |
Non-Patent Citations (2)
| Title |
|---|
| XIAO DONG ET AL.: "Design and Algorithms of the Tiangong-2 Interferometric Imaging Radar Altimeter Processor", 《2017 PROGRESS IN ELECTROMAGNETICS RESEARCH SYMPOSIUM—SPRING（PIERS）》 * |
| 全刚: "干涉SAR成象算法研究", 《中国优秀硕士学位论文全文数据库 信息科技辑》 * |
Cited By (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110764087A (zh) * | 2019-10-15 | 2020-02-07 | 中国科学院国家空间科学中心 | 一种基于干涉成像高度计的海面风向反加权反演方法 |
| CN110764087B (zh) * | 2019-10-15 | 2021-08-31 | 中国科学院国家空间科学中心 | 一种基于干涉成像高度计的海面风向反加权反演方法 |
| CN110907132A (zh) * | 2019-12-13 | 2020-03-24 | 中国人民解放军军事科学院国防科技创新研究院 | 一种波向探测方法、系统、设备、介质 |
| CN111505592A (zh) * | 2020-04-26 | 2020-08-07 | 中国科学院国家空间科学中心 | 一种利用底视差分干涉实现星载雷达横滚角测定的方法 |
| CN112308889A (zh) * | 2020-10-23 | 2021-02-02 | 香港理工大学深圳研究院 | 一种利用矩形及扁率信息的点云配准方法及存储介质 |
| CN113406628A (zh) * | 2021-05-12 | 2021-09-17 | 中国科学院国家空间科学中心 | 一种用于干涉成像高度计的数据重采样方法 |
| CN113589282A (zh) * | 2021-07-12 | 2021-11-02 | 中国科学院国家空间科学中心 | 一种基于图像域变换的星载干涉成像高度计去平地效应方法 |
| CN113589282B (zh) * | 2021-07-12 | 2023-08-08 | 中国科学院国家空间科学中心 | 一种基于图像域变换的星载干涉成像高度计去平地效应方法 |
| CN114359876A (zh) * | 2022-03-21 | 2022-04-15 | 成都奥伦达科技有限公司 | 一种车辆目标识别方法及存储介质 |
| CN118191741A (zh) * | 2024-03-25 | 2024-06-14 | 中国科学院国家空间科学中心 | 一种开环模式下合成孔径雷达高度计信号处理方法及系统 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109116321B (zh) | 2019-09-24 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109116321B (zh) | 一种星载干涉成像高度计的相位滤波方法及高度测量方法 | |
| CN107544095B (zh) | 一种地面三维激光点云与探地雷达图像融合的方法 | |
| US9417323B2 (en) | SAR point cloud generation system | |
| CN101738614B (zh) | 基于时空图像序列的isar目标转动估计方法 | |
| CN101078769B (zh) | 单次全极化合成孔径雷达图像反演数字高程模型的方法 | |
| US7796829B2 (en) | Method and system for forming an image with enhanced contrast and/or reduced noise | |
| CN112305510B (zh) | 一种基于dem匹配的合成孔径雷达影像几何定标方法 | |
| US8159384B2 (en) | Method for examining an ice region or dry region using radar echo sounding | |
| JP2003500658A (ja) | 市街地領域及び地滑り地帯の運動に関するレーダー測定のための手順 | |
| Morrison et al. | Laboratory measurement of the DInSAR response to spatiotemporal variations in soil moisture | |
| CN102288964A (zh) | 一种星载高分辨率合成孔径雷达的成像处理方法 | |
| CN103869316A (zh) | 基于稀疏表征的前视阵列sar超分辨成像方法 | |
| CN114545411B (zh) | 一种基于工程实现的极坐标格式多模高分辨sar成像方法 | |
| US6677885B1 (en) | Method for mitigating atmospheric propagation error in multiple pass interferometric synthetic aperture radar | |
| CN106526593A (zh) | 基于sar严密成像模型的子像素级角反射器自动定位方法 | |
| CN104062657A (zh) | 一种合成孔径雷达sar的广义极坐标成像方法 | |
| CN112986949B (zh) | 针对角反射器的sar高精度时序形变监测方法和装置 | |
| Martorella | Introduction to inverse synthetic aperture radar | |
| CN117148352B (zh) | 一种角度唯一性约束的阵列干涉sar三维成像方法 | |
| KR20120009186A (ko) | Ｓａｒ데이터를 이용하여 수치고도모형을 제작하기 위한 방법 | |
| CN110554377B (zh) | 基于多普勒中心偏移的单通道sar二维流场反演方法及系统 | |
| CN112346027B (zh) | 用于确定合成孔径雷达图像的散射特性的方法和系统 | |
| CN112179314B (zh) | 一种基于三维网格投影的多角度sar高程测量方法及系统 | |
| Parker et al. | Preliminary experimental results for RF tomography using distributed sensing | |
| Marghany | DEM reconstruction of coastal geomorphology from DINSAR |
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 |