基于分层搜索和距离空间投影的高精度闪电三维定位方法与流程

本发明属于闪电定位技术领域，更具体地，涉及到一种基于分层搜索和距离空间投影的高精度闪电三维定位方法。

背景技术：

闪电是一类发生在对流天气过程中的强放电现象，是自然界最常见的灾害性天气现象之一。随着微电子设备和电气设备的广泛应用，由雷电造成的直接和间接灾害日益严重，因此，对雷电的探测和预警成为了当前亟待解决的任务之一。

闪电三维定位技术不仅能对地闪进行很好地监测，也能对云闪等早期雷电电磁辐射发展过程进行观测，为雷电发生发展机理的研究和雷电的早期预警提供关键信息，因此，高精度的闪电三维定位技术成为了雷电探测研究领域的热点之一。现有闪电三维定位技术主要包括基于泰勒级数展开的迭代计算定位方法和最优搜索定位方法。基于泰勒级数展开的迭代计算定位方法根据多个观测站测量的一次闪电辐射电场所激发电磁信号的到达时间(toa)/到达时间差(tdoa)，选择其中5个观测站数据联立解算非线性方程组，获得闪电的空间定位信息和时间信息。由于观测站高度接近平面但又不完全相等，该方法在高度上具有较大的定位误差。同时，由于无法确定哪5个站点联立能得到最合适的解，通常采用遍历所有可能组合的方法来判断，因此极大地增加了计算负荷。

美国雷电探测网采用在地球表面逐步搜索雷电位置最优解的算法，全球雷电探测网采用“下山法”搜索雷电位置，这两个探测系统使用的定位算法都属于最优位置搜索算法。网格搜索法也是一种典型的位置搜索算法。网格搜索法采用一定的步长将探测网覆盖区域划分为网格，将网格单元对应的地理空间位置x带入以闪电发生时刻为适应函数的公式中，计算适应函数值，以令适应函数取值最小的网格位置作为最优解。当探测网覆盖范围大，而要求的计算精度高时，为了减小计算负荷，通常采用多级网格划分的策略，如文献[胡志祥，雷电定位算法和误差分析理论研究，2012.5]所述。这种方法存在两点不足之处：第一，选取的适应函数仅利用了时间信息，丢弃了观测站测量的闪电辐射电场所激发的电场值；第二，由于地理空间分辨率随着网格单元到雷电探测网距离的不同而不同，网格单元离探测网越近分辨率越高，离探测网越远分辨率越低，这意味着：同一步长划分下，如果离探测网近的网格单元大小满足定位精度要求，则在远离探测网的区域，由于分辨率的降低，多个网格单元代表同一地理位置，导致地理空间位置的不唯一性。因此，当闪电发生在定位站网外时，直接在地理空间上的网格搜索法不论二维平面误差还是高度垂直误差均达到公里级以上。

技术实现要素：

为了提高闪电三维定位精度，同时又能对闪电快速定位，本发明提供了一种基于分层搜索和距离空间投影的高精度闪电三维定位方法。该方法的基本思路是：采用由粗到细的分层搜索策略，逐层缩小搜索步长和搜索范围，将地理空间网格化后遍历网格单元，以实现对闪电的快速定位；采用距离空间投影算法将地理空间映射为距离空间，将观测站接收的闪电辐射电场值投影到距离空间中，从距离空间投影矩阵中提取闪电位置，利用距离空间分辨率不变的特性，实现对闪电的高精度三维定位。

该基于分层搜索和距离空间投影的高精度闪电三维定位方法，包括如下步骤：

步骤1)：观测站布站：设定一个观测站为中心站，3个最外围的观测站为辅助定位站，包括中心站的其他观测站位于辅助定位站之间连线构成的三角形内；

步骤2)：观测站记录闪电的闪电辐射电场值及其闪电触发时间；

步骤3)：利用初始搜索步长，对闪电观测站网探测范围进行地理空间网格单元的划分；

步骤4)：将地理空间网格单元转化到三维距离空间，获得距离空间网格单元；距离空间网格单元的三个坐标分量是将地理空间网格单元到三个辅助定位站和到中心站的距离差量化后所得值组成的单元；同时在三维距离空间上生成初值为零的距离空间投影矩阵；

步骤5)：根据步骤2)中获得的各观测站记录的闪电触发时间，将相应的闪电辐射电场值逐一投影到三维距离空间中，获得距离空间投影矩阵；

步骤6)：将距离空间投影矩阵中投影值最大者对应的地理空间网格单元的位置，确定为闪电估计位置；缩小搜索步长，对以闪电估计位置为中心的区域重新进行地理空间网格单元的划分、转化和投影操作，获得新的闪电估计位置，重复多次，直到满足结束条件，最终获得闪电定位结果；

步骤7)：输出闪电定位结果及闪电发生时刻。

其中，步骤1)中，观测站布站为：s个闪电观测站分散布置在不同地区，其中s≥5；以观测站1为中心站，辅助定位站选定观测站2、3、4，对观测站所在的地理空间建立三维笛卡尔直角坐标系，坐标原点位于中心站处，观测站的位置坐标已知，分别为rj＝[xj,yj,zj]^t，j＝1,2,...,s为站点编号，[]^t表示矩阵转置，其中，中心站位置坐标为r1＝[0,0,0]^t。

其中，所述步骤2)进一步针对当前记录的一次闪电数据，计算该次闪电到达各观测站与到达中心站的到达时间差，利用该到达时间差校正各观测站的闪电触发时间。

其中，步骤2)中，所述到达时间差采用波形互相关法计算，过程为：观测站记录一次闪电辐射电场所激发的电磁场信号在第j个观测站的触发时间t'＝{t'j}及其对应的辐射电场波形e＝{e(t'j,j)}，其中e(t'j,j)为闪电辐射电场值，j＝1,2,...,s，s为观测站总数，j＝1表示观测站1即中心站，j＝2,3,4表示观测站2,3,4即辅助定位站；根据观测站2,3,...,s接收的闪电辐射电场值e＝{e(t'i,i)}，i＝2,3,...,s和中心站接收的辐射电场值e＝{e(t'1,1)}，采用波形互相关法分别计算观测站2,3,...,s和中心站的波形互相关，获得一次闪电到达观测站2,3,...,s和到达中心站的时间差，记为δt＝{δti1},i＝2,3,...,s；

所述校正为：根据到达时间差校正观测站2,3,...,s接收闪电辐射电场的触发时间，获得校正后的闪电触发时间t＝{ti}，其中ti＝t'i-δti1，及与校正后闪电触发时间对应的闪电辐射电场波形e＝{e(ti,i)}，其中e(ti,i)＝e(t'i-ti1,i)，i＝2,3,...,s。

其中，所述步骤3)和步骤6)中地理空间网格单元的划分具体为：以kδ为搜索步长、rcen为三维地理空间中心点划分地理空间，生成nx×ny×nz个地理空间网格单元，每个地理空间网格单元的坐标为r＝[nx-nx/2,ny-ny/2,nz-nz/2]×kδ+rcen，其中k为地理空间分层搜索放大倍数，初始化k为设定值d，δ为设定的闪电定位精度，nx＝1,2,...,nx，ny＝1,2,...,ny，nz＝1,2,...,nz，nx、ny和nz均为正整数；其中，如果k＝d，rcen为闪电观测站网探测范围的中心位置坐标；如果k<d，rcen为步骤6)获取的对闪电估计位置的坐标。

其中，步骤4)中所述将地理空间网格单元转化到三维距离空间为：

将划分得到的地理空间网格单元的坐标r代入距离差公式ri＝||r-ri||2-||r-r1||2，其中，i＝2,3,4，表示观测站2、3和4即辅助定位站的位置坐标，下角标1表示观测站即中心站的位置坐标，||||2表示euclidean范数；通过所述距离差公式计算得到每个地理空间网格单元到观测站2、3和4和到中心站的距离差r2、r3和r4；以kδ为步长量化距离差r2、r3和r4，生成由观测站2、3、4决定的距离空间网格单元p＝{[n,m,l]}，其中，距离空间网格单元的三个坐标分量分别为k为地理空间分层搜索放大倍数，初始化k为设定值d，δ为设定的闪电定位精度，ceil(·)表示对括号内数字向上求整；

所述在三维距离空间上生成初值为零的距离空间投影矩阵为：根据距离空间网格单元p＝{[n,m,l]}生成距离空间投影矩阵i＝{i[n,m,l]}，i[n,m,l]为投影矩阵中距离空间网格单元[n,m,l]处的投影值，初始化投影矩阵投影值为零，即i＝{i[n,m,l]＝0}。

其中，步骤5)中所述投影为：

根据步骤4)获得的距离空间网格单元p＝{[n,m,l]}和已知的观测站坐标rj＝[xj,yj,zj]^t，j＝1,2,3,4，采用牛顿迭代法解非线性方程组，

获得与每个距离空间单元对应的地理空间坐标x；其中，j＝1表示观测站1即中心站，j＝2,3,4表示观测站2,3,4即辅助定位站

选取一个距离空间单元p1，针对观测站i＝2，将x代入到达时间差方程ti1＝(||x-ri||2-||x-r1||2)/c，计算x到观测站2和到中心站的时间差ti1，再将ti1与步骤2)校正后的第2个观测站的触发时间ti进行比较，如果|ti1-ti|>ε，则不对三维距离空间网格单元[n,m,l]处的投影值i[n,m,l]做任何操作；如果|ti1-ti|≤ε，则令距离空间网格单元[n,m,l]处的投影值i[n,m,l]等于其自身的值加上与时间ti对应的辐射电场值的绝对值|e(ti,i)|；令i遍历3～s，s为观测站总数，完成对所有观测站的累加操作，从而得到距离空间网格单元p1的投影值；其中，ε为设定门限；

针对每个距离空间网格单元计算地理空间坐标x并进行上述累加操作，完成对所有距离空间单元的投影。

其中，步骤6)中闪电估计位置的确定方式为：

找到距离空间投影矩阵i＝{i[n,m,l]}中最大值imax所对应的距离空间网格单元[nmax,mmax,lmax]，根据[nmax,mmax,lmax]找到与其对应的地理空间坐标xmax，即为闪电估计位置。

其中，步骤6)中是否满足结束条件的判定方式为：搜索步长设定为kδ，其中k为地理空间分层搜索放大倍数，k初始值为设定值d，且d为地理空间分层搜索缩小倍数d的幂，δ为设定的闪电定位精度；每次缩小搜索步长时，先判断k是否为1，如果是，则确定满足结束条件；否则，将k除以d。

其中，步骤7)中输出闪电定位结果及闪电发生时刻为：

步骤6)获得的与[nmax,mmax,lmax]对应的地理空间坐标xmax即为估计的闪电所在位置估计的闪电发生时刻为其中t'1为中心站测量到的一次闪电辐射电场所激发的电磁场信号的触发时间，c为电磁波传播速度。

有益效果：

1、与现有技术相比，本发明的创新点在于：针对基于网格搜索的闪电三维定位技术直接在地理空间上进行网格搜索，未考虑地理空间分辨率变化的情况，本发明采用距离空间投影算法将地理空间映射为距离空间，利用距离空间分辨率不变的特性，实现对闪电的高精度三维定位；针对以闪电发生时刻为适应函数，未利用观测站测量的闪电辐射电场激发的垂直电场信息的缺点，本发明将观测站接收的闪电辐射电场值投影到距离空间，在距离空间中积累多个观测站的测量信息，最后从距离空间投影矩阵中根据投影最大值提取闪电位置，充分利用了测量信息；同时，为了减小网格搜索的计算负荷，采用了由粗到细的分层搜索法，实现了对闪电快速、高精度地三维定位。

2、本发明的优点在于采用分层搜索和距离空间投影，分层搜索策略能显著减小计算量，提高定位速度；距离空间投影具有空间分辨率不变的特性，克服了地理空间投影过程中的空变特性带来的定位误差，因此，本发明能是实现对闪电快速地、高精度地三维定位。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明的流程图；

图2为基于分层搜索和距离空间投影的高精度闪电三维定位方法的系统模型；

图3为地理空间和距离空间中空间分辨率变化情况示意图；

图4为本发明实施例的定位误差图。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，一种基于分层搜索和距离空间投影的高精度闪电定位方法，包括如下步骤：

步骤1.观测站布站：设定一个观测站为中心站，3个最外围的观测站为辅助定位站，包括中心站的其他观测站位于辅助定位站之间连线构成的三角形内。

本实施例中，如图2所示，系统由安装在不同地区的s个闪电观测站组成，其中s＝7，以观测站1为中心站，观测站2、3、4为辅助定位站。对观测站所在的地理空间建立三维笛卡尔直角坐标系，坐标原点位于中心站处，观测站位置坐标已知，分别为rj＝[xj,yj,zj]^t，j＝1,2,...,7为站点编号，[]^t表示矩阵转置，其中，中心站位置坐标为r1＝[0,0,0]^t；观测站2、3和4分别位于r2＝[25,0,0]^tkm，r3＝[-25,0,0]^tkm和r4＝[0,43,0]^tkm处，观测站5、6和7位于观测站2、3和4站之间连线构成的三角形内。

其它已知的系统参数包括：地理空间分层搜索放大倍数k，初始化k＝d＝1000；地理空间分层搜索缩小倍数d＝10；闪电定位精度δ＝10m；电磁波传播速度c＝3×10⁸m/s；门限ε＝0.2×10^-6。

步骤2.观测站记录一次闪电的闪电辐射电场值及其闪电触发时间，并进行数据校正。

各观测站记录一次闪电辐射电场所激发的电磁场信号在第j个观测站的触发时间t'＝{t'j}及其对应的辐射电场波形e＝{e(t'j,j)}，其中e(t'j,j)为闪电辐射电场值，j＝1,2,...,7；根据观测站2,3,...,7接收的闪电辐射电场值e＝{e(t'i,i)}，i＝2,3,...,7和中心站接收的辐射电场值e＝{e(t'1,1)}，采用波形互相关法分别计算观测站2,3,...,7和中心站的波形互相关，获得一次闪电到达观测站2,3,...,7和到达中心站的时间差，记为t＝{ti1},i＝2,3,...,7。根据到达时间差校正观测站2,3,...,7接收闪电辐射电场的触发时间，获得校正后的闪电触发时间t＝{ti}，其中ti＝t'i-δti1，及与校正后时间对应的闪电辐射电场波形e＝{e(ti,i)}，其中e(ti,i)＝e(t'i-ti1,i)，i＝2,3,...,7。

步骤3.设定搜索步长，划分地理空间：

以kδ为搜索步长、rcen为三维地理空间中心点划分地理空间。本实施例中，搜索步长的初始值为kδ＝10km，生成100×100×100个地理空间网格单元，每个地理空间网格单元的坐标为r＝[nx-50,ny-50,nz-50]×10×k+rcen，其中nx＝1,2,...,100，ny＝1,2,...,100，nz＝1,2,...,100；rcen＝[xcen,ycen,zcen]表示被划分的地理空间中心处的笛卡尔直角坐标，其中，如果k＝d＝1000，rcen为观测站网探测范围的中心位置坐标，被划分的地理空间为闪电定位站网探测范围；如果k<1000时，rcen为步骤6获取的对闪电位置的估计结果，被划分的地理空间为以rcen为中心的区域。该区域所覆盖的总网格数不变，采用设定值；由于搜索步长变化，因此区域大小随搜索步长变化而变化。

步骤4.将地理空间网格单元转化到三维距离空间，获得距离空间网格单元；并生成距离空间投影矩阵：其中，所述距离空间网格单元是将地理空间网格单元到中心站和到三个辅助定位站的距离差量化后所得值组成的单元。

具体来说：将步骤3获得的r代入距离差公式ri＝||r-ri||2-||r-r1||2，其中，i＝2,3,4，||||2表示euclidean范数，计算每个地理空间网格单元到观测站2、3和4站和到中心站的距离差r2、r3和r4；以kδ为步长量化距离差r2、r3和r4，生成由观测站2、3、4决定的距离空间网格单元p＝{[n,m,l]}，其中，其中，ceil(·)表示对括号内数字向上求整；根据距离空间网格单元生成距离空间投影矩阵i＝{i[n,m,l]}，i[n,m,l]为投影矩阵中距离空间单元[n,m,l]处的投影值，此时只是生成了矩阵框架，投影矩阵的投影值被初始化为零，即i＝{i[n,m,l]＝0}。

步骤5.距离空间投影：根据步骤2中获得的各观测站记录的闪电触发时间，将相应的闪电辐射电场值逐一投影到三维距离空间中，获得距离空间投影矩阵。

本步骤的操作是一个循环累加的过程。

步骤51：从步骤4获得的距离空间p＝{[n,m,l]}中选取一个距离空间网格单元p1[n,m,l]，利用已知的观测站坐标rj＝[xj,yj,zj]^t，j＝1,2,3,4，采用牛顿迭代法解非线性方程组，

获得与距离空间网格单元[n,m,l]对应的地理空间坐标x；

步骤52：将x代入到达时间差方程ti1＝(||x-ri||2-||x-r1||2)/c，i＝2，计算x到观测站2和到中心站的时间差ti1(i＝2)，再将ti1与步骤2校正后的观测站2的触发时间ti(i＝2)进行比较，如果|ti1-ti|>ε，则不对距离空间单元[n,m,l]处的投影值i[n,m,l]做任何操作；如果|ti1-ti|≤ε，则将观测站2中与这个触发时间对应的绝对电场值进行投影，即令距离空间单元[n,m,l]处的投影值等于其自身的值加上与时间ti对应的辐射电场值的绝对值，即：

i[n,m,l]＝i[n,m,l]+e(ti,i)

其中，||表示求模运算。

步骤53：对同一个x，再代入时间差方程，计算它到中心站和到观测站3的时间差，完成对观测站3的累加；重复上述步骤52，直到与所有观测站都比较完，完成对所有观测站的累加操作，从而得到距离空间网格单元p1的投影值。然后进入步骤54。

步骤54：从距离空间p＝{[n,m,l]}中选取另一个距离空间网格单元[n,m,l]，执行步骤51～53，直到完成对所有距离空间网格单元的投影，获得投影矩阵i＝{i[n,m,l]}。

步骤6.从距离空间投影矩阵中提取闪电估计位置：

找到距离空间投影矩阵i＝{i[n,m,l]}中最大值imax所对应的距离空间网格单元坐标[nmax,mmax,lmax]，根据[nmax,mmax,lmax]找到与其对应的地理空间坐标xmax，估计的闪电地理空间坐标即为xmax。如果k＝1，则判定满足结束条件，转到步骤7；否则，令k＝k/d，rcen＝xmax，转到步骤3重复进行地理空间网格单元的划分、转化和投影操作。

步骤7.输出定位结果：步骤6获得的与[nmax,mmax,lmax]对应的地理空间坐标xmax即为估计的闪电所在位置估计的闪电发生时刻为其中t'1为中心站测量到的一次闪电辐射电场所激发的电磁场信号的触发时间。

本发明实施例主要采用仿真实验的方法进行验证，所有步骤、结论都在matlabr2012a上验证正确。

为了说明空间分辨率随网格单元到观测站网距离的变化情况，在相同的观测站网布站情况下，令闪电发生位置分别位于观测站网内和站网外，基于本发明实施例得到了图3，其中，图3(a)(d)(g)(j)(m)为闪电与观测站网几何位置图，x、y和z轴为三维地理空间位置坐标轴，圆圈表示观测站，星号表示闪电；图3(b)(e)(h)(k)(n)为地理空间中的一个分辨单元，x、y和z轴分别为一个分辨单元的三维坐标；图3(c)(f)(i)(l)(q)为距离空间中一个分辨单元，x、y和z轴分别为一个分辨单元的三维坐标。图3(a)闪电位于[0,0,10]km，图3(d)闪电位于[20,35,10]km，图3(b)和(c)分别是与(a)对应的地理空间和距离空间，图3(e)和(f)分别是与(d)对应的地理空间和距离空间，可以看出，当闪电位于观测站网内时，不论地理空间还是距离空间，空间分辨率均不变，即一个分辨单元等于一个网格单元；图3(g)闪电位于[50,80,10]km，图3(j)闪电位于[200,0,10]km，图3(m)闪电位于[200,200,10]km，图3(h)和(i)分别是与(g)对应的地理空间和距离空间，图3(k)和(l)分别是与(j)对应的地理空间和距离空间，图3(n)和(q)分别是与(m)对应的地理空间和距离空间，从图中可以看出，随着闪电离观测站网距离的增加，地理空间中的分辨率逐渐减小，一个分辨单元中包括了多个网格单元；而距离空间中，闪电处的分辨率几乎没有变化，如图3(i)、(l)和(q)所示。

图4为本发明实施例的定位误差等高线图，图中，三角形表示观测站点，连线为误差的等高线，等高线值在图中已标示，图4(a)为x-y方向误差在x-y面上的等高线投影图，图4(b)为z方向误差在x-y面上的等高线投影图。从图中可以看出，在发明实施例布站情况下，以中心站为坐标原点，x方向从-200km到200km，y方向从-200km到200km，z方向从5km到15km内，不论是在观测站网内还是在观测站网外，x-y方向平面上的水平误差最大为0.55km，z方向上的高度误差最大为1.1km。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。