屈進紅,姜作喜,周錫華,羅 鋒,李 芳

1. 自然資源部航空地球物理與遙感地質重點實驗室，北京 100083； 2. 中國自然資源航空物探遙感中心，北京 100083

航空重力測量以飛機為載體，融合了重力傳感器、慣性系統(tǒng)(INS)和差分全球定位系統(tǒng)(DGNSS)獲取近地空中重力加速度的一種動態(tài)重力測量方法[1-6]。在中、高頻重力場信息方面，彌補了地面重力、衛(wèi)星測高在獲取重力場信息方面的不足[7-8]，具有一定的優(yōu)越性。由于其動態(tài)性特點，同一觀測量無法進行重復觀測，因此測線上往往布設切割線形成縱橫交錯的網狀，根據交叉點重力值的符合程度來評估某個架次的測量質量[9]，直至評價整個測區(qū)的測量精度[10-15]。不僅如此，測線間存在的重力場水平差異通過與切割線上的交叉點重力差值進行水平調整，是數據處理中的關鍵性步驟[16]，對數據處理和地質解釋具有重要的意義[17]。交叉點不符值的搜索除用于航空重力外，還可用于船載和其他機載的地球物理場測量中，甚至衛(wèi)星測高中的平叉處理也涉及軌道交叉點的獲取[18-19]，研究該方法具有十分重要的意義。

航空重力測網設計密集且數據采樣率高，整體數據量龐大，如果采用逐一遍歷或局部遍歷搜索法，求取交叉點耗時長、效率低。目前常用的交叉點搜索方法主要是跳躍搜索法[20]、主副測線斜率法[21]、坐標平均法[22]、多項式擬合法[23]、滑動窗口法[24-25]及四叉樹分解法[26]等。其中，跳躍搜索法、主副測線斜率法和坐標平均法，都是根據測網計算近似交叉點，再設置一定范圍內進行遍歷搜索出實際交叉點。多項式擬合法通過擬合測網飛行軌跡后計算近似交叉點，再設置搜索范圍后遍歷搜索，準確率比前幾種方法有所提高，但搜索時間變長了。滑動窗口法對上述方法進行了改進，以一定搜索范圍在測線上進行滑動遍歷搜索，相交到切割線后確定搜索中心再以較小的搜索半徑進行遍歷搜索，準確率獲得了保障，但沒有提高搜索效率。四叉樹分解法，通過矩形窗口在測線和切割線上滑動遍歷鎖定相交矩形范圍后，再以四叉樹結構分解重合區(qū)域，自適應加密剖分出遍歷范圍，是一次提升遍歷搜索效率的嘗試，但對于大型測網交叉點的搜索效率，仍無法突破瓶頸。

針對以上幾種搜索方法的不足，本文提出一種交叉點的非遍歷逼近方法，不需對交叉點劃出遍歷范圍，直接迭代逼近至交叉點。試驗表明，本文方法優(yōu)化了搜索過程，還滿足了對不規(guī)則、非常規(guī)測網交叉點的搜索，提高了效率和準確率。

1 交叉點搜索方法

1.1 測網交叉點分析

航空重力測量布置測線和切割線相交叉的網狀測網，搜索交叉點就是獲取它們相交于水平面位置的坐標、場值差或其他物理量差值。在航空重力測量作業(yè)中，飛機基本保持直線飛行，受風向、人為控制不當等原因對飛行產生擾動而發(fā)生偏離航跡稱為偏航。因此，測網交叉點搜索方法設計時盡量考慮到極端偏航情況。

如圖1所示，假設航跡發(fā)生嚴重的偏航現象，A為偏航距。圖1模擬給出了測線line 1—line 5與切割線tie 1可能存在的5種位置情況：①line 1與tie 1正常相交；②line 2與tie 1真實相交，但沒有被切割線兩個端點相交到；③line 3與tie 1不相交，但端點延伸后被相交；④line 4與tie 1不相交；⑤line 5被tie 1的端點相交到，但實際沒有相交。示意圖中，實際相交的只有測線line 1、line 2，測線line 3、line 5也被誤相交到，后續(xù)需剔除。搜索測網交叉點的核心就是判斷測線與切割線是否真正相交，如果相交則獲取二者相交的準確坐標。圖1中虛線外擴端點，就是為了解決因偏航可能被遺漏的實際交叉點，把所有測線和切割線端點延伸，延長量為規(guī)范中偏航閾值80 m[27]的2～3倍，確保交叉點無遺漏。

圖1 測網中存在多種相交位置Fig.1 Multiple intersections at the boundary of the survey network

1.2 交叉點搜索思想

本文提出按照快速逼近和微調逼近組合思路，進行迭代式逼近交叉點的搜索方法設計算法，以最小計算量，快速準確定位交叉點。以每條測線和切割線上的平均點距為標尺進行2次或3次迭代式快速逼近后，測網中絕大部分近似交叉點無限逼近目標，只需小范圍內微調逼近1～2次便能搜索到交叉點。而個別交叉點處在不均勻、不連續(xù)的測網上，平均點距誤差會加大，甚至交叉點位于缺失數據的不連續(xù)測線上，所以在經過多次快速逼近后還沒有落到目標測點附近時，為避免程序陷入死循環(huán)也進行微調逼近。通過組合設計既保證交叉點搜索速度，同時兼顧精準度。

本文方法對測點間距適宜、分布均勻、測網軌跡規(guī)整的測網適用性好，搜索效率越高；對測網的測點過密、不均勻，軌跡不規(guī)整等，將增加微調逼近次數，搜索效率略低。本文搜索效率統(tǒng)計使用逼近次數概念，如快速逼近交叉點1次代表對測線line和切割線tie各進行了1次跳躍逼近，共計2次。為提高程序執(zhí)行效率，只對第1次和第2次快速逼近交叉點完成后進行判斷和逼近統(tǒng)計。若不相交則進入微調逼近循環(huán)，只要在line或tie線上逼近1次就累計1次，直至相交結束。

1.3 交叉點相交判斷方法

判斷線段是否相交包括行列式、投影法和面積法等。文獻[25]研究表明，3種方法的計算效率相當。本文使用了向量叉乘面積法，利用兩向量叉乘等于兩向量構成的平行四邊形(以兩向量的鄰邊)的面積，由于向量具有方向，所以面積也具有方向，構成三角形的面積以逆時針為正，順時針為負，令兩線段的4個端點坐標為T1(x1,y1)、T2(x2,y2)、L1(x3,y3)、L2(x4,y4)，其中x1-x4為經度，y1-y4為緯度，如圖2所示。

圖2 兩線段的向量叉乘Fig.2 Vector cross product of two line segments

圖2中兩線段4個端點構成4個三角形：L1L2T1、L1L2T2、T1T2L1和T1T2L2，它們的叉乘面積按式(1)—式(4)計算

(1)

(2)

(3)

(4)

滿足SL1L2T1·SL1L2T2≤0，切割線端點T1、T2被測線線段相交；不滿足，則位于測線同側。同理ST1T2L1·ST1T2L2≤0，測線端點L1、L2被切割線線段相交；不滿足，則位于切割線同側。當式(1)—式(4)同時滿足上述要求時，即存在交叉點。當切割線或測線的1個端點與測線或切割線兩端點處在同一條線上，此時共線下的兩向量叉乘面積為0，也存在交叉點。

符合交叉點存在條件后，計算交叉點位于切割線兩端點之間的位置系數α為

(5)

通過式(5)可獲得交叉點坐標P(x,y)，經度x、緯度y分別為

x=x1+α(x2-x1)

(6)

y=y1+α(y2-y1)

(7)

2 交叉點計算的程序實現

2.1 程序設計

按照上述搜索方法，設計了交叉點搜索流程，如圖3所示。裝載測網后，對測線和切割線延伸的4個遠端點組成線段后判斷是否相交，存在相交后第1次獲取近似交叉點坐標，再計算測線和切割線起始測點到近似交叉點的距離，根據各自的平均點距第1次快速逼近找出4個端點；若不相交接著把這4個端點再次進行線段相交，第2次獲取近似交叉點坐標，重復上述過程第2次快速逼近找出4個端點；若不相交，程序進入微調逼近循環(huán)直至相交。

圖3流程中，以微調逼近循環(huán)左側為例，先以tie線段為基準，判斷測線的相鄰端點L1、L2是否與之相交，若不相交測線兩個端點逐點向第2次近似交叉點坐標方向滑動，直至與tie線段相交；接著切割線的端點T1、T2，若不與line線段相交，兩端點逐點向最新的近似交叉點坐標方向滑動，直至與line線段相交；因切割線端點的滑動，為防止發(fā)生測線端點與tie線段又不相交，再次重復上述過程，直至真正相交。而圖3中的微調逼近循環(huán)右側，則用于在快速逼近后測線的兩端點已經相交于tie線段，但切測線的兩端點與line線段不相交的情況，實現過程與左側類同。在交叉點相交后，再搜索下一條測線或切割線，直至整個測網的交叉點搜索完成。程序通過快速逼近和微調逼近的組合進行迭代式逼近，搜索方法思路清晰，程序比較容易實現。

圖3 交叉點搜索流程Fig.3 Flow chart of intersections searching

2.2 實例分析

2.2.1 美國大地測量局EN01測網

2013年3月美國大地測量局(National Geodetic Survey,NGS)對外發(fā)布航空重力測量EN01數據塊Beta2版，此數據塊國內學者研究較多[25-26，28]，包含測點號、采樣時間、經度、緯度和場值，共100 830個數據點。表1給出了EN01測網的設計參數，圖4為測網交叉點的分布圖。

圖4 EN01數據塊測網及交叉點Fig.4 Survey network and cross point diagram of data block EN01

表1 EN01數據塊的設計參數Tab.1 Parameters of the data block EN01

表2中給出的搜索時間對比情況，逐一遍歷及優(yōu)化逐一遍歷法和本文搜索方法的運行環(huán)境：Intel Xeon E3-1505Mv6@3.0 GHZ，內存32 GB。其他方法的搜索時間都來自文獻[25-26]，四叉樹分解法使用的硬件信息：Intel Xeon E3-1230v5@3.4 GHZ，內存8 GB；其他幾個方法使用的硬件信息：Intel Core i3-2350M@2.3 GHz，內存4 GB。讀取和整理10萬個測點用時約0.1 s，本文搜索方法比其他方法的搜索速度提升了3～4個數量級。

表2 EN01用不同搜索方法的統(tǒng)計信息Tab.2 Statistical information with different methods for EN01

2.2.2 大型測網算例

使用美國大地測量局發(fā)布的航空重力測網AS08和塔里木盆地航空重力大型拼接測網TAZB，分別包含6 135 004個測點、310個交叉點，5 466 366個測點、29 517個交叉點，利用本文方法搜索并找到全部交叉點分別用時0.031 s和0.061 s。表3列出了與成熟商業(yè)軟件Oasis Montaj之間的交叉點搜索計時情況。讀取和整理AS08和TAZB文件(X、Y坐標和場值)分別約5.3 s和4.9 s，在Oasis Montaj測網數據庫中，沒有觀察到裝載數據用時。

表3 與商業(yè)地球物理軟件處理大型測網對比Tab.3 Compared with commercial geophysical software for large survey networks

3 交叉點搜索效率分析

3.1 快速逼近設計

根據多個測網驗證，常規(guī)測網可以設計成快速逼近循環(huán)直接逼近交叉點，但交叉點處在缺失數據的不連續(xù)測網中會使程序陷入死循環(huán)，另一方面，過多使用快速逼近也會影響到搜索效率，應視測網情況合理設置快速逼近的使用頻率。

取兩個國內航空重力測網為例，鞍山本溪直升機起伏飛行測網ANB的平均點距和標準差為：16.17±0.97 m，對測網分別進行2～10倍加密和2～10倍抽稀，形成平均點距1.6～161.7 m的共19個不同密度的測網。程序中分別設置快速逼近2、3和4次結合微調逼近后，統(tǒng)計這19個測網完成交叉點搜索時所用的逼近次數，以快速逼近2次的結果進行歸一化，觀察快速逼近3、4次的執(zhí)行情況，如圖5所示。圖中對平均點距在8.1 m及以上測網，設置3、4次的逼近次數歸一化大于1，不如2次時的執(zhí)行效率高。在平均點距8.1 m以下的加密測網，設置3、4次執(zhí)行效率要好于2次，但也沒有因為設置越高，執(zhí)行效率被提高。

圖5 測網ANB快速逼近次數歸一化Fig.5 Normalized graph of fast approximation times of ANB survey network

海南島周邊固定翼飛機測網HAIN的平均點距和標準差為：31.03±0.32 m，同樣進行2～10倍加密和2～10倍抽稀形成19個測網。然后進行快速逼近2、3、4次的逼近試驗，以設置2次的結果進行歸一化，發(fā)現設置2次時在各個測網密度中的執(zhí)行效率都是最佳的，并且隨著測網變疏它們之間的歸一化結果趨于收斂，如圖6所示。HAIN測網測點距的標準差偏小，即測點間分布比較均勻，搜索結論與ANB測網出現截然相反的現象，測網加密下沒有因為增加快速逼近的次數，而提高搜索效率反而出現下降。

圖6 測網HAIN快速逼近次數歸一化Fig.6 Normalized graph of fast approximation times of HAIN survey network

對比不同情況的測網，大量的試驗都得出類似結果。程序中設置成快速逼近2次和微調逼近組合，在常規(guī)測網中就可以發(fā)揮出較好的搜索效率；如果測點距標準差偏大，同時測點距小于8 m以下的測網，宜采用快速逼近3次甚至4次，再結合微調逼近進行搜索。

3.2 EN01測網搜索分析

本文逐一遍歷法：從切割線起算，在測線上遍歷找到交叉點后及時跳到下一條測線，計算工作量大約減半，即

(M-nL)×(N-nT)/2

式中，M為所有測線的測點數和；nL為測線條數；N為所有切割線的測點數和；nT為切割線條數。

優(yōu)化逐一遍歷法：對切割線相交測線的近似交叉點坐標進行排序，切割線逐點遍歷時只需遍歷第1條測線，找到交叉點后把它排除在外，剩余測線重新排序，以此類推。避開了大量的重復遍歷，缺點是對拼接測網、重復線混疊的時候不能保證100%搜索準確率。

EN01測網的測線有83 052個測點，切割線有17 778個測點，共61個交叉點。表4中為測網交叉點的搜索效率統(tǒng)計，本文方法統(tǒng)計出每個交叉點的平均逼近次數為6.56次，比逐一遍歷方式的逼近次數減少5～6個數量級，搜索速度則提升了4～5個數量級。

表4 EN01測網交叉點的搜索效率統(tǒng)計Tab.4 Search efficiency statistics for EN01 intersections

EN01測網中有兩條測線L1122和L1124對重力擾動假異常做了切除，分別丟失319和630個測點，造成約36 km和71 km數據缺失。本文方法受這兩條測線數據的不連續(xù)影響，測線中平均點距誤差加大導致交叉點搜索效率下降。在剔除這兩條測線后，每個交叉點平均逼近次數由原先6.56次下降到4.75次，搜索效率明顯提升，詳見表5。

3.3 多個測網交叉點搜索分析

表5中，美國大地測量局使用二軸阻尼穩(wěn)定平臺TAGS航空重力儀，因抗飛行顛簸性能差，人為切除了EN01、PN02和AS08測網中部分測線的擾動異常，從而造成測線中測點的不連續(xù)性。其中AS08測網切除了7條線的擾動異常，6個交叉點又位于剪切處，采樣率為20 Hz；剔除這7條線后為AS08′測網，平均點距及標準差為6.06±0.59 m。EN01和PN02測網的采樣率均為1 Hz，EN01′和PN02′為各剔除2條不連續(xù)測線后的測網。表5、表6中，除測網AS08、AS08′和YYQC使用快速逼近3次外，其余測網均使用2次，再采用微調逼近。

表5 多個測網交叉點的搜索效率統(tǒng)計Tab.5 Search efficiency statistics for more network intersections

測網EN01、PN02，平均點距的離散程度偏大，每個交叉點的平均逼近次數搜索效率略低；AS08測網因有6個交叉點位于不連續(xù)處，搜索效率更低。而較規(guī)則測網EN01′、PN02′、ANB、HAIN、TAZB的平均逼近次數在3.9～4.75之間；平均點距小于8 m的AS08′測網，快速逼近調高至3次后的平均逼近次數為6.64，調高至4次時沒有獲得提升。

某海域船載重力測網YYQC的平均點距及標準差為2.61±1.01 m，從標準差的偏離程度看出測點距比較離散或者測線之間的航行速度控制不穩(wěn)定，測網的船測航跡保持也欠規(guī)整，如圖7所示。使用快速逼近3次后的每個交叉點平均逼近次數為21.82，調高至4次后，平均逼近次數下降至13.13，獲得明顯提升。

圖7 遠洋船載重力測網Fig.7 Ocean-going ship gravity survey network

表6中，統(tǒng)計多個測網的交叉點逼近次數2～8次，分別對應交叉點搜索的完成情況。平均點距大于8 m，較規(guī)則測網中只有大型拼接測網TAZB在逼近8次后完成98.2%，測網EN01′、PN02′、ANB和HAIN都在8次逼近內完成。表5、表6中得出認識：測網測點連續(xù)性、航跡規(guī)整、平均點距大小和它的離散程度適宜，對提高交叉點的搜索效率至關重要。

表6 多個測網的逼近次數與搜索完成比例統(tǒng)計Tab.6 Statistics of the number of approximations and completion ratios formore survey networks

4 結 論

本文提出快速逼近和微調逼近組合進行迭代式非遍歷逼近交叉點的方法，克服了常規(guī)方法進行局部遍歷的不足，大幅度提升了交叉點的搜索效率和準確率。搜索速度較現有文獻方法提升3～4個數量級，也遠優(yōu)于著名商業(yè)地球物理軟件，對多個測網的應用效果十分顯著。通過測網交叉點搜索效率分析得出：

(1) 針對常規(guī)航空重力測網，交叉點搜索流程中快速逼近只需設置2次，就能實現快速搜索。

(2) 測網中測點不連續(xù)或者平均點距小于8 m，且標準差偏大，建議快速逼近調升至3次，發(fā)揮出較好的搜索效率。

(3) 對于航行速度慢又不穩(wěn)定，以及航跡不規(guī)整的船載重力測網，為提高搜索速度，建議快速逼近設置成4次，彌補船載測網的不足。