基于半牛顿迭代法的星载干涉合成孔径雷达基线估计方法

摘要

基于半牛顿迭代法的星载干涉合成孔径雷达基线估计方法，它不通过干涉复图像距离向频谱图来直接获得频率范围，而是求干涉复图像中每一像素点对应的距离向干涉条纹频率，一定距离向范围内的干涉复数据对应一定范围的频率，由此可得到所求的频率值，利用半牛顿迭代法可以精确地估计该频率范围值，然后结合方程组可以很精确地算出基线参数；采用本发明的方法基线长度的精确度可达到厘米级别，大大改善了利用干涉复图像估计基线的精度，从而极大地提高了最终数字高程的精度。

权利要求

1.基于半牛顿迭代法的星载干涉合成孔径雷达基线估计方法，其特征在于实现步骤如下：

第一步，对选定距离向范围内干涉复图像各像素点进行干涉条纹频率的粗估计，方法为：

(1.1)输入干涉复图像数据，对所述干涉复图像数据的幅度进行归一化处理，然后在干涉复图像数据上选取一定距离向范围的干涉复图像数据Ei，i＝1，2，...l，i为距离向的点数，选取的干涉复图像数据的距离向斜距值，即主星到地面点斜距值应大于需要的rmin，小于需要的rmax，rmin和rmax为最终需要的主星到地面点斜距值的最小值和最大值；

(1.2)在选取的距离向范围内的干涉复图像数据，除边界处干涉复图像数据外，对距离向第i，i＝N+1，N+2，...l-N个干涉复图像数据距离向前后各选取N个点，N为自然数，一共2N+1个点作为模板，设模板内数据点值为yim，m＝-N，…，N，yim∈Ei，求该选取的距离向范围内的干涉复图像各像素点处的条纹瞬时频率；边界处干涉复图像数据依情况选取点数作为模板，最后边界处的干涉复图像数据频率估值将被舍弃；

(1.3)在模板中的2N+1个点后补零至N1个点，N1＞2N+1，然后对补零后的数据进行距离向的快速傅里叶变换至频域，设为Y(n)，n＝1，2，…N1，n为频域的点数，求出频谱峰值点对应的条纹频率值，设在n0处频谱幅度最大，n0为频谱峰值对应的点数；

(1.4)在区间[n0-2，n0+2]上，再对频谱进行拉格朗日32倍插值，求得插值后频谱幅度最大值对应的条纹频率值，设为k′i，i＝1，2，...l，即获得干涉复图像上选定的距离向范围内第i个像素点对应的瞬时频率粗估计值k′i；

第二步，将瞬时频率粗估计值k′i进行半牛顿迭代计算，直至得到精确的条纹频率范围值，实现步骤为：(2.1)利用公式(16)，将瞬时频率粗估计值k′i作为迭代的初始值，并将H(k)代入公式(16)中，可得到利用半牛顿法一次迭代后的精确瞬时频率估计值ki，其中

Re{·}表示求实部的运算，k为一变量，表示条纹频率，H(k)代表关于变量k的函数，T表示距离向采样间隔，m为自然数，-N≤m≤N，k′i代表进行迭代前的频率粗估计值，H′k(k′i)表示H(k)对k求一次导数后的函数在k′i点处的取值，H″k(k′i)表示H(k)对k求两次导数后的函数在k′i点处的取值；

(2.2)将一次迭代后得到的ki作为第二次迭代的k′i值，再利用(16)式，可得到二次迭代后ki的值，依此规律迭代直至频率不再变化为止，此时获得比较精确的条纹频率值ki；

(2.3)经过半牛顿迭代后总体上干涉复图像条纹频率k的值随着距离向斜距值r取值的变大而减小，当基线较长时，频率估计值会有些微小的起伏波动，选取需要的斜距范围[rmin，rmax]内的频率值；为增加精确度，可以再进行曲线线性拟合，求得最终精确的频率范围值，经过拟合后频率k的值随着斜距r取值单调减小，得到kmin为干涉复图像距离向斜距rmax处的条纹频率值，kmax为干涉复图像在距离向斜距rmin处的条纹频率值，kmin和kmax分别为干涉复图像在距离向范围[rmin，rmax]内条纹频率的最小值和最大值；

第三步，将步骤(2.3)中求得的kmin和kmax处的值代入式(10)中，即

其中

r为对应于从干涉条纹图中选取的像元距离向的斜距值，即对应主星到地面点的斜距值；p(r)，g(r)为关于r的两个函数，H表示卫星到地面的高度；λ为信号波长，解得基线参数分量Bx和By，Bx＝B·cosα，By＝B·sinα，Bx和By分别为基线的水平分量和垂直分量；

则基线长度B和基线倾角α的值可通过式(11)和式(12)得到，

对于收发分置星载干涉合成孔径雷达系统，所述d取为0.5，对于收发同置星载干涉合成孔径雷达系统，所述d取为1。

2.根据权利要求1所述的基于半牛顿迭代法的星载干涉合成孔径雷达基线估计方法，其特征在于：第二步中需要迭代28-32次。

说明书

基于半牛顿迭代法的星载干涉合成孔径雷达基线估计方法

技术领域

[0001] 本发明涉及一种星载干涉合成孔径雷达(InSAR)基线的估计方法，属于信号处理技术领域。

背景技术

[0002] 星载干涉合成孔径雷达(InSAR)是在合成孔径雷达基础上发展起来的，以合成孔径雷达复数据提取的相位信息为信息源来获取地表的三维信息和变化信息。具有全天时、全天候获取大面积数字高程图的特点。InSAR利用两部天线获取同一区域的复雷达图像对，这两幅图像经过精确配准和复共轭相乘后得到干涉图像，干涉图像中的相位差值即为两次成像的相位差测量值，根据两次成像相位差与地面目标的三维空间位置之间存在的几何关系，利用飞行轨道的参数如辐射波长，天线视角，干涉基线，天线高度等参数既可测定地面目标的三维坐标，它可以用来提供大范围的高精度数字高程模型(DEM)。

[0003] 基线是InSAR中引入的重要概念，它定义为两个雷达传感器或两部天线之间的矢量几何关系。作为干涉合成孔径雷达的基础，基线决定着系统的干涉性能，它是InSAR计算地表高程或形变公式中的重要参数，其精度直接影响到最终提取的地表高程或形变的精度。因此精确的基线估计是极为重要和必需的。

[0004] 目前基线参数的估计主要有三种方式：利用卫星轨道状态矢量的基线估计；利用控制点来估计基线的长度；利用干涉相位信息的基线估计。

[0005] 第一种方式从基线的空间几何关系出发，根据星历表参数和雷达系统本身的参数来获得基线分量。但由于大气阻力和太阳万有引力的影响，卫星轨道常发生漂移，轨道参数具有一定程度的不确定性，卫星定轨精度不够高。由仪器测得的数据本身也会由于温度等影响存在一定的误差，使得该方法确定的基线参数精度较差，会给高程带来较大误差。第二种方式依赖几个地面控制点的高度信息和干涉相位并联立方程组解得基线参数，且基线精度依赖于控制点的精度。该方法对地面控制点的选取有一定的要求。而实际中控制点的布局及测量精度都很难满足要求。第三种方式直接利用已知的干涉相位数据获得基线参数，不需要轨道参数和设置控制点，因此更得到人们的青睐。目前基于干涉相位信息的基线估计方法主要有利用解缠相位进行估计和利用缠绕干涉相位进行估计。前者依赖于高精度的绝对干涉相位值，需要进行相位展开，而相位展开过程本身会带来较大的误差值。后者由Kuldip Singh等人于1997年第一次提出，直接利用从SAR干涉测量数据中提取的干涉条纹来估测基线。是目前基线估计方法里面最简洁的方法。该方法不需要对干涉相位进行解缠处理，仅仅对一定距离向范围内的干涉复图像数据在距离向进行傅里叶变换，得到干涉复图像距离向的频谱，由频谱图直接获得相应的条纹频率范围，并将该频率范围代入方程式解方程得到基线长度和基线倾角。该方法的关键是求得一定距离向范围内干涉复图像的条纹频率范围。对于星载InSAR，为最终得到精确的基线参数，需要有精度非常高的频率范围值，目前对于频率范围的选取主要是利用频谱图来获取频率值，因此估计精度很差。还没有找到合适的方法得到高精度的频率范围值，这也使得最终求得的基线值误差很大。利用干涉相位图的缠绕相位进行基线估计的方法因此受限。

发明内容

[0006] 本发明的技术解决方案：克服现有技术的不足，提供一种基于半牛顿迭代法的星载干涉合成孔径雷达基线估计方法，该方法使基线长度的精确度可达到厘米级别，大大改善了利用干涉相位估计基线的精度，极大地提高了最终数字高程的精度。

[0007] 本发明的技术解决方案：一种基于半牛顿迭代法的星载干涉合成孔径雷达基线估计方法，是利用半牛顿迭代法精确估计星载InSAR基线参数值，不需要对干涉相位进行解缠处理，直接通过对干涉复图像数据处理获得基线信息，从而实现精确估计地面目标的数字高程值。

[0008] 本发明的原理为：首先说明利用干涉相位信息估计基线的原理。由图1中InSAR基线几何图，可得

[0009] 对于干涉SAR系统，

[0010]

[0011] 对于收发分置星载InSAR系统，d取值为0.5，收发同置星载InSAR系统，d取值为1。φ1，φ2分别为主副天线由于往返路程接收的回波信号相位。φ为相位差。由几何关系，可得

[0012]

[0013] 若地面平坦，则有H＝rcosθ，

[0014]

[0015] 实际上测得的相位差是模糊相位值，φwrapped＝φ-2kπ。要得到真实相位差φ，需要对模糊相位进行解缠处理，解缠过程本身将带来一定的误差。本发明可以跳过解缠这一步，直接利用干涉相位值进行基线估计。对相位进行Fourier变换得：

[0016]

[0017] 选择一定距离向范围的干涉复图像数据，rmin和rmax对应于从干涉复图像中选取的边缘像元的斜距值。此时r∈[rmin，rmax]，ψ(k)为选定距离向干涉复图像的频谱，利用驻定相位原理进行近似处理，可得到(4)式中频谱ψ(k)是带限的。即干涉复图像条纹频率k存在于一定范围内。从图3中可以算出频率k的取值范围，但精确度非常低。同时由驻定相位原理，可得

[0018]

[0019] 由于

[0020] 代入(5)式，得到

[0021]

[0022] 可以看出干涉条纹频率k的值随着r取值单调线性减小。

[0023] 令Bx＝B·sinθ0＝B·cosα (8a)

[0024] By＝-B·cosθ0＝B·sinα (8b)

[0025]

[0026]

[0027] 代入(5)式则可得到方程组

[0028]

[0029] 将从频谱图中读出的kmin和kmax的值代入(10)式，可解得基线参数分量Bx和By。则基线长度B和基线倾角α的值可通过式(11)和式(12)得到。

[0030]

[0031]

[0032] 但由频谱图获取的kmin和kmax的值精度很低，使求得的基线值误差很大。

[0033] 本发明引入了半牛顿迭代法，可以精确求取kmin和kmax的值。由干涉相位信息估计基线的原理可知kmin和kmax的求取是基线估计的关键所在。由于k的值随着r取值单调线性减小，可以从选定距离向范围内干涉复图像数据的瞬时频率着手。

[0034] 本发明的实现步骤如下：

[0035] 第一步，对选定距离向范围内干涉复图像各像素点进行条纹频率的粗估计，方法为：

[0036] (1.1)输入干涉复图像数据，对所述干涉复图像数据的幅度进行归一化处理，然后在干涉复图像数据上选取一定距离向范围的复图像数据Ei，i＝1，2，...l，i为距离向的点数，选取的复图像数据的距离向斜距值，即主星到地面点斜距值应小于需要的rmin，大于需要的rmax，rmin和rmax为最终需要的主星到地面点斜距值的最小值和最大值；

[0037] (1.2)在选取的距离向范围内的复图像数据，除边界处复图像数据外，对距离向第i，i＝N+1，N+2，...l-N个复图像数据距离向前后各选取N个点，N为自然数，一共2N+1个点作为模板，设模板内像素点值为yim，m＝-N，…，N，yim∈Ei，求该复图像数据点的瞬时频率；边界处复图像数据依情况选取点数作为模板，最后边界处的复图像数据频率估值将被舍弃；

[0038] (1.3)为增大频域的分辨率，提高求取频率的精度，在模板中的2N+1个点后补零至N1个点，N1＞2N+1，然后对补零后的数据进行距离向的快速傅里叶变化(FFT)至频域设为Y(n)，求频谱峰值点对应的条纹频率值，设在n0处为频谱幅度最大，n0为频谱峰值点对应的点数；

[0039] (1.4)为更精确估计频谱峰值点，在区间[n0-2，n0+2]上，再对频谱进行拉格朗日32倍插值，计算求得插值后频谱幅度最大值对应的频率，设为k′i，i＝1，2，...l，即获得干涉复图像上选定的距离向范围内第i个像素点对应的瞬时频率粗估计值k′i。该步骤中获得的瞬时频率估计值k′i误差是比较大的，要精确地得到像素点的瞬时频率ki还需要进一步计算；

[0040] 第二步，将粗估计时得到的条纹瞬时频率值k′i进行半牛顿迭代计算，直至得到精确的频率范围值，步骤为：

[0041] 由于模板内像素点值可由下式表示： m＝-N，…，N，则求瞬时频率的运算可以表示为求

[0042]

[0043]

[0044] Re{·}表示求实部的运算。当k取ki时G(f)最大，(14)的运算等效于求的零点，可用牛顿迭代法来实现，即求

[0045]

[0046] 但由于 是未知的，所以引入半牛顿迭代法求零点

[0047]

[0048] 其中

[0049] (2.1)利用(16)式，将粗估计时得到的k′i作为迭代的初始值，并将H(k)代入该式中，可得到利用半牛顿法一次迭代后的ki的值；Re{·}表示求实部的运算，k为一变量，表示条纹频率，H(k)代表关于变量k的函数，T表示距离向采样间隔，m为自然数，-N≤m≤N，k′i代表进行迭代前的频率粗估计值，H′k(k′i)表示H(k)对k求一次导数后的函数在k′i点处的取值，H″k(k′i)表示H(k)对k求两次导数后的函数在k′i点处的取值；

[0050] (2.2)将一次迭代后得到的ki作为第二次迭代的k′i值，再利用(16)式，得到二次迭代后的fi的值。为增加精确度，一直迭代直至频率值不再变化为止，一般迭代28-32次左右即可。

[0051] (2.3)经过半牛顿迭代后得到的频率估计值为ki，ki为选定距离向范围内第i个复图像像素点对应的条纹频率值。总体上干涉条纹频率k的值随着距离向斜距r取值的变大而减小，当基线较长时，频率估计值会有些微小的起伏波动，选取需要的斜距范围[rmin，rmax]内的频率值；为增加精确度，再进行曲线线性拟合，求得最终精确的条纹频率范围值。经过拟合后频率k的值随着斜距r取值单调减小，得到kmin为干涉条纹在距离向rmax处的条纹频率值，kmax为干涉条纹在距离向rmin处的条纹频率值，kmin和kmax分别为干涉条纹在距离向范围[rmin，rmax]内条纹频率的最小值和最大值；

[0052] 第三步，将步骤(2.3)中求得的kmin和kmax处的值代入式(10)中，即[0053]

[0054] 其中

[0055]

[0056] r为对应于从干涉条纹图中选取的像元距离向的斜距值，即对应主星到地面点的斜距值；p(r)，g(r)为关于r的两个函数，H表示卫星到地面的高度，θ为下视角；λ为信号波长，解得基线参数分量Bx和By，Bx＝B·cosα，By＝B·sinα，Bx和By分别为基线的水平分量和垂直分量；

[0057] 则基线长度B和基线倾角α的值可通过式(11)和式(12)得到，

[0058]

[0059]

[0060] 本发明与现有技术相比的优点在于：

[0061] (1)本发明利用半牛顿迭代法来估计星载InSAR基线参数，基线长度的精确度可达到厘米级别(2-5厘米)，大大提高了基线估计的精确度，弥补了传统的利用干涉相位信息估计基线误差太大的缺点。实现了直接利用星载干涉复图像精确估计基线值。为利用干涉复图像进行基线估计的方法开辟了一条新思路。

[0062] (2)利用半牛顿迭代法估计基线参数时不需要对干涉相位进行解缠处理，直接利用干涉相位复图像。免去了解缠过程中各种误差的引入和解缠过程带来相位估计值精度下降的缺点。

[0063] (3)利用半牛顿迭代法估计基线参数不需要依赖轨道参数和设置控制点，只需求得干涉相位的精确瞬时频率值，再求解二元一次方程，方便简单。

附图说明

[0064] 图1所示为干涉SAR系统几何关系图；

[0065] 图2为本发明方法的实现流程图；

[0066] 图3给出了在具体参数下求得的干涉相位频谱ψ(k)的幅度图；

[0067] 图4为本发明在选定的距离向范围内像素点频率插值后的条纹频率粗估计值；

[0068] 图5为本发明选定的[rmin，rmax]范围内像素点经半牛顿迭代法第一次迭代后的条纹频率；

[0069] 图6为本发明在[rmin，rmax]范围内的像素点经半牛顿迭代法后得到的最后条纹频率值；

[0070] 图7为本发明对半牛顿迭代法得到的条纹频率值再经过线性拟合后得到的最终结果。

[0071] 图中符号表示如下：图1中A1，A2分别表示卫星上的两部雷达传感器，A1为主星，B表示基线的长度，r，R2分别表示两传感器到地面点p的斜距。rmin和rmax为最终需要的距离向斜距值(主星到地面点斜距值)的最小值和最大值，即rmin和rmax分别对应到地面不同点的r值，H表示主星到地面的高度。α为基线倾角，θ为下视角，h为点p到地面的高度。

具体实施方式

[0072] 下面利用由具体的收发分置星载InSAR参数得到的干涉复图像数据来验证半牛顿迭代法估计基线的有效性。

[0073] 表1给出了主要的星载InSAR参数：

[0074] 表1

[0075]InSAR参数 数值
主星高度H 514km
信号波长λ 0.031m
脉冲重复频率 150MHZ
下视角θ 40°
基线长度B 200m
基线倾角α 45°
距离向像素间隔 1.0m
波束中心处斜距R0 670979.346m

[0076] 尽量选取地面平坦地区，利用星载InSAR参数可获得相应的干涉复图像数据。

[0077] 在复图像上选取一定距离向范围的复数据，令斜距rmin为670487.346米，rmax为671470.346米，rmin和rmax之间共有984个数据点。由于选用的为收发分置InSAR系统，则在式(10)中d为0.5，将具体参数代入(9a)和(9b)中，需求解的方程组(10)变为：

[0078]

[0079] 为求得基线参数Bx和By的值，需要知道kmin和kmax的精确取值。

[0080] 若利用传统的方法估计kmin和kmax，kmin和kmax为干涉复图像上一定距离向范围内波数域频率的最小值和最大值，即直接利用频谱图求取kmin和kmax。为了增大频域的分辨率，对选取的rmin和rmax之间984个数据点后补零到16384点，再对该选取的复数据在距离向进行傅里叶变换，可以得到选定距离向范围内干涉复图像的频谱ψ(k)，ψ(k)的幅度谱如图3所示。可看到功率谱是带限的，且近似为一尖峰脉冲，kmin和kmax的值相差很小，从功率谱图像上很难求出精确的频率范围值。可看到频谱的最高值点在频率点0.07191处，对带宽的选取准则很难把握，若选用幅值降为一半时的频率点为带宽点，则可算得传统方法下kmin＝0.0692256，kmax＝0.0745944，代入式(18)中可解得Bx＝5275.04，By＝5975.84，从而得基线长度B＝7970.99米，与理论值200米相比，误差太大。利用缠绕干涉相位进行基线估计方法对频率精确度要求非常高，直接从频谱图上读出的频率范围值远远不能达到需要的精度。

[0081] 如图2所示，本发明对选定的复图像利用半牛顿迭代法估计基线时步骤如下：

[0082] 步骤(1.1)：输入干涉复图像数据，并对复数据幅度进行归一化，令各像素点幅值为1，考虑到估计瞬时频率时边界处数据不准的情况，在复图像上选取距离向范围670467.346米到671490.346米的复图像数据Ei，i＝1，2，...1024，而实际需要的rmin为
670487.346米，rmax为671470.346米，且rmin和rmax之间共有984个数据点。

[0083] 步骤(1.2)：在选取的数据范围内，除边界处数据，对距离向第i，i＝6，7，...1019个复图像数据前后各选取5个点，一共11个点作为模板设为yim，m＝-5，…，5，ym∈Ei，并求该点的条纹瞬时频率。边界处依情况选取点数，最后边界处的频率估值将被舍弃，不影响最后kmin和kmax的取值。

[0084] 步骤(1.3)：为增大频域的分辨率，对模板中的11个点后补零至16384个点，此处N1＝16384。然后对补零后的数据进行距离向的快速傅里叶变换(FFT)至频域设为Y(n)。求出频谱峰值点对应的条纹频率值，设在频域第n0点处为频谱幅度最大。

[0085] 步骤(1.4)：在频谱区间[n0-2，n0+2]上对频谱进行拉格朗日32倍插值，n0为频谱峰值处对应的点数，求得插值后频谱幅度最大值对应的条纹频率，设为k′i，i＝1，2，...1024。即获得干涉复图像上选定的距离向范围内第i个像素点对应的瞬时频率粗估计值。其取值见图4所示。横坐标为距离向上对应的点数。

[0086] 步骤(2.1)：将粗估计时得到的k′i作为迭代的初始值，并将H(k)代入式(16)中，可得到利用半牛顿法一次迭代后的ki的值，如图5所示。

[0087] 步骤(2.2)：将一次迭代后得到的ki作为第二次迭代的k′i值，再利用式(16)，可得到二次迭代后的fi的值。依此规律一直迭代至频率值不再变化为止，迭代30次后得到精确的条纹频率值ki，如图6所示。

[0088] 步骤(2.3)：由图6可看到经过半牛顿迭代后总体上干涉条纹频率k的值随着距离向斜距r取值减小，但有些微小的起伏波动，选取需要的[rmin，rmax]，rmin和rmax为需要的主星到地面点斜距值的最小值和最大值，在该范围内有984个像素点。为增加精确度，再进行曲线线性拟合，如图7所示。可得到最终精确度很高的条纹频率范围值。

[0089] 步骤3：由图7的数据可以得到在[rmin，rmax]范围内的第1个像素点，即rmin＝670487.346米处的条纹频率值为kmax＝0.0715580，在第984个像素点，即rmax＝671470.346米处的条纹频率值为kmin＝0.0719065。将求得的kmin和kmax处的值代入式(18)中，并解方程组，可得到基线参数分量Bx和By的值分别为Bx＝141.433m，By＝141.436m，由式(11)和式(12)，可求得基线长度B＝200.0191米和基线倾角α＝45.0004°，与理论基线长度仅相差2厘米左右，基线倾角α也仅与理论值相差0.0004°。基线参数误差很小。相比传统方法求得的B＝7970.99米，利用半牛顿迭代法进行星载InSAR基线估计时误差非常小。也证明了该方法的实用性及可行性。

[0090] 总之，本发明引入了一种新的方法-半牛顿迭代法来估计干涉复图像距离向的条纹频率范围，从而估计基线的参数值。它不通过干涉复图像距离向频谱图来直接获得频率范围，而是求干涉复图像中每一点对应的距离向干涉条纹频率，一定距离向范围内的干涉复数据对应一定范围的频率，由此可得到所求复图像像素点的频率值。利用半牛顿迭代法可以精确地估计该频率范围值。然后结合方程组可以很精确地算出基线参数。基线长度的精确度可达到厘米级别。大大改善了利用干涉相位估计基线的精度，从而极大地提高了最终数字高程的精度。

[0091] 本发明未详细阐述部分属于本领域公知技术。