楊金政， 邱崇濤， 田 宇， 高國(guó)林， 祁 程

(1.核工業(yè)航測(cè)遙感中心，石家莊 050002;2.中核集團(tuán)公司 鈾資源地球物理勘查技術(shù)中心(重點(diǎn)實(shí)驗(yàn)室)，石家莊 050002)

0 前言

航空地球物理勘探簡(jiǎn)稱航空物探，具有效率高、成本低、便于大面積工作、探測(cè)深度較大等優(yōu)點(diǎn)，是基礎(chǔ)性和公益性地質(zhì)調(diào)查、戰(zhàn)略性礦產(chǎn)勘查的重要手段，是地質(zhì)勘查現(xiàn)代化的標(biāo)志之一[1]。為了降低人員成本與風(fēng)險(xiǎn)，2011年，中國(guó)國(guó)土資源部試行了航空物探測(cè)量無(wú)人值守工作方式，2013年由中國(guó)核工業(yè)集團(tuán)公司和中國(guó)地質(zhì)研究院又成功進(jìn)行了無(wú)人機(jī)航空測(cè)量試驗(yàn)。無(wú)人值守的工作方式與傳統(tǒng)方式的一個(gè)顯著差異是飛機(jī)從起飛到落地，儀器使用同一個(gè)測(cè)線編號(hào)持續(xù)記錄數(shù)據(jù)，首先需要對(duì)航跡數(shù)據(jù)進(jìn)行測(cè)線識(shí)別、裁剪，目前，許多國(guó)內(nèi)、外著名的地球物理數(shù)據(jù)處理軟件(如Montaj Oasis、geoProbe等)，在此環(huán)節(jié)上，處理起來(lái)并不十分“舒暢”。筆者利用國(guó)內(nèi)地質(zhì)行業(yè)最常用的MapGis地理信息系統(tǒng)，對(duì)航測(cè)數(shù)據(jù)實(shí)現(xiàn)可視化的數(shù)據(jù)編輯與篩選，不僅適用于無(wú)人值守飛行數(shù)據(jù)，而且對(duì)于常規(guī)飛行數(shù)據(jù)，同樣具有針對(duì)性強(qiáng)、工作效率高的特點(diǎn)。

MapGis地理信息系統(tǒng)是武漢中地信息工程有限公司開(kāi)發(fā)的GIS基礎(chǔ)平臺(tái)軟件系統(tǒng)，集地圖輸入、數(shù)據(jù)庫(kù)管理及空間分析的綜合系統(tǒng)；

不僅圖形編輯功能強(qiáng)大，而且具有較完善的二次開(kāi)發(fā)功能[2]；是國(guó)內(nèi)地質(zhì)、物探領(lǐng)域中應(yīng)用最廣泛的軟件之一。屬性管理子系統(tǒng)是MapGis的重要組成部分，可實(shí)現(xiàn)對(duì)屬性數(shù)據(jù)管理[3]。本程序是利用MapGis圖元的屬性結(jié)構(gòu)(包括點(diǎn)、線和區(qū)屬性結(jié)構(gòu))，通過(guò)錄入、修改屬性數(shù)據(jù)，借助MapGis強(qiáng)大的二次開(kāi)發(fā)和圖形編輯功能，交互式地進(jìn)行數(shù)據(jù)的識(shí)別、篩選，最終獲得符合質(zhì)量要求的測(cè)量數(shù)據(jù)。

1 設(shè)計(jì)背景

1.1 航空物探數(shù)據(jù)簡(jiǎn)介

我中心承擔(dān)的航測(cè)任務(wù)中，大多以航磁、航放同步測(cè)量為主。航空數(shù)據(jù)采集是以飛行架次為單位。數(shù)據(jù)記錄為一個(gè)二進(jìn)制文件，先將二進(jìn)制文件轉(zhuǎn)換為3個(gè)ASCII文件(航磁、航放窗數(shù)據(jù)、航放全譜)；每個(gè)測(cè)量數(shù)據(jù)均包括文件頭和測(cè)量數(shù)據(jù)兩個(gè)部分；文件頭主要記錄儀器型號(hào)、參數(shù)設(shè)置、測(cè)量時(shí)間、地點(diǎn)、架次等參數(shù)。測(cè)量數(shù)據(jù)記錄了該架次所有實(shí)測(cè)數(shù)據(jù)，包括早晚測(cè)試(能譜)、早、晚基線(磁、放)、測(cè)線(磁、放)等測(cè)量數(shù)據(jù)。圖1為無(wú)人機(jī)航放窗數(shù)據(jù)ASCII文件的數(shù)據(jù)格式。從圖1中可看出，測(cè)量數(shù)據(jù)包括測(cè)線號(hào)、基點(diǎn)號(hào)、經(jīng)緯度、坐標(biāo)值、高度值、溫度值、放射性元素測(cè)量值、日期和時(shí)間等內(nèi)容。

1.2 數(shù)據(jù)質(zhì)量控制

航測(cè)數(shù)據(jù)的采集受諸多因素的影響，在確保天氣、儀器性能、地面磁日變站等正常前提下，通過(guò)測(cè)區(qū)范圍、飛行速度、離地飛行高度、偏航距等指標(biāo)衡量數(shù)據(jù)質(zhì)量；在物探項(xiàng)目中，項(xiàng)目設(shè)計(jì)測(cè)線是野外數(shù)據(jù)采集最基本的依據(jù)。離地飛行高度(離地高度)，在航空伽瑪能譜測(cè)量中是一項(xiàng)十分重要的參數(shù)。一般情況下，飛機(jī)應(yīng)沿地形起伏飛行，高度保持在80 m～ 120 m。由于放射性計(jì)數(shù)率隨高度呈指數(shù)衰減，離地高度過(guò)高，導(dǎo)致計(jì)數(shù)率偏低，將會(huì)使地面弱異常被遺漏，同時(shí)增大統(tǒng)計(jì)漲落誤差；離地高度過(guò)低，無(wú)法將測(cè)線間區(qū)域全覆蓋[4]，同時(shí)飛行安全指數(shù)大大降低。偏航距指實(shí)際飛行航跡與設(shè)計(jì)測(cè)線的偏移距離，依據(jù)測(cè)量比例尺確定最大偏航距[5]。在實(shí)際飛行測(cè)量中，一條測(cè)線結(jié)束后，飛機(jī)轉(zhuǎn)入下一條測(cè)線時(shí)，測(cè)量數(shù)據(jù)雖處于項(xiàng)目設(shè)計(jì)之外，但在實(shí)際數(shù)據(jù)處理時(shí)，應(yīng)盡量多地保留測(cè)區(qū)外符合質(zhì)量要求的數(shù)據(jù)。在測(cè)區(qū)中，受高大建筑物、地形等因素影響，往往也會(huì)出現(xiàn)“超高”、“偏離”的現(xiàn)象，此時(shí)必須及時(shí)了解具體情況，做出下一步飛行方案。

2 程序?qū)崿F(xiàn)

為了盡可能地保留測(cè)區(qū)外的有效數(shù)據(jù)，同時(shí)及時(shí)準(zhǔn)確掌握各架次的飛行質(zhì)量，可視化的數(shù)據(jù)篩選顯得十分必要。

2.1 設(shè)計(jì)流程

總體上，利用此程序進(jìn)行數(shù)據(jù)編輯遵循“數(shù)據(jù)—圖形—數(shù)據(jù)”變換過(guò)程，分為5個(gè)步驟，除第4步為手工操作外，其余均由程序自動(dòng)完成(圖2)。

1)設(shè)計(jì)測(cè)線繪制。根據(jù)設(shè)計(jì)坐標(biāo)，生成具有屬性結(jié)構(gòu)的點(diǎn)、線、區(qū)三個(gè)文件，將各自設(shè)計(jì)測(cè)線信息錄入其中，作為自動(dòng)識(shí)別依據(jù)。

圖2 程序設(shè)計(jì)流程圖Fig.2 The flow chart of program

圖3 無(wú)人機(jī)某一架次物探作業(yè)航跡圖Fig.3 Map of a sortie of flight track of airborne geophysical survey using drone

2)航跡線的繪制。將實(shí)測(cè)數(shù)據(jù)繪制成圖，線文件為航跡線，包含“測(cè)線編號(hào)”的屬性。點(diǎn)文件包括：①子圖點(diǎn)文件，表明“測(cè)點(diǎn)位置”和“高度等級(jí)”；②注釋點(diǎn)文件，包括“測(cè)線編號(hào)”、“飛行高度”信息(圖3)。

3)依據(jù)設(shè)計(jì)測(cè)線，將航跡線自動(dòng)分離、識(shí)別并補(bǔ)充錄入“偏航距”信息。

4)在MapGis中，依據(jù)測(cè)點(diǎn)信息(點(diǎn)文件)，手工編輯航跡線，剔除不合格數(shù)據(jù)。

5)依據(jù)編輯后的航跡線端點(diǎn)坐標(biāo)，篩選實(shí)測(cè)數(shù)據(jù)。

6)換名存儲(chǔ)航跡線，作為下一架次數(shù)據(jù)編輯的參照。

2.2 程序?qū)崿F(xiàn)

此程序基于MapGis軟件提供的COM組件[6]，使用C#語(yǔ)言編寫(xiě)。有關(guān)MapGis的二次開(kāi)發(fā)環(huán)境、組件注冊(cè)與引用以及點(diǎn)、線、區(qū)的寫(xiě)入等技術(shù)[7-8]。筆者以2015年在新疆喀什地區(qū)開(kāi)展的無(wú)人機(jī)航空物探綜合站測(cè)量技術(shù)研發(fā)與應(yīng)用示范項(xiàng)目[9-10]中的航測(cè)數(shù)據(jù)為依據(jù)，得到簡(jiǎn)述程序?qū)崿F(xiàn)過(guò)程。

2. 2.1 設(shè)計(jì)測(cè)線成圖

1)準(zhǔn)備工作。在程序運(yùn)行前，建立帶有屬性結(jié)構(gòu)的點(diǎn)、線、區(qū)模板文件。由于該程序設(shè)計(jì)僅針對(duì)于飛行質(zhì)量，并不對(duì)測(cè)量值進(jìn)行修正、處理，所以線、區(qū)文件只添加“測(cè)線編號(hào)”屬性；點(diǎn)文件添加“基點(diǎn)號(hào)”和“測(cè)線編號(hào)”兩項(xiàng)屬性(圖4)，其他數(shù)值未寫(xiě)入屬性結(jié)構(gòu)中。上述三個(gè)模板文件均為空文件，隨后生成的MapGis文件均基于此模板完成[3]。設(shè)置方法可參見(jiàn)參考文獻(xiàn)[3]。

為提高運(yùn)算效率，坐標(biāo)位置也未寫(xiě)入屬性結(jié)構(gòu)中，而是通過(guò)圖中繪制的測(cè)線、測(cè)點(diǎn)的坐標(biāo)值換算獲得。

2)設(shè)計(jì)測(cè)線的繪制。根據(jù)基線、設(shè)計(jì)測(cè)線的端點(diǎn)坐標(biāo)，在MapGis中繪制設(shè)計(jì)測(cè)線圖(圖5)，包括測(cè)線編號(hào)(點(diǎn)文件)、設(shè)計(jì)測(cè)線(線文件)、測(cè)線識(shí)別范圍(區(qū)文件)。測(cè)線識(shí)別范圍是指單條設(shè)計(jì)測(cè)線兩側(cè)外擴(kuò)接近1/2線距，兩端外擴(kuò)1 km～ 2 km的區(qū)域。由于線、區(qū)文件均是基于上述模板文件生成，所以它們中的圖元均具有“測(cè)線編號(hào)”的屬性結(jié)構(gòu)。線文件，屬性錄入代碼為：

myLinArea.Load(modelWL);//加載線模板文件

Record att2 = new Record(); //屬性定義并實(shí)例化

myLinArea.att.Get(LineCounts, out att2); //獲取原有屬性

att2.Value[2] = LineNo; //獲取屬性(線號(hào))

myLinArea.att.Write(LineCounts, att2);//寫(xiě)入屬性

……

//另存線文件

myLinArea.Save(WorkArea23 + @"設(shè)計(jì)測(cè)線.wl");

//區(qū)屬性定義與線屬性類似，(略)

另外，基于測(cè)區(qū)范圍外1 km～2 km，小于飛機(jī)轉(zhuǎn)彎區(qū)繪制一個(gè)范圍框，稱之為“航跡分離框”，用于實(shí)測(cè)測(cè)線的分離(圖4)。

圖4 點(diǎn)、線、區(qū)文件屬性結(jié)構(gòu)參數(shù)Fig.4 Parameters of attribute structure of point, line and region files(a)點(diǎn)文件屬性結(jié)構(gòu);(b)線文件屬性結(jié)構(gòu);(c)區(qū)文件屬性結(jié)

圖5 在MapGis中航跡編輯界面圖Fig.5 The edit interface of flight track on MapGis

2.2.2 實(shí)測(cè)數(shù)據(jù)成圖

本次測(cè)量中，航磁測(cè)量采樣率為10 Hz，航放測(cè)量為1 Hz[9]。為了提高運(yùn)算效率，優(yōu)先選用數(shù)據(jù)量較小的能譜數(shù)據(jù)，利用線號(hào)、基點(diǎn)號(hào)、坐標(biāo)位置(x，y)、離地高度等數(shù)據(jù)列，繪制具有屬性結(jié)構(gòu)的測(cè)點(diǎn)(子圖)、注釋點(diǎn)、航跡(線)(圖3、圖5)。依據(jù)離地高度等級(jí)，將測(cè)點(diǎn)和注釋點(diǎn)賦予不同的顏色，直觀顯示數(shù)據(jù)質(zhì)量。測(cè)點(diǎn)及注釋點(diǎn)的繪制代碼如下：

private void writingCedian(PntArea PntAreaCedian, PntArea PntAreaCDnote, DotInfo5 myDot2, Pnt_Info PntInfoNor, Pnt_Info PntInfoNote) //測(cè)點(diǎn)、注釋點(diǎn)子函數(shù)

{

double fontSize = Convert.ToDouble(textBox11.Text);

Record att3 = new Record();//屬性定義并實(shí)例化

D_Dot pntXY3cd = new D_Dot();

pntXY3cd.x = myDot2.x; //測(cè)點(diǎn)X坐標(biāo)

pntXY3cd.y = myDot2.y; //測(cè)點(diǎn)Y坐標(biāo)

D_Dot pntXY3cdAlt = new D_Dot();

pntXY3cdAlt.x = myDot2.x+1; //注釋點(diǎn)X坐標(biāo)

pntXY3cdAlt.y = myDot2.y - fontSize*0.75; //注釋點(diǎn)X坐標(biāo)

PntAreaCedian.Append(pntXY3cd, "a", PntInfoNor);

int DotCounts = PntAreaCedian.count - 1;

PntAreaCedian.att.Get(DotCounts, out att3); //取點(diǎn)屬性

att3.Value[1] = myDot2.fn;//基點(diǎn)號(hào)賦值

att3.Value[2] = myDot2.lineNo;//測(cè)線號(hào)賦值

att3.Value[3] = myDot2.ralt;//雷達(dá)高度號(hào)賦值

PntAreaCedian.att.Write(LineCounts, att3);//寫(xiě)入屬性值

//寫(xiě)入測(cè)點(diǎn)

PntAreaCDnote.Append(pntXY3cd,myDot2.fn,PntInfoNote);

//寫(xiě)入注釋點(diǎn),高度值,偏航距為零,測(cè)線識(shí)別后,再寫(xiě)入

PntAreaCDnote.Append(pntXY3cdAlt,myDot2.ralt.ToString() + "-" + myDot2.dis, PntInfoNoteNor);

}

程序中myDot2、PntInfoNor和PntInfoNote分別為測(cè)點(diǎn)、注釋點(diǎn)的參數(shù)值，以不同的點(diǎn)參數(shù)值標(biāo)示飛行質(zhì)量等級(jí)。

2.2.3 測(cè)線識(shí)別

首先，利用航跡分離框?qū)⒈炯艽魏桔E線分離為獨(dú)立的航跡線段(圖5)。代碼如下：

double m_Radiu = 0.0001;//模糊半徑

MAPGISBASCOM1Lib.IAnalysis pClip = new MAPGISBASCOM1Lib.Analysis();

//myFrame-分離框rawLinArea-原始線文件desLinArea-結(jié)果線文件

pClip.ClipLin(myFrame, rawLinArea, desLinArea, m_Radiu, Enum_Clip_Type.gisOVLY_INCLIP);

然后，遍歷各航跡線段，利用測(cè)線識(shí)別范圍(區(qū)文件)，逐一賦予“測(cè)線編號(hào)”屬性。傳遞屬性關(guān)鍵代碼如下：

for (int intK = 1; intK < myLinArea.count; intK++)

{ //遍歷所有實(shí)測(cè)測(cè)線

//獲取實(shí)測(cè)測(cè)線上的點(diǎn)集合

int i = myLinArea.Get(intK, out dataSet, out inf, out dima);

…….(在測(cè)線尋找3個(gè)合適的點(diǎn)坐標(biāo)，代碼較長(zhǎng)，略)

myLinArea.att.Get(intK, out att2); //獲取實(shí)測(cè)測(cè)線屬性

for (int intN = 1; intN < myRegArea.count; intN++)

{ //遍歷所有設(shè)計(jì)測(cè)線范圍

//計(jì)算實(shí)測(cè)測(cè)點(diǎn)與設(shè)計(jì)測(cè)線范圍的最大距離

double minDis1 = myRegArea.MaxDistOfPntToReg(dot1,intN); double minDis2 =myRegArea.MaxDistOfPntToReg(dot2,intN);

if (minDis1 == 0 && minDis2 == 0)

{ //最大距離為0時(shí)，表明測(cè)線在設(shè)計(jì)測(cè)線范圍內(nèi)

//獲取設(shè)計(jì)測(cè)線范圍的屬性，即線號(hào)

myRegArea.RegAtt.Get(intN, out myRecord);

string LineNo = myRecord.Value[3]; //賦值

att2.Value[2] = LineNo; // 傳遞給線屬性(線號(hào))

myLinArea.att.Write(intK, att2); // 錄入線屬性

break;

}

}

}

最后，根據(jù)設(shè)計(jì)測(cè)線位置，計(jì)算每個(gè)航跡點(diǎn)的偏航距并繪制更新偏航距的標(biāo)示(圖5)。

2.2.4 編輯航跡

借助MapGis軟件平臺(tái)，使用不同的顏色、符號(hào)標(biāo)示出每個(gè)測(cè)點(diǎn)數(shù)據(jù)質(zhì)量等級(jí)，將實(shí)測(cè)數(shù)據(jù)以圖形展現(xiàn)，可以很清晰了解本架次的數(shù)據(jù)質(zhì)量。參照以往架次的航跡線和測(cè)點(diǎn)標(biāo)識(shí)，利用MapGis圖形編輯功能，逐一修剪本架次航跡線，確定需保留數(shù)據(jù)范圍。

2.2.5 數(shù)據(jù)篩選與輸出

數(shù)據(jù)編輯基于MapGis中確定的航跡線端點(diǎn)，從源數(shù)據(jù)文件中提取質(zhì)量合格的數(shù)據(jù)。

讀取每條航跡線的端點(diǎn)坐標(biāo)和航跡點(diǎn)(點(diǎn)文件)屬性結(jié)構(gòu)，從航跡點(diǎn)中獲取“測(cè)線編號(hào)”和“基點(diǎn)號(hào)”。以“基點(diǎn)號(hào)”為依據(jù)，比對(duì)、篩選航磁、航放源數(shù)據(jù)，更換“測(cè)線編號(hào)”，獲得最終數(shù)據(jù)。從航跡線中提取信息的代碼如下：

for (int intL = 1; intL < myLinArea3.count; intL++)

{

//獲取圖元位置

myLinArea3.Get(intL, out lineSet, out linA, out dim);

myLinArea3.att.Get(intL, out attLine); //獲取圖元屬性

for (int intI = 0; intI < lineSet.count; intI++)

{ //遍歷線圖元節(jié)點(diǎn)

tempXY.x = Convert.ToDouble (lineSet[intI].x.ToString ("0.00"));//提取X坐標(biāo)

tempXY.y = Convert.ToDouble (lineSet[intI].y.ToString ("0.00"));//提取Y坐標(biāo)

tempXY.lineNo =attLine.Value[2]; //提取線號(hào)

lineLinInfo.Add(tempXY); //寫(xiě)入泛型集合

}

}

從航跡點(diǎn)中提取信息的代碼如下：

for (int intK = 1; intK < myPntArea3.count; intK++)

{ //遍歷點(diǎn)圖元

myPntArea3.att.Get(intK, out myAtt);//獲取圖元屬性

myPntArea3.GetPos(intK, out myXY); //獲取圖元位置

//提取X坐標(biāo)

myDot1.x = Convert.ToDouble(myXY.x.ToString("0.00"));

//提取Y坐標(biāo)

myDot1.y = Convert.ToDouble(myXY.y.ToString("0.00"));

myDot1.lineNo =myAtt.Value[2]; //提取線號(hào)

myDot1.Fn =myAtt.Value[1]; //提取基點(diǎn)號(hào)

//提取離地高度值

myDot1.ralt = Convert.ToDouble(myAtt.Value[3]);

myDot2.Add(myDot1); //寫(xiě)入泛型集合

}

圖6為數(shù)據(jù)篩選后的結(jié)果數(shù)據(jù)。從圖6中可以看出，數(shù)據(jù)格式未變，測(cè)線編號(hào)已將原來(lái)的“005”變?yōu)榱藢?shí)際測(cè)線號(hào)，并刪除了不符合質(zhì)量要求的數(shù)據(jù)。

另外，該程序還進(jìn)行了本架次的數(shù)據(jù)質(zhì)量、工作量等信息計(jì)算、統(tǒng)計(jì)。

將本架次的航跡線更換不同顏色，以飛行日期和架次命名保存此次航跡線，供下一架次的數(shù)據(jù)裁截做參照對(duì)比。

2.3 有人值守航測(cè)數(shù)據(jù)的處理

該程序同樣適用于有人值守航測(cè)數(shù)據(jù)篩選。在正常情況下，由于每個(gè)架次的測(cè)線編號(hào)已事先錄入，可直接手工編輯航跡線，再根據(jù)線號(hào)、基點(diǎn)號(hào)比對(duì)，篩選實(shí)測(cè)數(shù)據(jù)。如若在儀器發(fā)生故障、操作員記錄錯(cuò)誤等情況下，容易圈定錯(cuò)誤范圍，利用航跡線屬性修改、圖形編輯篩選數(shù)據(jù)，可明顯提高工作效率。

3 結(jié)語(yǔ)

MapGis軟件是國(guó)內(nèi)地質(zhì)行業(yè)使用范圍最廣的軟件之一，基于其屬性結(jié)構(gòu)進(jìn)行二次開(kāi)發(fā)，具有開(kāi)發(fā)效率高，操作簡(jiǎn)單快捷等優(yōu)點(diǎn)。利用程序提供的圖形界面，交互式地編輯圖形，獲取合格的實(shí)測(cè)數(shù)據(jù)，不僅提高了數(shù)據(jù)裁切的針對(duì)性，而且可以熟知飛行數(shù)據(jù)質(zhì)量和工作進(jìn)度。同時(shí)，在數(shù)據(jù)截取后，對(duì)飛行數(shù)據(jù)質(zhì)量進(jìn)行統(tǒng)計(jì)，避免了重復(fù)工作，提高了工效。此程序已應(yīng)用于新疆喀什地區(qū)無(wú)人機(jī)試驗(yàn)項(xiàng)目(2015年)、黑龍江省完達(dá)山—太平嶺地區(qū)航測(cè)項(xiàng)目(2016年)和黑龍江省龍江—嫩江地區(qū)航測(cè)項(xiàng)目(2016年)中，經(jīng)不斷完善，實(shí)用效果較好。