发明名称 一种基于曲面投影的InSAR成像方法
摘要 本发明公开了一种基于曲面投影的InSAR成像方法,该方法通过在低精度的高程曲面上进行反向投影成像,将不同雷达天线的测量数据投影到以低精度高程曲面建立的InSAR成像空间中,提高了InSAR成像精度和干涉相位保持精度,同时提取干涉相位以后,陡变地形区域的干涉条纹稀疏,后续相位解缠的难度降低,相位解缠的精度提高,使得高程反演的精度提高,为InSAR高质量干涉相位生成提供了一种新方法。
申请公布号 CN103018740B 申请公布日期 2014.07.16
申请号 CN201210475166.X 申请日期 2012.11.21
申请人 电子科技大学 发明人 张晓玲;付涛
分类号 G01S13/90(2006.01)I 主分类号 G01S13/90(2006.01)I
代理机构 电子科技大学专利中心 51203 代理人 曾磊
主权项 一种基于曲面投影的InSAR成像方法,其特征是它包括以下几个步骤:步骤1、初始化InSAR成像系统参数InSAR成像空间由InSAR成像空间中的三个相互正交的坐标基确定,定义与雷达平台速度方向平行并在地平面内的单位向量作为InSAR成像空间的第一个坐标基,记做<img file="FDA0000485737230000011.GIF" wi="77" he="69" />该坐标基方向为方位向;定义在地平面内,并与InSAR成像空间的第一个坐标基<img file="FDA0000485737230000012.GIF" wi="47" he="69" />垂直的单位向量作为InSAR成像空间的第二个坐标基,记做<img file="FDA0000485737230000013.GIF" wi="78" he="69" />该坐标基方向为距离向;定义垂直于地平面向上的单位向量作为InSAR成像空间的第三个坐标基,记做<img file="FDA0000485737230000014.GIF" wi="79" he="69" />该坐标基方位为高度向;InSAR雷达平台包含两组天线,即主天线和副天线,两组天线之间的距离为基线长,记做B<sub>l</sub>,主天线发射脉冲信号,经过T<sub>d</sub>时间的延迟,主天线和副天线同时接收回波延迟信号;雷达平台主天线接收的回波数据,记做<img file="FDA0000485737230000015.GIF" wi="86" he="70" />雷达平台副天线接收的回波数据,记做<img file="FDA0000485737230000016.GIF" wi="75" he="70" />其中<img file="FDA0000485737230000017.GIF" wi="52" he="70" />和<img file="FDA0000485737230000018.GIF" wi="56" he="70" />均为二维矩阵,第一维均对应方位向,第二维均对应距离向,即二维矩阵<img file="FDA0000485737230000019.GIF" wi="48" he="70" />和<img file="FDA00004857372300000110.GIF" wi="52" he="70" />的行存储的是方位向数据,二维矩阵<img file="FDA00004857372300000111.GIF" wi="54" he="70" />和<img file="FDA00004857372300000112.GIF" wi="50" he="70" />的列存储的是距离向数据;初始化InSAR成像系统参数包括:雷达系统工作的信号波长,记做λ,雷达平台主天线发射信号带宽,记做B,雷达平台主天线发射脉冲时宽,记做T<sub>r</sub>,雷达平台接收系统采样频率,记做F<sub>s</sub>,雷达系统脉冲重复频率,记做PRF,雷达平台一个合成孔径长度内的慢时刻个数,记做N<sub>l</sub>,雷达平台速度矢量,记做<img file="FDA00004857372300000113.GIF" wi="73" he="64" />雷达平台主天线初始位置矢量,记做<img file="FDA00004857372300000114.GIF" wi="153" he="81" />雷达平台副天线初始位置矢量,记做<img file="FDA00004857372300000115.GIF" wi="149" he="80" />场景参考点位置矢量,记做<img file="FDA00004857372300000116.GIF" wi="82" he="73" />雷达系统距离向采样点数,记做Nr,雷达系统方位向采样点数,记做N<sub>a</sub>,场景距离向散射点间隔,记做d<sub>r</sub>,场景方位向散射点间隔,记做d<sub>a</sub>,场景参考点到雷达平台主天线各慢时刻天线相位中心的最短距离,记做R<sub>mc</sub>,场景参考点到雷达平台副天线各慢时刻天线相位中心的最短距离,记做R<sub>sc</sub>,场景参考点在雷达平台主天线回波数据和雷达平台副天线回波数据中的距离门相同,距离门位置记做I<sub>c</sub>;对于场景散射点P(a,r),a表示散射点位于场景方位向的第a个位置,a=1,...,s<sub>a</sub>,s<sub>a</sub>为场景方位向总的散射点数,r表示散射点位于场景距离向的第r个位置,r=1,...,s<sub>r</sub>,s<sub>r</sub>为场景距离向总的散射点数,场景散射点P(a,r)的低精度高程,记做H(a,r);上述参数中,雷达系统工作的信号波长λ,雷达平台主天线发射的信号带宽B,雷达平台主天线发射的脉冲时宽T<sub>r</sub>,雷达平台接收系统的采样频率F<sub>s</sub>,雷达系统的脉冲重复频率PRF,两组天线之间的基线长B<sub>l</sub>以及接收系统接收波门相对于发射信号发射波门的延迟T<sub>d</sub>在InSAR雷达系统设计过程中已经确定;雷达平台一个合成孔径长度内的慢时刻个数N<sub>l</sub>,雷达平台速度矢量<img file="FDA0000485737230000021.GIF" wi="73" he="64" />雷达平台主天线初始位置矢量<img file="FDA0000485737230000022.GIF" wi="168" he="80" />雷达平台副天线初始位置矢量<img file="FDA0000485737230000023.GIF" wi="142" he="80" />场景参考点位置矢量<img file="FDA0000485737230000024.GIF" wi="76" he="73" />雷达系统距离向采样点数N<sub>r</sub>,雷达系统方位向采样点数N<sub>a</sub>,场景距离向散射点间隔d<sub>r</sub>,场景方位向散射点间隔d<sub>a</sub>,场景参考点到雷达平台主天线各慢时刻天线相位中心的最短距离R<sub>mc</sub>,场景参考点到雷达平台副天线各慢时刻天线相位中心的最短距离R<sub>sc</sub>,场景参考点在雷达平台主天线回波数据和雷达平台副天线回波数据中的距离门位置I<sub>c</sub>以及场景散射点P(a,r)的低精度高程H(a,r),a表示散射点位于场景方位向的第a个位置,a=1,...,s<sub>a</sub>,s<sub>a</sub>为场景方位向总的散射点数,r表示散射点位于场景距离向的第r个位置,r=1,...,s<sub>r</sub>,s<sub>r</sub>为场景距离向总的散射点数,在InSAR雷达成像观测方案设计中已经确定;根据InSAR雷达系统方案和InSAR雷达成像观测方案,以上基于曲面投影的InSAR成像方法需要的初始化成像系统参数均为已知;步骤2:InSAR原始回波数据进行距离压缩采用传统的合成孔径雷达标准距离压缩方法对雷达平台主天线距离向回波数据<img file="FDA0000485737230000025.GIF" wi="62" he="72" />进行压缩,得到平台主天线距离压缩后数据,记做<img file="FDA0000485737230000026.GIF" wi="92" he="79" />采用传统的合成孔径雷达标准距离压缩方法对雷达平台副天线距离向回波数据<img file="FDA0000485737230000031.GIF" wi="58" he="70" />进行压缩,得到平台副天线距离压缩后数据,记做<img file="FDA0000485737230000032.GIF" wi="78" he="79" />步骤3、计算InSAR成像空间中散射点的距离史采用公式<img file="FDA0000485737230000033.GIF" wi="482" he="72" />计算得到雷达平台主天线第n个慢时刻的天线相位中心矢量<img file="FDA0000485737230000034.GIF" wi="140" he="72" />采用公式<img file="FDA0000485737230000035.GIF" wi="464" he="72" />计算得到雷达平台副天线第n个慢时刻的天线相位中心矢量<img file="FDA0000485737230000036.GIF" wi="138" he="72" />其中n表示第n个慢时刻,n=1,...,Na,N<sub>a</sub>为步骤1初始化得到的雷达系统方位向采样点数,PRF为步骤1初始化得到的雷达系统脉冲重复频率,<img file="FDA0000485737230000037.GIF" wi="43" he="63" />为步骤1初始化得到的雷达平台速度矢量,<img file="FDA0000485737230000038.GIF" wi="136" he="81" />为步骤1初始化得到的雷达平台主天线初始位置矢量,<img file="FDA0000485737230000039.GIF" wi="124" he="80" />为步骤1初始化得到的雷达平台副天线初始位置矢量;采用公式<maths num="0001" id="cmaths0001"><math><![CDATA[<mrow><mover><mi>P</mi><mo>&OverBar;</mo></mover><mrow><mo>(</mo><mi>a</mi><mo>,</mo><mi>r</mi><mo>)</mo></mrow><mo>=</mo><msub><mover><mi>P</mi><mo>&OverBar;</mo></mover><mi>c</mi></msub><mo>+</mo><mrow><mo>(</mo><mi>r</mi><mo>-</mo><mn>1</mn><mo>)</mo></mrow><mo>&CenterDot;</mo><msub><mi>d</mi><mi>r</mi></msub><mo>&CenterDot;</mo><msub><mover><mi>&zeta;</mi><mo>&OverBar;</mo></mover><mi>u</mi></msub><mo>+</mo><mrow><mo>(</mo><mi>a</mi><mo>-</mo><mn>1</mn><mo>)</mo></mrow><mo>&CenterDot;</mo><msub><mi>d</mi><mi>a</mi></msub><mo>&CenterDot;</mo><msub><mover><mi>&zeta;</mi><mo>&OverBar;</mo></mover><mi>v</mi></msub><mo>+</mo><mi>H</mi><mrow><mo>(</mo><mi>a</mi><mo>,</mo><mi>r</mi><mo>)</mo></mrow><mo>&CenterDot;</mo><msub><mover><mi>&zeta;</mi><mo>&OverBar;</mo></mover><mi>h</mi></msub></mrow>]]></math><img file="FDA00004857372300000310.GIF" wi="1323" he="76" /></maths>计算得到场景散射点P(a,r)的位置矢量,其中a表示散射点位于场景方位向的第a个位置,a=1,...s,<sub>a</sub>,s<sub>a</sub>为场景方位向总的散射点数,r表示散射点位于场景距离向的第r个位置,r=1,...,s<sub>r</sub>,s<sub>r</sub>为场景距离向总的散射点数,d<sub>r</sub>为步骤1初始化得到的场景距离向散射点间隔,d<sub>a</sub>为步骤1初始化得到的场景方位向散射点间隔,<img file="FDA00004857372300000311.GIF" wi="63" he="73" />为步骤1初始化得到的场景参考点位置矢量,H(a,r)为步骤1初始化得到的场景散射点P(a,r)的低精度高程,<img file="FDA00004857372300000312.GIF" wi="53" he="69" />为步骤1中定义的InSAR成像空间的第一个坐标基,<img file="FDA00004857372300000313.GIF" wi="45" he="69" />为步骤1中定义的InSAR成像空间的第二个坐标基,<img file="FDA00004857372300000314.GIF" wi="57" he="69" />为步骤1中定义的InSAR成像空间的第三个坐标基;i为雷达平台距离散射点P(a,r)前后半个合成孔径内的慢时刻i的取值满足条件:i=a‑round(N<sub>l</sub>/2),...,a+round(N<sub>l</sub>/2),若a‑round(N<sub>l</sub>/2)<1,则a‑round(N<sub>l</sub>/2)=1,若a+round(N<sub>l</sub>/2)>N<sub>a</sub>,则a+round(N<sub>l</sub>/2)=N<sub>a</sub>,其中round(·)为近似取整函数,N<sub>l</sub>为步骤1初始化得到的雷达平台一个合成孔径长度内的慢时刻个数,N<sub>a</sub>为步骤1初始化得到的雷达系统方位向采样点数,即i只取雷达系统采样范围内的雷达平台距离散射点P(a,r)前后半个合成孔径内的慢时刻;采用公式<img file="FDA0000485737230000041.GIF" wi="504" he="92" />计算得到雷达平台主天线i时刻处散射点P(a,r)的距离史<img file="FDA0000485737230000042.GIF" wi="210" he="72" />采用公式<img file="FDA0000485737230000043.GIF" wi="480" he="92" />计算得到雷达平台副天线i时刻处散射点P(a,r)的距离史<img file="FDA0000485737230000044.GIF" wi="193" he="72" />其中||·||<sub>2</sub>为L2范数;步骤4、距离压缩后数据插值重采样<img file="FDA0000485737230000045.GIF" wi="165" he="70" />为雷达平台主天线i时刻处散射点P(a,r)的距离史<img file="FDA0000485737230000046.GIF" wi="176" he="72" />对应的距离门,计算公式为<maths num="0002" id="cmaths0002"><math><![CDATA[<mrow><msub><mover><mi>I</mi><mo>&OverBar;</mo></mover><mi>m</mi></msub><mrow><mo>(</mo><mi>i</mi><mo>;</mo><mi>a</mi><mo>,</mo><mi>r</mi><mo>)</mo></mrow><mo>=</mo><msub><mi>I</mi><mi>c</mi></msub><mo>+</mo><mi>round</mi><mrow><mo>(</mo><mrow><mo>(</mo><msub><mover><mi>R</mi><mo>&OverBar;</mo></mover><mi>m</mi></msub><mrow><mo>(</mo><mi>i</mi><mo>;</mo><mi>a</mi><mo>,</mo><mi>r</mi><mo>)</mo></mrow><mo>-</mo><msub><mi>R</mi><mi>mc</mi></msub><mo>)</mo></mrow><mo>/</mo><msub><mi>d</mi><mi>r</mi></msub><mo>)</mo></mrow><mo>,</mo><msub><mover><mi>I</mi><mo>&OverBar;</mo></mover><mi>s</mi></msub><mrow><mo>(</mo><mi>i</mi><mo>;</mo><mi>a</mi><mo>,</mo><mi>r</mi><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000485737230000047.GIF" wi="980" he="72" /></maths>为雷达平台副天线i时刻处散射点P(a,r)的距离史<img file="FDA0000485737230000048.GIF" wi="176" he="72" />)对应的距离门,计算公式为<maths num="0003" id="cmaths0003"><math><![CDATA[<mrow><msub><mover><mi>I</mi><mo>&OverBar;</mo></mover><mi>s</mi></msub><mrow><mo>(</mo><mi>i</mi><mo>;</mo><mi>a</mi><mo>,</mo><mi>r</mi><mo>)</mo></mrow><mo>=</mo><msub><mi>I</mi><mi>c</mi></msub><mo>+</mo><mi>round</mi><mrow><mo>(</mo><mrow><mo>(</mo><msub><mover><mi>R</mi><mo>&OverBar;</mo></mover><mi>s</mi></msub><mrow><mo>(</mo><mi>i</mi><mo>;</mo><mi>a</mi><mo>,</mo><mi>r</mi><mo>)</mo></mrow><mo>-</mo><msub><mi>R</mi><mi>sc</mi></msub><mo>)</mo></mrow><mo>/</mo><msub><mi>d</mi><mi>r</mi></msub><mo>)</mo></mrow><mo>,</mo></mrow>]]></math><img file="FDA0000485737230000049.GIF" wi="775" he="72" /></maths>其中round(·)为近似取整函数,<img file="FDA00004857372300000410.GIF" wi="180" he="72" />为步骤3中得到的雷达平台主天线i时刻处散射点P(a,r)的距离史,<img file="FDA00004857372300000411.GIF" wi="179" he="72" />为步骤3中得到的雷达平台副天线i时刻处散射点P(a,r)的距离史,R<sub>mc</sub>为步骤1初始化得到的场景参考点到雷达平台主天线各慢时刻天线相位中心的最短距离,R<sub>sc</sub>为步骤1初始化得到的场景参考点到雷达平台副天线各慢时刻天线相位中心的最短距离,I<sub>c</sub>为步骤1初始化得到的场景参考点在雷达平台主天线回波数据和雷达平台副天线回波数据中的距离门位置,d<sub>r</sub>为步骤1初始化得到的场景距离向散射点间隔,i为步骤3中得到的雷达平台距离散射点P(a,r)前后半个合成孔径内的慢时刻,a表示散射点位于场景方位向的第a个位置,a=1,...,s<sub>a</sub>,s<sub>a</sub>为场景方位向总的散射点数,r表示散射点位于场景距离向的第r个位置,r=1,...,s<sub>r</sub>,s<sub>r</sub>为场景距离向总的散射点数;从步骤2中得到的平台主天线距离压缩后数据<img file="FDA00004857372300000412.GIF" wi="64" he="79" />的第i行中取第<img file="FDA00004857372300000413.GIF" wi="164" he="70" />个数据的前W<sub>0</sub>个数据和后W<sub>0</sub>个数据,采用标准辛格插值方法对这组数据进行插值,得到插值重采样后的数据<img file="FDA00004857372300000414.GIF" wi="200" he="71" />从步骤2中得到的平台副天线距离压缩后数据<img file="FDA00004857372300000415.GIF" wi="49" he="79" />的第i行中取第<img file="FDA00004857372300000416.GIF" wi="167" he="70" />个数据的前W<sub>0</sub>个数据和后W<sub>0</sub>个数据,采用标准辛格插值方法对这组数据进行插值,得到插值重采样后的数据<img file="FDA00004857372300000417.GIF" wi="195" he="71" />其中W<sub>0</sub>为标准辛格插值的半个窗长;步骤5、插值重采样后数据相干求和根据补偿相位因子的计算公式<maths num="0004" id="cmaths0004"><math><![CDATA[<mrow><mover><mi>F</mi><mo>&OverBar;</mo></mover><mrow><mo>(</mo><mi>i</mi><mo>;</mo><mi>a</mi><mo>,</mo><mi>r</mi><mo>)</mo></mrow><mo>=</mo><mi>exp</mi><mo>{</mo><mi>j</mi><mo>&CenterDot;</mo><mn>4</mn><mo>&CenterDot;</mo><mi>&pi;</mi><msub><mover><mi>R</mi><mo>&OverBar;</mo></mover><mi>m</mi></msub><mrow><mo>(</mo><mi>i</mi><mo>;</mo><mi>a</mi><mo>,</mo><mi>r</mi><mo>)</mo></mrow><mo>/</mo><mi>&lambda;</mi><mo>}</mo><mo>,</mo></mrow>]]></math><img file="FDA0000485737230000051.GIF" wi="791" he="98" /></maths>得到散射点P(a,r)在慢时刻i应补偿的相位因子<img file="FDA0000485737230000052.GIF" wi="208" he="85" />其中j为虚数单位(即‑1开根),<img file="FDA0000485737230000053.GIF" wi="214" he="76" />为步骤3中得到的雷达平台主天线i时刻处散射点P(a,r)的距离史,i为步骤3中得到的雷达平台距离散射点P(a,r)前后半个合成孔径内的慢时刻,λ为步骤1初始化得到的雷达系统工作的信号波长,a表示散射点位于场景方位向的第a个位置,a=1,...,s<sub>a</sub>,s<sub>a</sub>为场景方位向总的散射点数,r表示散射点位于场景距离向的第r个位置,r=1,...,s<sub>r</sub>,s<sub>r</sub>为场景距离向总的散射点数;将步骤4中得到的插值重采样后的数据<img file="FDA0000485737230000054.GIF" wi="170" he="76" />与散射点P(a,r)在慢时刻i应补偿的相位因子<img file="FDA0000485737230000055.GIF" wi="190" he="85" />相乘,得到相位补偿后的数据<img file="FDA0000485737230000056.GIF" wi="233" he="86" />将步骤4中得到的插值重采样后的数据<img file="FDA0000485737230000057.GIF" wi="164" he="72" />与散射点P(a,r)在慢时刻i应补偿的相位因子<img file="FDA0000485737230000058.GIF" wi="193" he="90" />相乘,得到相位补偿后的数据<img file="FDA0000485737230000059.GIF" wi="232" he="87" />把雷达平台主天线距离散射点P(a,r)前后半个合成孔径内的所有慢时刻的相位补偿后的数据<img file="FDA00004857372300000510.GIF" wi="204" he="86" />相加,得到雷达平台主天线场景散射点P(a,r)成像后的数据<img file="FDA00004857372300000511.GIF" wi="229" he="82" />即<img file="FDA00004857372300000512.GIF" wi="515" he="117" />把雷达平台副天线距离散射点P(a,r)前后半个合成孔径内的所有慢时刻的相位补偿后的数据<img file="FDA00004857372300000513.GIF" wi="194" he="88" />相加,得到雷达平台副天线场景散射点P(a,r)成像后的数据<img file="FDA00004857372300000514.GIF" wi="176" he="76" />即<maths num="0005" id="cmaths0005"><math><![CDATA[<mrow><msub><mover><mi>B</mi><mo>&OverBar;</mo></mover><mi>s</mi></msub><mrow><mo>(</mo><mi>a</mi><mo>,</mo><mi>r</mi><mo>)</mo></mrow><mo>=</mo><munder><mi>&Sigma;</mi><mi>i</mi></munder><msub><mover><mi>T</mi><mo>&OverBar;</mo></mover><mi>s</mi></msub><mrow><mo>(</mo><mi>i</mi><mo>;</mo><mi>a</mi><mo>,</mo><mi>r</mi><mo>)</mo></mrow><mo>;</mo></mrow>]]></math><img file="FDA00004857372300000515.GIF" wi="489" he="118" /></maths>步骤6、干涉相位提取将步骤5中得到的雷达平台主天线场景散射点P(a,r)成像后的数据<img file="FDA00004857372300000516.GIF" wi="148" he="75" />和雷达平台副天线场景散射点P(a,r)成像后的数据<img file="FDA00004857372300000517.GIF" wi="144" he="72" />进行共轭相乘,得到InSAR成像的干涉相位<img file="FDA00004857372300000518.GIF" wi="153" he="71" />完成基于曲面投影的InSAR成像处理,其中a表示散射点位于场景方位向的第a个位置,a=1,...,s<sub>a</sub>,s<sub>a</sub>为场景方位向总的散射点数,r表示散射点位于场景距离向的第r个位置,r=1,...,s<sub>r</sub>,s<sub>r</sub>为场景距离向总的散射点数。
地址 611731 四川省成都市高新区(西区)西源大道2006号