本发明涉及地层岩性、地质构造地球物理勘探领域，特别涉及一种利用自然闪电为刺激源的大尺度地质构造探测方法。
背景技术：
地球物理勘探通过研究和观测各种地球物理场的变化来探测地层岩性、地质构造等地质条件。由于组成地壳的不同岩层介质往往在密度、弹性、导电性、磁性、放射性以及导热性等方面存在差异，这些差异将引起相应的地球物理场的局部变化。通过量测这些物理场的分布和变化特征，结合已知地质资料进行分析研究，就可以达到推断地质性状的目的。地球物理勘探是以岩石、矿石(或地层)与围岩的物理性质差密度、磁化性质、导电性、放射性差异为基础，利用探测到的人工或天然数据信息进行反演，并结合信息类型与地质参数之间的关系获取地层结构。但信号的发射传输随尺度的增加迅速衰减，且反演过程中信息的解译方法仍具有局限性，因此利用传统的地球物理勘探方法进行大尺度条件下地层结构的精细刻画时还存在局限与不足。
水力层析扫描是近年来兴起的一种有效的刻画地层三维非均质构造的探测手段刻画含水介质水力参数三维空间分布的创新方法。借助医学及地球物理领域的层析成像原理，水力层析扫描技术利用自然界或者人为对渗流场施加刺激，同时在不同地点观测该刺激带来的响应。将多次刺激过程所搜集的有价值信息通过数据融合方法吸收，逐步提高含水介质的水力参数空间分辨率。水力层析扫描的本质是通过观测值反推出非均质参数分布规律的数据融合与反分析算法。它基于贝叶斯地统计学原理解译监测数据，以连续线性估计算法吸收含水层水力参数反演中的非线性过程，以观测数据和待估参数的协相关关系为主要依据更新参数后验均值和协方差，可以方便吸收各种类型的数据。连续线性估计算法以迭代的方法对数据的信息进行线性化，在迭代过程中，在时间上采用同时吸收模式，可以有效解决非线性问题，避免传统数据同化算法等可能出现的参数-观测不一致问题。从虚拟理想算例，到室内试验验证，再到野外试验验证，都印证了该方法可以有效地刻画地质结构的三维空间分布。
然而，水力层析扫描技术对于小尺度含水层非均质特性有良好的刻画，但是在实际的野外探测中研究的区域往往都是大尺度，再以水力作为刺激源会存在刺激强度低、响应慢、要求监测井数量多等问题。
因此，亟需开发一种刺激强度高、响应快的刺激源的大尺度地质构造三维层析扫描方法。
技术实现要素：
本发明的目的是提供一种以闪电为刺激源的大尺度地质构造三维层析扫描方法，以解决现有技术中存在的问题。
为实现本发明目的而采用的技术方案是这样的，一种以闪电为刺激源的大尺度地质构造三维层析扫描方法，包括以下步骤：
1)在待探测区域内布设网状电势探测井，并在每个探测井内间隔布置电势测量电极。其中，所述电势测量电极与电势测量装置相连。
2)测量待探测区域上空的气压值。
3)发射引雷火箭。所述引雷火箭到达雷雨云后激发对地闪电，并将对地闪电引至待探测区域。其中，所述引雷火箭尾部线轴上拖带的导线与待探测区域内预设的地电极相连接。
4)测量对地闪电的电压和能量参数。测量待探测区域在对地闪电作用下不同平面位置及不同深度处的电势响应。
5)重复步骤3)和4)，在探测区域的其它位置依次发射引雷火箭，并测量相应数据。
6)建立探测区域三维地质正、反分析模型，利用连续线性估计算法融合闪电刺激作用下获得的电势响应数据，进行探测区域地层三维结构反演。
7)利用电阻率在地层中的三维分布规律，对地层的构造特性。
进一步，步骤6)具体包括以下步骤：
6.1)根据前期勘探数据，给出地层电阻率的均值、方差和相关尺度参数。
6.2)将电阻率的均值赋予分析区域内每一个网格内。结合边界条件，利用探测区域三维电场正分析模型，计算每一次闪电刺激下每个探测位置在估计参数条件下的电压电流值。
6.3)利用连续线性估计算法对探测区域内待估点的电阻率值进行估计。其中，未知电阻率的迭代计算公式如式(1)所示。权重系数矩阵ω的计算公式如式(2)所示。
式中，r为迭代次数。t为转置符号。uc为待估的探测区域电阻率参数向量。
[εdd+λdiag(εdd)]ω＝εdu(2)
式中，εdd为观测数据之间的协方差矩阵。εdu为观测数据与参数之间的协方差矩阵。λ为levenberg-marquardt算法动态乘子。diag()函数用于构造一个对角矩阵。
6.4)重复式(1)的迭代过程，直至得到的参数估计值条件下得到的计算值与探测值之差小于误差允许值或达到一定的迭代步数。
进一步，步骤6.3)中，协方差矩阵由敏感度得到：
式中，jdu是电压电流观测数据对探测区域电阻率变化的敏感度矩阵。参数协方差矩阵εuu在r＝0时由先验地质信息给出，随后每次迭代按照式(4)逐步更新。
进一步，步骤7)之后，还具有利用未知地层结构参数与电阻率之间的关系，得到地层地质参数的三维分布的相关步骤。
本发明的技术效果是毋庸置疑的：利用电极孔测得的电阻率变化数据可以进行高效的地下含水层电阻率层析扫描，得到的结果可以吸收多种的探测数据源信息，并实现利用观测数据对未知电阻率参数的非线性估计，对提升大尺度含水层地质参数模拟的精度效果明显。
附图说明
图1为工作流程图；
图2求解域及网格划分；
图3电阻率层析扫描迭代过程曲线；
图4不同迭代次数下不确定性的变化情况；
图5不同迭代次数下渗透系数反演刻画情况；
图6渗透系数真实值与反演值对比。
具体实施方式
下面结合实施例对本发明作进一步说明，但不应该理解为本发明上述主题范围仅限于下述实施例。在不脱离本发明上述技术思想的情况下，根据本领域普通技术知识和惯用手段，做出各种替换和变更，均应包括在本发明的保护范围内。
实施例1：
参见图1，本实施例公开一种以闪电为刺激源的大尺度地质构造三维层析扫描方法，包括以下步骤：
1)在待探测区域内，按一定的间距布设网状电势探测井，并在每个探测井内间隔布置电势测量电极。其中，所述电势测量电极与电势测量装置相连。
2)利用探测气球测量待探测区域上空的气压值。利用雷达站用于追踪雨层云和积雨云的位置，进而预测闪电发生位置以及闪电能量。
3)在待探测区域内，发射引雷火箭。所述引雷火箭到达雷雨云后激发对地闪电，并将对地闪电引至待探测区域。其中，所述引雷火箭尾部线轴上拖带的导线与待探测区域内预设的地电极相连接。
4)测量对地闪电的电压和能量参数。测量待探测区域在对地闪电作用下不同平面位置及不同深度处的电势响应。包括电压响应与电流响应等。
5)重复步骤3)和4)，在探测区域的其它位置依次发射引雷火箭，并测量相应数据。将对地闪电引入区域内的不同预设地电极上。
6)建立探测区域三维地质正、反分析模型，利用连续线性估计算法融合闪电刺激作用下获得的电势响应数据，进行探测区域地层三维结构反演。
6.1)根据前期勘探数据，给出地层电阻率的均值、方差和相关尺度参数。
6.2)将电阻率的均值赋予分析区域内每一个网格内，结合边界条件，利用探测区域三维电场正分析模型，计算每一次闪电刺激下每个探测位置在估计参数条件下的电压及电流值。
6.3)利用连续线性估计算法对探测区域内待估点的电阻率值进行估计，未知电阻率的迭代计算公式为：
式中，uc为待估的探测区域电阻率参数向量。t为转置符号。上式左端为参数向量uc在第r+1次的条件估计值，r＝0的参数来源于步骤6.2)正分析的结果以及先验资料。u的估计值在每个迭代步骤逐步由每次闪电作用下电压电流的观测值d*以及利用三维电场正分析模型的模拟值d(r)的差异进行调整，权重系数矩阵ω计算公式为：
[εdd+λdiag(εdd)]ω＝εdu(2)
其中，εdd是观测数据之间的协方差矩阵，而εdu是观测数据与参数之间的协方差矩阵。λ是levenberg-marquardt算法动态乘子，而diag()运算符代表取对角阵。协方差由敏感度得到：
式中，jdu是电压电流观测数据对探测区域电阻率变化的敏感度矩阵，参数协方差矩阵εuu在r＝0时由先验地质信息给出，随后每次迭代按照下式逐步更新：
6.4)重复式(1)的迭代过程，直至得到的参数估计值条件下得到的计算值与探测值之差小到一定的误差允许值，或达到一定的迭代步数。此时的参数估计值即为层析扫描的结果。
7)利用电阻率在地层中的三维分布规律，对地层的构造特性。利用未知地层结构参数与电阻率之间的关系，得到地层地质参数的三维分布。
实施例2：
本实施例选取一块50km×25km的含水层区域作为研究区域并对其进行数值建模，利用层析扫描方法对含水层的渗透系数分布进行刻画。假设研究区域南北两侧为水头边界，在边界围成的求解域内进行单元网格划分，网格为正方形，尺度为1km×1km，共得到有限计算网格1250个。假设在求解域中共布置了40个电极孔，划分网格后的求解域如图2所示。
为了能够使用电阻率层析扫描这一技术，需要首先猜测地层渗透性的大体范围，猜测值越接近真实值，则层析扫描所需要的时间越短。假设根据研究区域若干电极孔的电阻率试验结果已经得到了地层渗透系数的平均值为10m/day，代入数值模拟软件的反分析模块分析研究区域的渗透系数。扫描过程的迭代曲线如图3所示。可以看出，在迭代计算80次后计算收敛，即得到了符合误差条件的反分析结果。
在电阻率层析扫描过程中，计算过程中得到的不确定性数据反应了反分析误差的大小。图4为80次迭代计算过程中计算区域内不确定性分布规律的变化。4a为第1次迭代，4b为第10次迭代，4c为第20次迭代，4d为第40次迭代，4e为第60次迭代，4f为第80次迭代。由图4可知随着迭代次数的增加，研究区域的不确定性呈显著下降趋势，而在分析结束时，即第80次迭代计算后，几乎整个研究区域的不确定性都很小，即说明整体计算结果的收敛性较好。
不确定性的不断降低说明在层析过程中渗透系数值逐渐从猜测值向真实值的进化。如图5所示，在80次迭代计算中，计算区域内的渗透系数值逐渐由初始统一的猜测值(即10m/day)向明显的非均质去变化，到第80次迭代结束时，渗透系数的分布已经可以较为清楚地展示出渗透系数地空间变异特性。5a为第1次迭代，5b为第20次迭代，5c为第40次迭代，5d为第80次迭代。
图6中，6a为真实渗透系数分布，6b为反演得到的渗透系数分布。对比渗透系数真实分布与利用层析扫描得到的渗透系数分布，可以看出，反分析得到的渗透系数在中间部分具有高值区，与实际渗透系数分布具有很好的一致性。
值得说明的是，电阻率层析扫描技术是利用有限的监测数据反演得到高解析度含水层特性的先进手段。通过对比分析，可以看出利用电极孔测得的电阻率变化数据可以进行高效的地下含水层电阻率层析扫描，得到的结果对提升大尺度含水层地质参数模拟的精度效果明显。