一种基于最小二乘克希霍夫偏移的探月雷达数据处理方法与流程

本发明属于月球及深空探测领域，具体地讲，是涉及一种基于最小二乘克希霍夫偏移的探月雷达数据处理方法。

背景技术：

2013年我国嫦娥三号航天器成功着陆，着陆地点位于月海盆地的北部，月球地理坐标为44.12°n，19.50°w。随后，玉兔巡视器从嫦娥三号航天器中着陆月球地面。玉兔巡视器沿着特定路径行驶，探月雷达搭载在巡视器上，发射和接收电磁信号，并将接收到的回波信号传回地球。玉兔巡视器共行驶了大约114m，中途共有17个航点，起点为n101点，终点为n209点。巡视器沿着特定路径行驶，相邻航点之间间隔一小段距离，目的是为了定期刻画出玉兔巡视器的进度。当玉兔巡视器处于固定位置时，探月雷达开始采集信号，并将信号以数据块的形式进行保存，并传回地球。主频为60mhz和500mhz的两组天线同时进行采集，接收信号的时间窗口长度分别为640ns和10240ns。

近年来，国内外科学家围绕探月雷达数据开展了大量的处理和解释工作。其中有代表性的研究工作包括：肖龙等在science期刊上首次对嫦娥三号谈阅历雷达数据开展了基本处理并得到了月球内部的分层结构。张金海，法文哲等先后采用fk偏移，时频分析等技术手段对探月雷达开展精细处理，在获得原有认识的基础上对月球内部分层结构以及岩屑、岩块的散射干扰得到了不错的应用效果。随后，张领等国内学者采用一系列的去噪方法技术对探月雷达双通道数据开展了不同形式的去噪和信号识别研究。

在地质结构比较完整，速度比较简单，没有明显的横向空间变化时可以考虑使用stolt偏移或者fk偏移即可，在地质条件比较复杂的区域，stolt偏移或者fk偏移则不能精确的计算，但嫦娥三号数据是共偏移距数据，基于波动方程的逆时偏移等高精度算法不能应用于探月雷达数据，而常规的stolt偏移或者fk偏移也无法实现数据的处理，因此如何实现在探月雷达数据进行高精准的处理，是本领领域技术人员亟需解决的问题。

技术实现要素：

本发明的目的在于提供一种基于最小二乘克希霍夫偏移的探月雷达数据处理方法，主要解决现有技术中存在的常规的stolt偏移或者fk偏移也无法实现探月数据处理的问题。

为了实现上述目的，本发明采用的技术方案如下：

一种基于最小二乘克希霍夫偏移的探月雷达数据处理方法，包括如下步骤：

(s1)读入探月雷达数据，并对数据进行预处理；

(s2)进行速度谱分析，确定月球内部主要层结构的位置；

(s3)开展克希霍夫的探月雷达数据偏移处理；

(s4)加入最小二乘迭代进行数据偏移算法；

(s5)输出探月雷达数据成像结果。

进一步地，所述步骤(s1)中数据预处理包括以下过程：

(1)数据读取：根据标准存储格式存储9段数据，并逐卷读取道集及位置信息；

(2)数据拼接：将读取到的9段数据拼接成一组完整剖面数据；

(3)道集调整：对得到的剖面数据根据连续同相轴进行道集调整；

(4)道集择取：将调整后的道集进行筛选，选择数据质量好的道集；

(5)时间差调整：发射天线与接收天线启动时间不一致，因此需要对择取道集后的数据进行时间差调整，如果发射天线与接收天线启动时间一致，则跳过该步骤直接执行下一步；

(6)无用数据删除：将调整后的数据信噪比低的无用数据部分删除；

(7)带通滤波：进行带通滤波，去除数据中的低频及高频噪声；

(8)背景去除：由于月球表面反射强烈，存在水平连续的干扰波，因此将滤波去噪后的数据减去平均道，突出异常信息；

(9)自动增益：为突出目标信息，进行自动增益。

进一步地，所述步骤(s2)中进行速度谱分析，确定月球内部主要层结构的位置是通过转换数据由共偏移距道集为中心点道集开展动校正，在水平迭代处理过程中开展速度谱分析，通过时频分析技术，计算出月球内部主要雷达信号反射层的位置。

进一步地，所述步骤(s3)中进行克希霍夫的探月雷达数据偏移处理的主要流程如下：

线性正演模拟算子l～满足

p＝l～m(1)

其中，p表示建模数据的矢量，m表示模型反射率模型矢量，克希霍夫偏移使用公式(1)中的正演模型算子的转置：

mk＝l～^tp(2)

将公式(1)带入公式(2)，可得

mk＝l～^tl～m(3)

其中，l～^tl～是hessian矩阵，对应l～^tl～的积分算子的显式形式表达公式为：

k(x,x′)＝∫dshs(s)asxasx′×∫drhr(r)axrax′rr(τsx+τxr-τsx′-τx′r)(4)

其中，x和x′分别表示与mk和m相关的位置，τsx表示从源s处到x处反射器的波传播时间，τxr表示从反射器到r处接收器的波传播时间，asx和axr表示是振幅项，r表示时间(τsx+τxr-τsx′-τx′r)延迟的相关时间交叉系数，hs(s)和hr(r)分别表示源和接收器采样的函数。

进一步地，所述l～^tl～的积分算子的显式形式中主对角线的表达公式为：

其中，k(x0,x0)表示公式(4)中在x0点的表达形式，asx和axr表示是振幅项，r(0)表示零时间延迟的相关时间交叉系数，hs(s)和hr(r)分别表示源和接收器采样的函数。

具体地，所述步骤(s4)中加入最小二乘迭代进行数据偏移算法，其中，克希霍夫偏移将在偏移中产生偏移畸变，并通过下列公式求解偏移畸变量m：

p(m)＝||l～m-p0||²+ε²||c～m-c～mapr||²(6)

其中，||l～m-p0||²表示数据失配函数，ε²||c～m-c～mapr||²表示正则化项，ε²表示正则化权重，c～表示作用于m的线性算子；

最小化公式(6)的模型表达公式为：

m＝(l～^tl～+ε²c～^tc～)^-1(l～^tp0+ε²c～^tc～mapr)(7)

求解公式(7)：

(l～^tl～+ε²c～^tc～)m＝l～^tp0+ε²c～^tc～mapr(8)

最后采用共轭梯度方案求解关于m的公式(8)，求解的m即是最后的克希霍夫偏移成像输出结果。

与现有技术相比，本发明具有以下有益效果：

(1)本发明克希霍夫(kirchhoff)偏移的主要优势是计算效率高，能够适合于探月雷达的数据观测系统，并且对输入的雷达数据没有特殊的要求，不依赖于速度模型，处理方式方便灵活，适合于目标成像。同时采用基于最小二乘迭代方式的克希霍夫偏移可以提高目标成像分辨率和振幅保真度，在探月雷达数据处理和目标成像定位等方面具有很好的应用价值。

(2)本发明的数据预处理部分采用完整的探月雷达数据处理流程，包括去除重复数据，滤波，增益，去坏道集等过程。保证下一步过程克希霍夫偏移输入探地雷达数据有较高的信噪比和信号强度，同时为提高探月雷达数据偏移成像精度，进行了克希霍夫偏移成像技术，它的主要优势是：第一，能够适应不同的观测系统，对输入地震数据没有特殊要求，处理方式方便灵活，非常适于做目标成像；其次，计算效率较高。对于复杂构造，横向起伏较大的地层分布，常规探地雷达fk偏移方法难以获得准确的成像结果，在起伏较大的构造边缘存在绕射干扰难以消除等问题。而本发明采用克希霍夫偏移可以很好的解决上述问题。

附图说明

图1为本发明的流程图。

图2为本发明探月数据预处理流程图。

图3中a)为理论探月雷达数据介电常数模型；b)为常规fk偏移成像结果。

图4中a)为克希霍夫偏移成像结果；b)为最小二乘克希霍夫偏移成像结果。

图5中a)为月壤结果概念模型；b)为克希霍夫偏移成像结果；c)最小二乘克希霍夫偏移成像结果。

图6中a)为常规fk偏移方法偏移成像结果；b)探月雷达数据克希霍夫偏移成像结果；c)探月雷达数据最小二乘克希霍夫偏移成像结果。

具体实施方式

下面结合附图和实施例对本发明作进一步说明，本发明的实施方式包括但不限于下列实施例。

实施例

如图1至图6所示，一种基于最小二乘克希霍夫偏移的探月雷达数据处理方法，包括如下步骤：

对来源于中科院天文台雷达数据进行读入，其格式采用了航空航天界的标准存储格式(.psd)，数据共分为9段，根据其存储格式，逐卷读取道集及位置信息；将9段数据拼接成一组数；由于卷集不同以及启始采集时间存在误差等问题，个别道集之间存在错位现象，因此根据连续同相轴进行道集调整；众所周知，搭载lpr的“玉兔号”月球车并非匀速行走，在某些航点上，会停下采集其他科学数据(如地形地貌相机等数据)，与此同时，雷达并未停止采集，导致在同一位置上重复采集多道数据，需将同一位置上的重复道集累加求平均，在探月雷达启动之始，由于参数设置问题，第一卷及第二卷道集信噪比极低，亦需要去除从而实现道集择取；由于发射天线与接收天线启动时间不一致，从而进行时间差调整；删除信息极弱且信噪比较低的无用数据部分；然后进行带通滤波，去除低频及高频噪声；由于月球表面反射强烈，存在水平连续的干扰波，需要减去平均道，突出异常信息；同时为了突出下部信息，需要进行自动增益(agc)。

通过转换数据由共偏移距道集为中心点道集，通过开展动校正，水平迭代处理过程开展速度谱分析，通过时频分析技术，计算出月球内部主要雷达信号反射层的位置。

在开展克希霍夫的探月雷达数据偏移处理时，其克希霍夫偏移的主要流程如下：

当线性正演模拟算子l～满足

p＝l～m(1)

其中，p是建模数据的矢量，m是地球反射率模型矢量。观测数据p0由p0＝l～0m0描述，其中m0是真正的地球反射率模型向量，而l～0是实际地球模型的正演模拟算子，除非另有说明，否则假设l～＝l～0；克希霍夫偏移使用等式1中的正演模型算子的转置：

mk＝l～^tp(2)

将公式1带入公式2，可得

mk＝l～^tl～m(3)

矩阵l～^tl～是hessian矩阵，并将mk定义为m的l～^tl～滤波版本，克希霍夫kirchhoff偏移算子将正确地重建实际的地球模型向量，通常情况下l～^tl～不是单位矩阵，并具有主对角线元素是不均匀且不统一，同时主对角线之外的元素是非零的特征，对应l～^tl～的积分算子的显式形式：

k(x,x′)＝∫dshs(s)asxasx′×∫drhr(r)axrax′rr(τsx+τxr-τsx′-τx′r)(4)

其中x和x′分别表示与mk和m相关的位置，τsx是从s处的源到x处的反射器的波传播时间，τxr是从反射器到r处的接收器的波传播时间，asx和axr是振幅项(根据传播方程估算)，w＝w(t)表示子波源随时间的变化，r是时间(τsx+τxr-τsx′-τx′r)延迟的相关的时间交叉系数，而hs(s)和hr(r)分别是源和接收器采样函数。

可以从表达式(4)中观察到：

主对角线的形式为

该表达式表明，主对角线元素(对于给定的子波源)的值取决于扩散损失和由hs(s)和hr(r)决定的波场的不连续采样。kirchhoff偏移的scheme通常会被增加以弥补扩散损失；bleistein(1984)等的分析反演公式由此产生。通过在每个cdp箱中的数字命中计数对k(x0，x0)归一化，可以减轻不连续采样的影响。该补偿方案可以被解释为将对角矩阵(l～^tl～)-1应用为预处理矩阵。

在建模和偏移之后，l～^tl～每列都是对偏移部分mk的m中固定元素的响应，对于一道，偏移响应会沿着由x描述满足τsx+τxr＝τsx′+τx′r的偏移椭圆，对于完整数据，偏移响应是对各道的总响应(以各道相对应的加权系数相加)，并且对于x′处的每个衍射点，偏移产生在x＝x′附近的聚焦图像，在这种情况下，l～^tl～是对角占优矩阵l～^t≈l-1，对于不完整的数据，偏移响应不只围绕x＝x′聚焦，因为缺失的数据导致偏移畸变的不完全消除，因此在l～^tl～中存在明显的非对角元素。在这种情况下，l～^t/≈l-1，需要某种矩阵求逆方法。

在加入最小二乘迭代进行数据偏移算法时，其中，kirchhoff偏移将在偏移中产生偏移畸变，为了最小化畸变，通过最小化以下目标函数来求解m：

p(m)＝||l～m-p0||²+ε²||c～m-c～mapr||²(6)

等式(6)右边的第一项是数据失配函数，第二项是正则化项，c～表示作用于m的线性算子，使用先验信息为每项指定c～，模型的先验信息可以与mapr项合并，其中ε²是正则化权重，表达式(6)当c-滤波模型残差和数据残差分布是高斯分布，则最小化等式(6)的模型是：

m＝(l～^tl～+ε²c～^tc～)^-1(l～^tp0+ε²c～^tc～mapr)(7)

接着解正规方程：

(l～^tl～+ε²c～^tc～)m＝l～^tp0+ε²c～^tc～mapr(8)

最后采用共轭梯度方案求解关于m的公式(8)，求解的m即是最后的克希霍夫偏移成像输出结果。

其中，公式4是为了计算l～^tl～是hessian矩阵，其显示形式用k表示(即公式4)，公式4是其一般形式，在偏移方法中为了简化，即只计算公式5中的x0点的k值，公式6是计算最后偏移输出量m的表达式，7和8是求解的过程，而其中所需的l～^tl～是hessian矩阵是通过公式5求得。

并且附图3是张金海等(pnas，2014)采用fk偏移方法对理论模型偏移成像结果，从中我们可以看出存在很多绕射干扰。附图4是采用本申请的克希霍夫偏移a)和最小二乘克希霍夫偏移b)方法对应附图3中相同模型的偏移成像结果，结合附图3和附图4可以得出本发明方法有更高的成像精度，对于各个目标体周围的干扰能很好的消除。

附图5中a)是采用法文哲(grl，2015)论文中的月壤结构解释模型测试本申请中克希霍夫偏移(图b))和最小二乘克希霍夫偏移(图c))的效果，通过月壤概念模型也验证了本发明方法能对月壤内部结构得到很好的成像结果。

附图6是采用本申请的方法对嫦娥三号探月雷达实测500m数据开展的偏移成像，其中图a)是张金海等(pnas，2014)的结果。图b)和图c)是本申请采用的两种方法结果，可以看出本发明的结果在实测数据应用中也有相当较好的偏移成像结果。

上述实施例仅为本发明的优选实施例，并非对本发明保护范围的限制，但凡采用本发明的设计原理，以及在此基础上进行非创造性劳动而做出的变化，均应属于本发明的保护范围之内。