基于粗数字高程图的InSAR绝对相位模糊估计方法与流程

本发明属于信号处理
技术领域：
，特别涉及一种干涉合成孔径雷达InSAR的绝对相位模糊估计方法，可用于以获取有高程起伏地面场景的数字高程图为目的的雷达实时信号处理系统。
背景技术：
：干涉合成孔径雷达InSAR的应用极其广泛，最直接的应用就是获取数字高程图数字高程图。单航过InSAR获取地面场景的数字高程图依靠的是两个天线的雷达波束的波程差，这波程差反映在两个天线雷达回波的相位上。从雷达回波中解出整幅场景的绝对相位是获取整幅场景数字高程图的前提条件。在无地面精确控制点时，研究实时高效的绝对相位模糊估计方法对实时获取地面场景的数字高程图至关重要，也是InSAR实时信号处理技术的一个研究热点。缠绕相位是对两幅天线的SAR图像进行干涉处理得到的，是折叠在(-π,π]区间内的。解缠相位是对缠绕相位进行相位展开得到的，是一个连续的曲面。解缠后的干涉相位与其绝对干涉相位之间仍然相差2kπ，k为整数，这个k就称作绝对相位模糊数。基于粗数字高程图的InSAR绝对相位模糊估计是对于有高程起伏的场景，在无地面精确控制点时，利用粗数字高程图来解算出整幅地面场景的绝对相位模糊数k，结合解缠后的连续相位，得到整幅地面场景的绝对相位，进而反演出整幅地面场景的数字高程图。现有的绝对相位模糊估计方法大都是基于地面的精确控制点来计算模糊数k的。通常控制点以地固系下的三维坐标形式给出，并且控制点足够精确，将控制点反定位出图像中对应的坐标，确定出方位时刻，进而确定出绝对相位，再通过对应坐标位置的解缠相位确定出整幅地面场景的绝对相位模糊数k。这种方法的使用具有局限性，不适合于实时信号处理系统。陈立福等在文章“一种用于机载双天线InSAR系统的实时DEM生成方法及实现”(中国科学院大学学报,2014,31(5):678-683)中采用基于外部粗精度数字高程图数据的初始相位偏置实时估计算法。此种方法基于控制点和定标场，主要用于机载平台，以及较平坦的地形场景。而且上述方法进行实时估计时需要进行天线下视角的反演，需要进行反余弦计算，且反余弦的结果应该在[-π/2,π/2]之间，因而上述方法的不足之处在于：反余弦计算的正确结果需要进行判断获取，否则反演时会出现问题，而且此种方法不适用于运动状态变化剧烈的雷达平台，以及场景高度起伏较大的情况。技术实现要素：本发明的目的在于提出一种基于粗数字高程图的InSAR绝对相位模糊估计方法，以在无地面精确控制点且高度起伏较大的地面场景下正确估计出场景的绝对相位模糊数。为实现上述目的，本发明的技术方案包括如下：(1)根据粗数字高程图确定场景的参考平均高度h；(2)获取主天线回波场景中心参考点斜距R1：(2a)根据雷达发射时间与接收波门时间的关系，结合场景的参考平均高度，估计出场景中心所在的采样单元i：i＝round[(R1′-Rn)/Rs]<1>其中，为中间变量，β是主天线的下视角，H是主天线平台高度，h是参考平均高度；是主天线回波的最近斜距，是距离向采样间隔，t0是波门打开时间相对雷达发射时间的延时，c是光速，fs是距离向采样频率；round表示四舍五入取整操作；(2b)由(2a)所得参数估计主天线回波场景中心斜距R1：R1＝Rn+i·Rs；<2>(3)将整幅场景解缠之后的相位以及对应解缠相位的解缠标志矩阵作为输入数据，计算场景中心参考点的解缠相位(4)根据InSAR成像的几何关系，计算场景中心参考点的绝对相位ψ：ψ=2Qπ(R1-R2)λ---<3>]]>其中是辅天线回波的场景中心斜距，为中间变量，B是基线长度，α是基线倾角，R1是主天线回波场景中心斜距，β是主天线的下视角；Q为系数，一发双收模式时Q＝1，自发自收模式时Q＝2，λ是载波波长；(5)根据场景中心参考点的解缠相位和绝对相位ψ，计算整幅场景的模糊数k：其中，round表示四舍五入取整操作。(6)利用整幅场景的绝对相位模糊数k进行后续的干涉相位真实值求解及场景数字高程图反演。本发明与现有的技术相比具有以下优点：第一，本发明中利用了雷达发射机打开时间与接收波门打开时间之间的关系，以及场景的参考平均高度来估计主天线场景回波的中心斜距，克服了脉冲重复频率过高引起的距离模糊问题。第二，本发明由于进行绝对相位模糊数估计的步骤比较简明，没有大量的运算，不需要大量的数据，能够提高雷达平台实时处理效率。第三，本发明中整幅场景的绝对相位模糊数是通过估计的主天线场景回波中心斜距得到的场景中心绝对相位，结合场景中心的解缠相位估计得到，不需要定标场对地面控制点进行精确位置定标，简化了计算流程。附图说明图1为本发明的实现流程图；图2为本发明的几何原理示意图；图3为本发明仿真时输入的解缠相位图；图4为本发明仿真时输入的解缠标志矩阵图；图5为用本发明生成仿真环境下回波数据所使用的基准数字高程图；图6为本发明用图5描述的基准数字高程图的回波数据所生成的实时数字高程图；图7为用本发明生成半实物仿真环境下回波数据所使用的基准数字高程图；图8为本发明用图7描述的基准数字高程图的回波数据所生成的实时数字高程图；具体实施方式下面结合附图对本发明做进一步的描述。参照附图1，本发明在基于现场可编程门阵列FPGA与多核心数字信号处理板TIDSPTMS320C6678的实时系统上实现，其具体实施步骤如下：步骤1，根据仿真回波用的基准数字高程图估算成像场景的参考平均高度h。参照图5和图7，所述的基准数字高程图为回波录取场景的数字高程图，根据该基准数字高程图估算成像场景的参考平均高度h，通过下式进行：h＝mean2(HDEM)；<5>其中，HDEM是基准数字高程图所对应的高度信息矩阵，mean2为2维求平均操作。步骤2，估计场景中心参考点P所在的采样单元i。2a)将波门打开时间相对雷达发射时间的延时t0与光速c相乘，所得商值除以2，得到主天线回波场景的最近斜距Rn。Rn=ct02;---<6>]]>2b)参照附图2，将雷达平台高度H与场景参考平均高度h作差，所得差值除以主天线下视角β的余弦值，得到中间变量R1′：R1=H-hcosβ;---<7>]]>2c)将光速c除以距离向采样频率fs，所得商除以2得到距离向采样间隔Rs：Rs=c2fs;---<8>]]>2d)由(2a)-(2c)所得结果，结合下式估计出场景中心参考点P所在的采样单元i，i＝round[(R1′-Rn)/Rs]；<9>2e)将计算好的i值作为常量固化在DSP的内存中，通过define或者int定义i的值，本实例将i在数字信号处理芯片DSP中设为常量。步骤3，估计主天线回波场景的中心斜距R1，并将其固化在FPGA中。根据由步骤2所得的参数，估计得到主天线回波场景中心参考点斜距R1；R1＝Rn+i·Rs；<10>将估计得到的R1在FPGA中进行固化，即通过硬件编程语言中assign函数、parameter函数或者非阻塞赋值方式赋给其常数值R1。步骤4，输入通道数据。4a)分别将解缠相位与解缠标志矩阵通过代码调试软件CCS导入到数字信号处理芯片DSP的第3代双倍速率动态随机存储器DDR3中；4b)FPGA通过1553B总线接口接收惯导参数信息，即雷达主天线的航向位置、天向位置、侧向位置、航天速度、天向速度、侧向速度、载体俯仰角、偏航角、横滚角。4c)将除4b)参数以外的下视角、基线倾角、波长这些系统参数固化在FPGA中。步骤5，将FPGA中参数发送给TIDSPTMS320C6678。FPGA将接收的1553B总线发送过来的惯导参数与其它的系统参数进行组合，通过RapidIO接口发送给多核心数字信号处理板TIDSPTMS320C6678，该TIDSPTMS320C6678将接收的参数存入共享内存中。步骤6，根据i值选取解缠相位和解缠标志矩阵，并将选取的解缠相位和解缠标志矩阵从第3代双倍速率动态随机存储器DDR3中搬移到核0的第二层快速存储器L2中。现有技术在InSAR数据处理时，是利用增强型直接内存访问数据传输模块EDMA3将InSAR数据搬移到每个核的L2中进行处理，但由于L2内存空间有限，数据处理时需要将中间处理结果实时存放到DDR3，这样会带来额外的运算时间。本发明绝对相位模糊数估计的步骤比较简明，没有大量的运算，不需要大量的数据，所以为了减少运算时间，本发明直接在核0的L2中进行数据操作，不再将中间处理结果存放到DDR3，即利用EDMA3将DDR3中的所有距离单元为i的解缠相位和解缠标志矩阵搬移到核0的L2中。步骤7，确定场景中心已解缠的参考点P位置。根据步骤6中搬入到核0的L2中的解缠标志矩阵，沿方位向寻找已解缠的点，从而确定参考点的方位向位置Na0。具体操作步骤为：7a)设置标志信号列矩阵F，初始化F为全零矩阵大小为1×N，N为方位单元长度，以表明此时未找到已解缠点；初始化中间变量A＝0；7b)利用for循环对L2中解缠标志矩阵沿方位向寻找，循环变量num值随着for循环循环次数的增加由1按照步长1增加到N；7c)由7b)设置完循环变量num后，进入for循环体，在for循环体中判断F[num]和解缠标志矩阵标志位的值，如果F[num]等于0，并且解缠标志矩阵的标志位大于2，则此时num的值便是已解缠点的方位向位置，同时用当前的循环变量num对中间变量A进行更新即A＝A+num，然后将F[num]赋值为1，表明该点为已解缠点；for循环体经N次循环之后，终止for循环，得到更新后的最终的中间变量A′和标志矩阵F′。7d)由7c)得到的最终标志矩阵F′、最终的中间变量A′，求出参考点方位向位置Na0：Na0=round(A′sum(F′));---<11>]]>sum为矩阵求和操作，round为四舍五入操作；步骤8，估计场景中心参考点P的绝对相位ψ。8a)参照附图2，通过RapidIO接收的主天线回波场景中心参考点斜距R1，在核0的第二层快速存储器L2中计算出场景中心参考点的绝对相位ψ：ψ=2Qπ(R1-R2)λ;---<12>]]>其中是辅天线回波的回波场景的中心斜距，为中间变量，B是基线长度，α是基线倾角，R1是主天线回波场景中心斜距，β是主天线的下视角；Q为系数，对于一发双收模式，其Q＝1，对于自发自收模式，其Q＝2，λ是载波波长；8b)将上述绝对相位ψ存于第二层快速存储器L2的地址(*(LL2Input+4096))中。步骤9，取出场景中心参考点P的解缠相位根据步骤7得到的参考点的方位向位置Na0和步骤2得到的距离采样单元i，在解缠相位中取出参考点的解缠相位并存于多核心数字信号处理板TIDSPTMS320C6678中核0的第二层快速存储器L2的地址(*(LL2Input+4097))中。步骤10，估计场景的绝对相位模糊数k。根据步骤8得到的场景中心参考点P的绝对相位ψ，步骤9得到的场景中心参考点P的解缠相位估计场景的绝对相位模糊数k。round为四舍五入取整操作，在数字信号处理芯片DSP中实现round四舍五入取整操作是利用floor向下取整操作完成；步骤11，计算绝对模糊相位差并存于共享内存中。11a)根据步骤10得到的绝对相位模糊数k，在多核心数字信号处理板TIDSPTMS320C6678中计算绝对模糊相位差：其中，为绝对模糊相位差；11b)将11a)得到的存储到共享内存中。步骤12，生成实时数字高程图。根据步骤11的共享内存中的绝对模糊相位差按照现有数字高程图生成方法生成实时数字高程图，即由图5生成的回波数据通过本发明处理得到的实时数字高程图，如图6所示；由图7生成的回波数据通过本发明处理得到的实时数字高程图，如图8所示。下面结合仿真实验对本发明的效果做进一步的说明。1、仿真条件：表1仿真参数参数符号数值单位基线长度B0.3m下视角Β22.8°脉冲宽度τ10μs波长λ0.008m飞行高度H32000m脉冲重复周期PRT35.2μs航向速度V3500m/s基线倾角α22.8°参照表1所示的仿真参数设计，并将通过1553B总线接收的惯导参数在MATLAB中按照要求写成.dat文件，由个人计算机PC通过以太网接口缓存到现场可编程门阵列FPGA的第3代双倍速率动态随机存储器DDR3中。解缠相位和解缠标志矩阵如附图3和附图4，通过代码调试软件CCS导入到数字信号处理芯片DSP的第3代双倍速率动态随机存储器DDR3中。由上述仿真条件分别进行仿真环境和半实物仿真环境下的实时数字高程图反演：仿真实验1：在仿真环境下对图5进行实时数字高程图反演，结果如图6，其中图6的横坐标是方位向像素单元，纵坐标是距离向像素单元，右侧的色度条表示图像显示的大小范围，即场景高度变化范围对比图5与图6，可以看出：基准数字高程图与实时数字高程图的右侧色度条显示范围是一致的，说明实时数字高程图得到的高程值是正确的。从图像上看，基准数字高程图与实时数字高程图基本一致，说明绝对模糊估计的正确性。对图5所示的基准数字高程图与图6所示的实时数字高程图与按照下式进行相对干涉高程精度Δhrel的评估：Δhrel=1CN2Σi=1N-1Σj=iN(ΔZi,j-ΔZi,j′)2---<15>]]>其中，ΔZi,j表示第i个和第j个点目标的真实高度差，ΔZi,j′表示实时处理得到的第i个和第j个点目标高度差，得到相对干涉高程精度Δhrel值为3.2556米，精度较高，能满足大多数场合的应用。表明本发明方法可正确的实时估计绝对相位模糊数，得到正确的实时数字高程图。仿真实验2：在半实物仿真环境下对图7进行实时数字高程图反演，结果如图8，其中图8横坐标是方位向像素单元，纵坐标是距离向像素单元，右侧的色度条表示图像显示的大小范围，即场景高度变化范围。对比图7与图8，可以看出：实时数字高程图与基准数字高程图的右侧色度条显示范围是一致的，说明实时数字高程图得到的高程值是正确的。从图像上看，实时数字高程图与基准数字高程图基本一致，说明绝对模糊估计的正确性。将图7与图8进行相对干涉高程精度计算，得相对干涉高程精度Δhrel为4.7439米，精度较高，能满足大多数场合中的应用要求。表明本发明方法可正确的实时估计绝对相位模糊数，得到正确的实时数字高程图。当前第1页1 2 3