航空重力測網(wǎng)交叉點的非遍歷逼近方法

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

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

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

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

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

1 交叉點搜索方法

1.1 測網(wǎng)交叉點分析

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

如圖1所示，假設(shè)航跡發(fā)生嚴重的偏航現(xiàn)象，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ù)需剔除。搜索測網(wǎng)交叉點的核心就是判斷測線與切割線是否真正相交，如果相交則獲取二者相交的準確坐標。圖1中虛線外擴端點，就是為了解決因偏航可能被遺漏的實際交叉點，把所有測線和切割線端點延伸，延長量為規(guī)范中偏航閾值80 m[27]的2～3倍，確保交叉點無遺漏。

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

1.2 交叉點搜索思想

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

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

1.3 交叉點相交判斷方法

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

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

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

(1)

(2)

(3)

(4)

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

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

(5)

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

x=x1+α(x2-x1)

(6)

y=y1+α(y2-y1)

(7)

2 交叉點計算的程序?qū)崿F(xiàn)

2.1 程序設(shè)計

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

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

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

2.2 實例分析

2.2.1 美國大地測量局EN01測網(wǎng)

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

圖4 EN01數(shù)據(jù)塊測網(wǎng)及交叉點Fig.4 Survey network and cross point diagram of data block EN01

表1 EN01數(shù)據(jù)塊的設(shè)計參數(shù)Tab.1 Parameters of the data block EN01

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

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

2.2.2 大型測網(wǎng)算例

使用美國大地測量局發(fā)布的航空重力測網(wǎng)AS08和塔里木盆地航空重力大型拼接測網(wǎng)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測網(wǎng)數(shù)據(jù)庫中，沒有觀察到裝載數(shù)據(jù)用時。

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

3 交叉點搜索效率分析

3.1 快速逼近設(shè)計

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

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

圖5 測網(wǎng)ANB快速逼近次數(shù)歸一化Fig.5 Normalized graph of fast approximation times of ANB survey network

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

圖6 測網(wǎng)HAIN快速逼近次數(shù)歸一化Fig.6 Normalized graph of fast approximation times of HAIN survey network

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

3.2 EN01測網(wǎng)搜索分析

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

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

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

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

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

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

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

3.3 多個測網(wǎng)交叉點搜索分析

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

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

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

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

圖7 遠洋船載重力測網(wǎng)Fig.7 Ocean-going ship gravity survey network

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

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

4 結(jié) 論

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

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

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

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