本发明属于航空重力勘探，涉及一种航空重力测网交叉点搜索方法。
背景技术：
航空重力勘探是一种广泛采用的航空物探方法。实际工作中，通常采用布设垂直交叉的测线和切割线的方法，形成航空重力测网。开展航空重力测量精度评定和测线平差之前，首先需要快速、准确地搜索出航空重力测网的交叉点。
《航空重力测量测线网平差的理论与方法》(李海.2002.[硕士论文].郑州:中国人民解放军信息工程大学.p18-21.)、《航空重力测量理论方法及应用研究》(孙中苗.2004.[博士论文].郑州:中国人民解放军信息工程大学.p114.)、《航空重力测量网络平差方法研究》(蔡劭琨.2009.[硕士论文].长沙:国防科技大学.p33-36.)等文献公开了一种广泛采用的方法，首先根据航空重力测网计算出近似交叉点，然后在设置范围内遍历搜索出实际交叉点，最后计算获得交叉点信息。该方法只能在测线与切割线在南北、东西向垂直分布的条件下使用，不能用于其他沿任意方向布设测线的情况；同时由于该方法采用遍历搜索方式，搜索速度较慢，只能满足小型的试验区计算要求，对于几万、十几万公里的航空重力测区计算，工作效率十分低下。
技术实现要素：
本发明的目的在于克服现有技术的缺陷，提供一种计算过程简洁、计算效率高、计算结果准确度高的航空重力测网交叉点非遍历搜索方法。
为实现上述目的，本发明采用了如下技术方案：
一种航空重力测网交叉点非遍历搜索方法，包括以下步骤：加载测网数据；逐一计算交叉点，直至完成所有交叉点的计算；输出交叉点的信息。
进一步，所述加载测网数据的方法包括以下步骤：读取测网数据；提取测线和切割线。
进一步，所述计算交叉点的方法包括如下步骤：初步计算交叉点位置；寻找测线、切割线上距离交叉点位置最近的4个目标测点；通过获得的4个目标测点，采用跳跃逼近法第二次计算交叉点位置；更新测线、切割线上距离交叉点位置最近的4个目标测点；判断目标测点之间是否存在交叉点；若否，则微调4个目标测点，并再次判断，直至4个目标测点之间存在交叉点；计算并缓存交叉点位置信息和重力场值信息。
进一步，所述初步计算交叉点位置的方法为：对测线和切割线两端进行外扩；通过测线两端点的投影坐标构建第一直线；通过切割线两端点的投影坐标构建第二直线；计算第一直线和第二直线的相交点，即为交叉点位置。
进一步，所述寻找距离交叉点位置最近的4个目标测点的方法为：通过测线两端点的投影坐标，计算测线总长度及测点之间的平均点距；通过测线一个端点至交叉点位置之间的距离，选定测线上距离交叉点位置最近的2个测点；通过切割线两端点的投影坐标，计算切割线总长度及测点之间的平均点距；通过切割线一个端点至交叉点位置之间的距离，选定切割线上距离交叉点位置最近的2个测点；所述测线上距离交叉点位置最近的2个测点和切割线上距离交叉点位置最近的2个测点，即为4个目标测点。
进一步，所述采用跳跃逼近法第二次计算交叉点位置的方法为：通过选定的测线上距离交叉点位置最近的2个目标测点的投影坐标，构建第三直线；通过选定的切割线上距离交叉点位置最近的2个目标测点的投影坐标，构建第四直线；计算第三直线和第四直线的相交点，即为采用跳跃逼近法第二次计算得到的交叉点位置。
进一步，所述判断目标测点之间是否存在交叉点的方法为：判断选定的测线上距离交叉点位置最近的2个目标测点，是否真正位于切割线两侧；若否，则判定4个目标测点之间不存在交叉点；若是，则进一步判断选定的切割线上距离交叉点位置最近的2个目标测点，是否真正位于测线两侧；若否，则判定4个目标测点之间不存在交叉点；若是，则判定4个目标测点之间存在交叉点。
进一步，所述微调4个目标测点的方法为：根据测线或切割线上的目标测点与交叉点位置之间的距离及方向，采取+1或-1个测点位置的方式，滑动微调测线或切割线上2个相邻的目标测点。
进一步，所述判断选定的测线上距离交叉点位置最近的2个目标测点是否真正位于切割线两侧的方法为：记选定的测线上距离交叉点位置最近的2个测点分别为l1、l2，记选定的切割线上距离交叉点位置最近的2个测点分别为t1、t2；计算△t1t2l1和△t1t2l2的向量面积的正、负符号；若△t1t2l1和△t1t2l2的向量面积符号相反，则判定选定的测线上距离交叉点位置最近的2个测点真正位于切割线两侧。
进一步，所述判断选定的切割线上距离交叉点位置最近的2个目标测点是否真正位于测线两侧的方法为：记选定的测线上距离交叉点位置最近的2个测点分别为l1、l2，记选定的切割线上距离交叉点位置最近的2个测点分别为t1、t2；计算△l1l2t1和△l1l2t2的向量面积的正、负符号；若△l1l2t1和△l1l2t2的向量面积符号相反，则判定选定的切割线上距离交叉点位置最近的2个测点真正位于测线两侧。
本发明一种航空重力测网交叉点非遍历搜索方法，首先计算获得近似交叉点，并在测线和切割线上寻找离交叉点最近的测点；其次根据寻找到的目标测点，采用跳跃逼近法再次计算近似交叉点，第二次逼近交叉点；然后判断目标测点是否真正相交，如不相交又则对其进行微调滑动至完全相交。本发明避开了传统的遍历搜索方式，使计算效率大幅提升。例如某测网数据包含约16万测线公里、546万个测点，数据读取整理时间约9.6s，29517个交叉点搜索及信息提取耗时不足0.09s，有效缩短了航空重力测网交叉点搜索的耗时。
附图说明
图1是本发明一种航空重力测网交叉点非遍历搜索方法的流程示意图；
图2是本发明一种航空重力测网交叉点非遍历搜索方法中计算交叉点的方法的流程示意图；
图3是实施例3中两向量叉乘示意图；
图4是实施例3中交叉点二次逼近后的两个双重循环搜索流程图。
具体实施方式
以下结合附图1至4，进一步说明本发明一种航空重力测网交叉点非遍历搜索方法的具体实施方式。本发明一种航空重力测网交叉点非遍历搜索方法不限于以下实施例的描述。
实施例1：
本实施例给出一种航空重力测网交叉点非遍历搜索方法的整体流程。
如图1所示，一种航空重力测网交叉点非遍历搜索方法，包括以下步骤：加载测网数据；逐一计算交叉点，直至完成所有交叉点的计算；输出交叉点的信息。
其中，所述加载测网数据的方法包括以下步骤：读取测网数据；提取测线和切割线。
其中，所述计算交叉点的方法如图2所示，包括如下步骤：初步计算交叉点位置；寻找测线、切割线上距离交叉点位置最近的4个目标测点；通过获得的4个目标测点，采用跳跃逼近法第二次计算交叉点位置；更新测线、切割线上距离交叉点位置最近的4个目标测点；判断目标测点之间是否存在交叉点；若否，则微调4个目标测点，并再次判断，直至4个目标测点之间存在交叉点；计算并缓存交叉点位置信息和重力场值等相关信息。
其中，所述初步计算交叉点位置的方法为：对测线和切割线两端进行外扩；通过测线两端点的投影坐标构建第一直线；通过切割线两端点的投影坐标构建第二直线；计算第一直线和第二直线的相交点，即为交叉点位置。对测线和切割线两端进行外扩的目的，是避免由于飞行或测量误差导致的测线和切割线上靠近端点位置的交叉点丢失。优选的，所述测线和切割线的两端各外扩3个测点距离。
其中，所述寻找距离交叉点位置最近的4个目标测点的方法为：通过测线两端点的投影坐标，计算测线总长度及测点之间的平均点距；通过测线一个端点至交叉点位置之间的距离，选定测线上距离交叉点位置最近的2个测点；通过切割线两端点的投影坐标，计算切割线总长度及测点之间的平均点距；通过切割线一个端点至交叉点位置之间的距离，选定切割线上距离交叉点位置最近的2个测点；所述测线上距离交叉点位置最近的2个测点和切割线上距离交叉点位置最近的2个测点，即为4个目标测点。由于飞行速度不是完全精准，可能会导致测点的点距或疏密度不同，因此所述4个目标测点有可能不会真正的相交，需要在后续步骤中进行验证和修正。
其中，所述采用跳跃逼近法第二次计算交叉点位置的方法为：通过选定的测线上距离交叉点位置最近的2个目标测点的投影坐标，构建第三直线；通过选定的切割线上距离交叉点位置最近的2个目标测点的投影坐标，构建第四直线；计算第三直线和第四直线的相交点，即为采用跳跃逼近法第二次计算得到的交叉点位置。采用跳跃逼近法替代传统的遍历搜索方法，可以更加迅速地找到交叉点，减小计算工作量。在实际应用中，通过进行一次跳跃逼近，便可以直接找到真正交叉点，或者定位到距离交叉点1至2个点距附近。
其中，所述判断目标测点之间是否存在交叉点的方法为：判断选定的测线上距离交叉点位置最近的2个目标测点，是否真正位于切割线两侧；若否，则判定4个目标测点之间不存在交叉点；若是，则进一步判断选定的切割线上距离交叉点位置最近的2个目标测点，是否真正位于测线两侧；若否，则判定4个目标测点之间不存在交叉点；若是，则判定4个目标测点之间存在交叉点。
其中，所述微调4个目标测点的方法为：根据测线或切割线上的目标测点与交叉点位置之间的距离及方向，采取+1或-1个测点位置的方式，滑动微调测线或切割线上2个相邻的目标测点。
其中，所述判断选定的测线上距离交叉点位置最近的2个目标测点是否真正位于切割线两侧的方法为：记选定的测线上距离交叉点位置最近的2个测点分别为l1、l2，记选定的切割线上距离交叉点位置最近的2个测点分别为t1、t2；计算△t1t2l1和△t1t2l2的向量面积的正、负符号；若△t1t2l1和△t1t2l2的向量面积符号相反，则判定选定的测线上距离交叉点位置最近的2个测点真正位于切割线两侧。
其中，所述判断选定的切割线上距离交叉点位置最近的2个目标测点是否真正位于测线两侧的方法为：记选定的测线上距离交叉点位置最近的2个测点分别为l1、l2，记选定的切割线上距离交叉点位置最近的2个测点分别为t1、t2；计算△l1l2t1和△l1l2t2的向量面积的正、负符号；若△l1l2t1和△l1l2t2的向量面积符号相反，则判定选定的切割线上距离交叉点位置最近的2个测点真正位于测线两侧。
实施例2：
本实施例给出一种对实施例1所述航空重力测网交叉点非遍历搜索方法的改进。
根据实施例1中所述计算交叉点的方法，在采用跳跃逼近法第二次计算交叉点位置、更新4个目标测点的步骤之后，增加以下步骤：
判断目标测点之间是否存在交叉点；
若不存在交叉点，则进一步判断目标测点与交叉点之间的距离；
若距离超过2个点距，则采用跳跃逼近法第三次计算交叉点位置、再次更新4个目标测点；本步骤可以多次迭代重复进行；
若距离不超过2个点距，则采用滑动微调的方式，微调4个目标测点。
本实施例采用了三次或多次迭代使用跳跃逼近法寻找交叉点位置，可以进一步提高计算效率，特别是针对飞行速度变化大、点距分布不均匀的测网数据，具有更快的计算速度。
实施例3：
本实施例给出实施例1所述航空重力测网交叉点非遍历搜索方法的计算机程序实现。本实施例使用c#编程语言，visualstudio2013编译和调试。
本实施例设计了两个交叉点计算函数，满足交叉点搜索时的不同需求。
第一个为iscrosspoint()函数。该函数可以返回测线线段和切割线线段交叉点位置分布情况中的两个判断和最终的交叉点坐标，返回3种结果；第二种为getcross()函数，只返回两条线段的交叉点坐标，交叉点可以分布在线段内外。
利用两向量叉乘等于两向量构成的平行四边形(以两向量为邻边)的面积，由于向量是有方向的，所以面积也是有方向的，通常我们以逆时针为正，顺时针为负数。三角形l1l2t1、l1l2t2面积公式分别为：
sl1l2t1＝((l1.x–t1.x)*(l2.y–t1.y)-(l1.y–t1.y)*(l2.x–t1.x))/2(1)sl1l2t2＝((l1.x–t2.x)*(l2.y–t2.y)-(l1.y–t2.y)*(l2.x–t2.x))/2(2)
如图3所示，如果“测线线段l1l2和切割线测点t1构成的三角形面积”与“测线线段l1l2和切割线测点t2构成的三角形面积”构成的三角形面积的正负符号相异，那么切割线测点t1、t2位于测线线段l1l2的两侧；反之，切割线测点t1、t2位于测线线段l1l2的同侧。
iscrosspoint(t1,t2,l1,l2)函数设计了3种返回结果，通过输入切割线和测线线段的4个端点，先以切割线线段t1t2为基准，判断测线l1、l2端点是否相交，如果两个端点在同一侧，代表与线段l1l2不相交则返回常数c_tie，函数退出；如果l1、l2端点分布在两侧，代表线段l1l2被相交函数继续进行，此时以l1l2为基准，判断t1、t2端点分布情况，同理线段t1t2不被相交时返回常数c_lie，函数退出；如果线段t1t2被相交时，则返回交叉点位置系数，见公式(5)，函数退出。
st1t2l1＝((t1.x–l1.x)*(t2.y–l1.y)-(t1.y–l1.y)*(t2.x–l1.x))/2(3)st1t2l2＝((t1.x–l2.x)*(t2.y–l2.y)-(t1.y–l2.y)*(t2.x–l2.x))/2(4)iscp＝sl1l2t1/(st1t2l2-st1t2l1)(5)
第二个为getcross()函数。getcross(t1,t2,l1,l2)交叉点函数，是一个常规的交叉点计算返回函数，输入两条线段的4个端点后，不区分两线段是否相交情况，返回的交叉点坐标可以在线段内外。
本实施例采用的方法具体包括以下步骤：
1.读取测网数据到测线类中，首先需要建立一个具有盛放测线相关信息的测线类数据类型，具体信息有：飞行架次、飞行日期、测线类型、测线号、数据信息组的名称、测点数据信息组集合(含坐标x、y和场值数据)；又对数据信息建立了测点数据结构体(含坐标位置和场值)和坐标位置结构体(含坐标x、y)；为了节省资源又便于测线集合的收集和测点集合的收集，这里使用了测线类泛型，测点结构体泛型；每条测线和切割线都存储在测线类中，数据体中包含上述列出的概况信息和数据。
2.继承测线类扩展出运算测线类，添加信息有：测线外扩首尾两测点坐标、测线平均点距、测线方向趋势(用于区分偏向x、y方向)、交叉点信息组集合；建立交叉点结构体，内有：坐标、切割线线号+架次号、切割线上点位、切割线上场值、切割线上梯度、测线线号+架次号、测线上点位、测线上场值、测线上梯度、交叉点不符场值；丰富运算测线类中能够暂时获取到的数据信息，再从原测线类中提取出测线和切割线，使之成为测线和切割线的运算测线类，同样使用运算测线类泛型。
3.使用iscrosspoint()函数计算外扩端点后的切割线与测线之间的交叉点，不相交的情况iscrosspoint()函数可以很便捷进行排除；获取交叉点后，测线和切割线通过各自的平均点距寻找离交叉点最近的各2个测点；从刚获取的测线和切割各2个测点，使用getcross()函数再次计算交叉点，同样取出离交叉点最近的各2个测点，第二次逼近已经很逼近交叉点，也可能已经符合交叉情况。
4.通过iscrosspoint()函数对第二次获得的测线和切割线各2个测点进行交叉点验证，如果返回交叉点系数则搜索结束；如果返回常数c_tile，表明测线上两个测点没有相交于切割线，进入双重循环while(c_tile){while(c_line){}的第一个循环中进行搜索，使用getcross()函数计算交叉点位置进行判断，超过2个点距采取第三次逼近交叉点，然后使用微调加1或减1进行滑动逼近交叉点，再使用iscrosspoint()函数进行判断，同样返回交叉点系数则搜索结束，返回c_tile则继续上面的循环直至搜索结束，返回常数c_line进入双重循环的第二个while(c_line){}中，这次是切割线没有相交于测线，过程与上面类似返回三种结果后依旧在双重循环中进行判断直至搜索结束，此时第一个双重循环结束；如果返回常数c_line，直接进入第二个双重循环while(c_line){while(c_tile){}}中进行搜索直至结束。详细流程如图4所示。
5.确定四个相交端点后，提取交叉点坐标(见公式(6)、(7))，测线和切割线线号，飞行架次信息，测线和切割线上的测点位置插值，测线和切割线上的场值插值，测线和切割线上的场值梯度，测线和切割线上的交叉点场值不符值，将这些交叉点信息和统计报告输出。
p.x＝t1.x+iscp*(t2.x-t1.x)(6)
p.y＝t1.y+iscp*(t2.y-t1.y)(7)
实施例4：
本实施例给出实施例3所述方法的一种实际应用结果。
某测区由tazb01、tazb02两个子区组成，各约8万测线公里，合并区16万测线公里(546万多个测点数据)。本方法读取合并区数据及整理所需时间约9.6s，29517个交叉点搜索及信息提取耗时不足0.09s，整个搜索过程约10s。下表列出了本方法与现有成熟软件之间交叉点搜索的耗时情况，本方法大幅度缩短了航空重力测网交叉点搜索的运行时间。本实施例运行环境如下：联想m8300t，intel酷睿i72600，内存4g，2011年出厂。
表1某测区16万测线公里的测网交叉点搜索耗时统计
以上内容是结合具体的优选实施方式对本实施例所作的进一步详细说明，不能认定本实施例的具体实施只局限于这些说明。对于本实施例所属技术领域的普通技术人员来说，在不脱离本实施例构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本实施例的保护范围。