基于Python 編程的航空伽馬全譜分析方法的應(yīng)用研究

馮春圓，王會(huì)波，周子陽，蔡文軍，楊春雷，朱林

（1.核工業(yè)航測遙感中心 河北省航空探測與遙感技術(shù)重點(diǎn)實(shí)驗(yàn)室 中核鈾資源地球物理勘查技術(shù)中心（重點(diǎn)實(shí)驗(yàn)室），河北 石家莊 050002；2.中核地礦科技集團(tuán)有限公司，北京100013；3.華業(yè)（北京）科技有限公司，北京100013；4.中國冶金地質(zhì)總局一局，河北 三河 065201）

傳統(tǒng)的航空伽馬能譜數(shù)據(jù)處理方法主要為國際原子能機(jī)構(gòu)（1991）［1］推薦的“標(biāo)準(zhǔn)三窗法”，簡稱“三窗法”（TWM），即所謂的窗數(shù)據(jù)處理法，這種方法在地表勘查和淺部鈾礦床勘探過程中曾發(fā)揮了重要作用，而且仍在地質(zhì)礦產(chǎn)勘查等許多領(lǐng)域得到較為廣泛的應(yīng)用。該方法也是眾多數(shù)據(jù)處理軟件常用的數(shù)據(jù)處理方法，如oasis montaj、AGRS GeoProbe 等軟件，均提供標(biāo)準(zhǔn)三窗法數(shù)據(jù)處理模塊。隨著航空伽馬能譜測量向放射性本底調(diào)查、環(huán)境測量、核事故應(yīng)急航空伽馬能譜測量、地下水及溫泉調(diào)查以及大氣氡評價(jià)和修正等方面的應(yīng)用［2-9］。數(shù)據(jù)處理方法也從單一的窗數(shù)據(jù)處理走向譜數(shù)據(jù)處理發(fā)展。

航空伽馬全能譜分析方法，簡稱“全譜法”（FSA），是基于全能譜數(shù)據(jù)進(jìn)行的一種數(shù)據(jù)處理方法。就是利用航空伽馬能譜測量中所有能譜數(shù)據(jù)，而不是常規(guī)的窗數(shù)據(jù)，直接進(jìn)行飛機(jī)和儀器本底、宇宙輻射本底、航測高度等的各項(xiàng)修正，而將大氣氡視為變量，與地面輻射源的鉀、鈾、釷元素含量以及其他核污染源（人工核素）的濃度或活度一起進(jìn)行計(jì)算的分析方法［10］。該方法能有效解決窗數(shù)據(jù)處理中大氣氡修正引入較大不確定度問題，前人已在原理和應(yīng)用方面作了大量的研究［11-14］。核工業(yè)航測遙感中心于2003年在Windows XP 操作系統(tǒng)上，基于Microsoft.NET Framework 模塊通過C#調(diào)用MATLAB軟件進(jìn)行了針對航空伽馬全能譜分析方法的編程。受限于MATLAB 與Microsoft.NET Framework的通信、MATLAB與C#通信過程中的數(shù)據(jù)轉(zhuǎn)換等影響，以致軟件整體的運(yùn)行速度較慢、可移植性和兼容性較低。

Python 是一種面向?qū)ο蟮慕忉屝缘目缙脚_高級計(jì)算機(jī)語言，具有語法簡單、可移植性強(qiáng)的特點(diǎn)。同時(shí)，它是一種通用性語言并具有極為豐富的類庫，可應(yīng)用于數(shù)據(jù)庫、多媒體、科學(xué)計(jì)算、網(wǎng)絡(luò)、游戲等諸多領(lǐng)域。

作者從航空伽馬全能譜分析方法的原理出發(fā)，以計(jì)算天然放射性核素為例，利用Python 的科學(xué)計(jì)算庫，實(shí)現(xiàn)了航空伽馬全能譜分析方法計(jì)算天然放射性核素。通過對云南昆明某試驗(yàn)區(qū)航空伽馬能譜測量數(shù)據(jù)的處理，發(fā)現(xiàn)利用Python 編程處理數(shù)據(jù)，全過程無需調(diào)用MATLAB、SCILAB 等大型數(shù)值計(jì)算軟件，即可實(shí)現(xiàn)了航空伽馬全能譜分析方法的全過程，數(shù)據(jù)處理方法可行、結(jié)果可靠，處理效率大幅提升，為航空伽馬全能譜分析方法計(jì)算以及其他全譜數(shù)據(jù)處理提供了一種可實(shí)現(xiàn)的方式。

1 航空伽馬全能譜分析方法原理

對于以天然放射性核素為目的的航空伽馬能譜測量，測量源項(xiàng)包括：宇宙射線、飛機(jī)本底、大氣氡及陸地天然放射性核素40K、鈾系和釷系等。測量能譜的計(jì)數(shù)率可由下式表示：

式中：N—一個(gè)（1×256）觀測能譜計(jì)數(shù)率數(shù)據(jù)矩陣；B—一個(gè)（1×256）飛機(jī)和儀器輻射本底矩陣；C—一個(gè)（1×256）宇宙射線本底矩陣；Q—一個(gè)（1×4）放射性核素（天然放射性K、U、Th 及大氣氡Rn）的含量（或質(zhì)量活度）矩陣；M—一個(gè)（256×4）單位能譜矩陣（又稱靈敏度矩陣），MT為M的轉(zhuǎn)置矩陣，M-T為MT的逆矩陣。其中，計(jì)數(shù)率數(shù)據(jù)矩陣N為實(shí)際測量值，屬于已知項(xiàng)，飛機(jī)和儀器輻射本底矩陣B、宇宙射線本底矩陣C、單位能譜矩陣U為修正系數(shù)，需在實(shí)際飛行前通過高高度測量、動(dòng)態(tài)帶飛行測量、航模模型刻度等獲取，各測量要求詳見航空伽馬能譜測量規(guī)范［15］，矩陣Q為天然放射性K、U、Th 及大氣氡（Rn）含量或質(zhì)量活度的待求矩陣。

航空伽馬全能譜分析方法是基于全譜來分析計(jì)算，在譜分析之前須檢查實(shí)測譜數(shù)據(jù)和雷達(dá)、溫度等數(shù)據(jù)，以確保用于全能譜分析方法計(jì)算的準(zhǔn)確性。主要的檢驗(yàn)內(nèi)容有數(shù)據(jù)檢驗(yàn)、峰漂修正等。數(shù)據(jù)檢驗(yàn)包括雷達(dá)、溫度、氣壓等數(shù)據(jù)是否完整、正常，譜數(shù)據(jù)是否有峰漂現(xiàn)象等；峰漂修正是針對有峰位漂移的譜進(jìn)行的一種修正，使得全譜信息更符合實(shí)際測量情況。

2 航空伽馬全能譜分析的程序設(shè)計(jì)

根據(jù)航空伽馬全能譜分析方法的原理，程序?qū)崿F(xiàn)過程分為數(shù)據(jù)讀取與整理、預(yù)處理、擬合計(jì)算、結(jié)果輸出4 個(gè)部分（圖1）。

圖1 航空伽馬全譜分析法Python 編程設(shè)計(jì)簡圖Fig.1 The design diagram of Python programming for airborne gamma-ray full spectrum analysis

數(shù)據(jù)讀取與整理包括本底譜數(shù)據(jù)、單位譜（K、U、Th、Rn）數(shù)據(jù)以及實(shí)測譜數(shù)據(jù)的讀取，并對實(shí)測譜數(shù)據(jù)進(jìn)行初步數(shù)據(jù)檢查。檢查內(nèi)容主要包括雷達(dá)、溫度、氣壓數(shù)據(jù)等是否齊全和可靠，如有跳點(diǎn)、假值等情況時(shí)需進(jìn)行數(shù)據(jù)剔除和插值，尤其是雷達(dá)數(shù)據(jù)的連續(xù)性較差時(shí)，可采用5 點(diǎn)平滑濾波等方法進(jìn)行處理。預(yù)處理主要包括死時(shí)間修正、能譜濾波、峰漂修正、計(jì)算STP 等效高度、本底扣除等，最后將飛機(jī)儀器本底譜、宇宙射線本底譜和單位譜數(shù)據(jù)分別整理成矩陣B、C、M（K、U、Th、Rn 單位譜及修正系數(shù)）。其中死時(shí)間修正、能譜濾波、峰漂修正、計(jì)算STP 等效高度均是對實(shí)測譜數(shù)據(jù)做的處理。擬合計(jì)算將凈計(jì)數(shù)率譜矩陣N 按矩陣表達(dá)式（3）進(jìn)行擬合即可得到K、U、Th、Rn 結(jié)果。結(jié)果輸出根據(jù)數(shù)據(jù)讀取的內(nèi)容，將坐標(biāo)信息（x，y）、雷達(dá)高度、K、U、Th、Rn 等重新整理匹配后輸出。

3 航空伽馬全能譜分析的編程

本次Python 編程主要引用了Csv（csv 數(shù)據(jù)格式）、NumPy（數(shù)組、矩陣計(jì)算）、Pandas（數(shù)據(jù)整理）、SciPy（工程科學(xué)計(jì)算庫）等第三方類庫。

3.1 數(shù)據(jù)讀取及整理

航空伽馬全能譜分析方法計(jì)算所用數(shù)據(jù)分為兩類，一類為實(shí)測譜數(shù)據(jù)，一般包含了測線、基點(diǎn)號、測量日期、時(shí)間、測量活時(shí)間、溫度、雷達(dá)高度、氣壓高度、K、U、Th、cos 窗數(shù)據(jù)以及0～255 道下測譜數(shù)據(jù)等；其次為參數(shù)文件，包括本底譜數(shù)據(jù)和K、U、Th、Rn 單位譜數(shù)據(jù)等，本次以“csv”數(shù)據(jù)格式為例。宇宙射線系數(shù)與飛機(jī)儀器本底參數(shù)合并為一個(gè)文件，包含道址（Char）、宇宙射線系數(shù)（cos_xs）、飛機(jī)儀器本底（air_bd）；單位譜數(shù)據(jù)包括道址（Char）、K 單位譜（K_COUNT）和衰減系數(shù)（Kb）、U 單位譜（U_COUNT）和衰減系數(shù)（Ub）、Th 單位譜（TH_COUNT）和衰減系數(shù)（THb）；大氣氡單位譜文件包括道址（Char）、大氣氡單位譜（Rn_unit）和衰減系數(shù)（Rnb）。

3.2 數(shù)據(jù)預(yù)處理

預(yù)處理包括死時(shí)間修正、能譜濾波、峰漂修正、計(jì)算STP 等效高度、宇宙射線本底計(jì)算和本底扣除以及STP 等效高度下單位譜計(jì)算等內(nèi)容。

1）死時(shí)間修正：利用每道計(jì)數(shù)率與活時(shí)間對每道進(jìn)行修正。

2）能譜濾波：由航空伽馬全能譜分析方法的原理可知，該方法是基于全譜進(jìn)行數(shù)據(jù)處理，那么譜線的噪聲水平如何，對于最終的計(jì)算結(jié)果有很大的影響。國內(nèi)外學(xué)者對能譜濾波已有較多的研究，如多項(xiàng)式擬合法、最小二乘法、小波濾波、Kalman（卡爾曼濾波）及NASVD（奇異值分解）方法等。其中小波濾波、Kalman（卡爾曼濾波方法）及NASVD 方法是目前較為流行的能譜降噪技術(shù)［16-22］，本次譜線濾波選擇了NASVD 方法，基于前人的研究結(jié)果對原始譜線性濾波處理，通過程序編程完成了濾波處理，核心部分編程如下：

其中，np.linalg.svd 為奇異值分解函數(shù)，np.sqrt（）為平方根計(jì)算函數(shù)，np.dot（）為矩陣的乘積函數(shù)，diag（）為創(chuàng)建對角矩陣函數(shù)，A_NA 為利用實(shí)測譜數(shù)據(jù)計(jì)算所得的單位方差矩陣，U_p 為奇異值分解產(chǎn)生的m×m 階酉矩陣，T_p 為奇異值分解產(chǎn)生的m×n 階對角矩陣，V_pT 是奇異值分解產(chǎn)生的n×n階酉矩陣；St_unit為實(shí)測譜數(shù)據(jù)歸一化行向量，V_p 為V_pT 的轉(zhuǎn)置矩陣。

3）峰漂修正：能譜峰漂一直是航空伽馬譜測量中最為關(guān)鍵的一項(xiàng)指標(biāo)，本次將能譜峰漂大于±0.5 道的譜均進(jìn)行修正。此次修正方法基于前人關(guān)于峰漂修正的理論編程完成［23-24］。

4）STP 等效高度計(jì)算：當(dāng)讀入一個(gè)實(shí)測譜數(shù)據(jù)后，利用溫度、雷達(dá)和氣壓數(shù)據(jù)根據(jù)公式（4），計(jì)算該測點(diǎn)對應(yīng)的STP 等效高度。

式中：Hstp—某一測點(diǎn)對應(yīng)的STP 等效高度，m；Hr—實(shí)測雷達(dá)高度，m；T是溫度，℃；Hb是氣壓高度，m。

5）計(jì)算宇宙射線本底：根據(jù)實(shí)測宇宙射線窗計(jì)數(shù)與宇宙射線本底修正系數(shù)計(jì)算宇宙射線本底譜。計(jì)算公式如下：

式中：Cos—某一測點(diǎn)對應(yīng)宇宙射線本底譜；cos—某一測點(diǎn)實(shí)測宇宙射線窗計(jì)數(shù)率，cos_xs—經(jīng)過高高度飛行得到的宇宙射線修正系數(shù)譜。

6）本底扣除：將實(shí)測譜數(shù)據(jù)減掉宙射線本底C和飛機(jī)儀器本底B，最后得到凈計(jì)數(shù)率譜。

7）STP 等效高度下單位譜計(jì)算：根據(jù)IAEA 推薦的公式（6），將地面1 m 處的單位譜（K、U、Th）根據(jù)其隨高度的變化系數(shù)換算至任一STP 等效高度上，具體公式如下：

式中：N?—經(jīng)換算修正至任一STP 等效高度后的單位譜；N0—地面上1 m 處的單位計(jì)數(shù)率譜；μ—K、U、Th、Rn 的單位譜隨高度變化的衰減系數(shù)譜；h—某一測點(diǎn)的STP 等效高度，m。

3.3 擬合計(jì)算

將扣除本底后的凈計(jì)數(shù)率譜數(shù)據(jù)與STP 高度下的單位譜進(jìn)行矩陣擬合計(jì)算，即可得到與K、U、Th、Rn 元素含量的結(jié)果矩陣，具體擬合函數(shù)定義如下：

其中，spe 為實(shí)測凈計(jì)數(shù)率的譜數(shù)據(jù)（只有下測0-255 道）；M 為某測點(diǎn)STP 等效高度上的單位譜矩陣；np.dot（）為求取兩個(gè)矩陣的乘積函數(shù)；.T 為求取某矩陣的轉(zhuǎn)置矩陣；np.linalg.inv（）為求取某矩陣的逆矩陣；rs為結(jié)果矩陣；包含K、U、Th、Rn 4個(gè)元素的含量。

最后從讀取的數(shù)據(jù)中整理提取x坐標(biāo)、y坐標(biāo)、雷達(dá)高度、窗數(shù)據(jù)等，從結(jié)果矩陣中提取K、U、Th、Rn 元素含量形成新數(shù)組，寫入到結(jié)果文件中。

4 應(yīng)用效果分析

根據(jù)程序設(shè)計(jì)進(jìn)行了軟件編程，其中核心代碼不到100 行，即完成了數(shù)據(jù)預(yù)處理和擬合計(jì)算兩大部分。本次選擇云南昆明-昆陽磷礦試驗(yàn)區(qū)進(jìn)行了數(shù)據(jù)處理。

4.1 程序運(yùn)行效果分析

試驗(yàn)區(qū)長、寬約30 km，線距為1 km，共計(jì)2.3萬余組能譜數(shù)據(jù)（256 道）。在i5CPU、8GBRAM 配置的Win 8系統(tǒng)中，該程序耗時(shí)約1 min 30 s左右；在相同配置的Win7、Win10 系統(tǒng)中，數(shù)據(jù)處理耗時(shí)誤差均不超過15 s，并在各系統(tǒng)上運(yùn)行穩(wěn)定，無退出、死機(jī)等情況。相比以往通過調(diào)用MATLAB、SCILAB 等數(shù)值計(jì)算軟件來擬合計(jì)算的編程，該程序運(yùn)行效率有大幅提升，且兼容性較好。

該程序在不同版本的windonws 系統(tǒng)上安裝簡單、便捷，無需單獨(dú)安裝不同版本的Microsoft.NET Framework 模 塊 以 及MATLAB 或SCILAB 等大型軟件，僅安裝開源的Anaconda 編程開發(fā)工具即可運(yùn)行程序，具有較好的移植性。

4.2 數(shù)據(jù)處理效果分析

將本程序計(jì)算的全譜法鉀、鈾、釷含量與傳統(tǒng)三窗法計(jì)算的鉀、鈾、釷含量對比分析，結(jié)果如見表1 和圖2?？傮w上，兩種方法的鉀、鈾、釷分布特征及變化規(guī)律基本一致，鉀元素分布特征最為接近，釷次之，鈾元素變化略大。如在縣街鎮(zhèn)西部以及海口鎮(zhèn)以西一帶，三窗法航放鈾含量偏高暈、高值暈規(guī)模較大、連續(xù)性較好，全譜法航放鈾含量連續(xù)性較差；其次在安寧市南部7 km 處，三窗法航放鈾含量呈南北向帶狀分布的偏高暈，而全譜法鈾含量則為偏低值特征。結(jié)合飛行高度數(shù)據(jù)分析，三窗法的航放鈾含量偏高是由航測飛行超高引起。

表1 試驗(yàn)區(qū)兩種數(shù)據(jù)處理方法結(jié)果統(tǒng)計(jì)Table 1 Statistics of processed result of test area by two processing methods

圖2 試驗(yàn)區(qū)兩種數(shù)據(jù)處理結(jié)果對比Fig.2 Result comparison of two processing methods in the test area

為進(jìn)一步驗(yàn)證本程序數(shù)據(jù)處理結(jié)果的可靠性以及與傳統(tǒng)三窗法的差別，選取了L2033 線中南段約14 km 測線和L1891 測線數(shù)據(jù)進(jìn)行了對比分析。兩測線全譜法與三窗法計(jì)算的鉀、鈾、釷含量結(jié)果如表2所示，圖3、圖4分別為L2033測線和L1891 測線鉀、鈾、釷含量結(jié)果對比圖。

L2033 線穿過昆陽磷礦區(qū)（164），鉀和釷含量變化趨勢一致（圖3）；全譜法鈾含量整體略大于三窗法鈾含量，并在最小值上體現(xiàn)較為明顯，但均值相差較小，相對誤差僅為6.1%（表2）。

圖3 L2033 線三窗法（TWM）與全譜法（FSA）計(jì)算結(jié)果對比Fig.3 Results comparison of L2033 Calculated by TWM and FSA

L1891 線穿過了觀音山磷礦區(qū)（114），兩種方法計(jì)算的結(jié)果基本一致（圖4），僅在171-203測點(diǎn)（灰線區(qū)間）內(nèi)，鉀、鈾含量結(jié)果變化較大，鉀含量最大值相差可達(dá)44.4%，鈾含量最大值相差為25.9%，且全譜法計(jì)算結(jié)果均低于三窗法（表2）。經(jīng)分析該區(qū)段為陡峭溝谷，飛行高度均大于140 m，局部高達(dá)425.8 m，下探測器接收到的計(jì)數(shù)極少，三窗法剝離及修正處理中產(chǎn)生了局部“假值”，而全譜法計(jì)算結(jié)果低于三窗法計(jì)算結(jié)果，綜合分析更符合該區(qū)地質(zhì)情況。

表2 L2033 與L1891 測線兩種方法計(jì)算結(jié)果統(tǒng)計(jì)Table 2 Result of measuring line L2033 and L1891 calculated by the two methods

圖4 L1891 線三窗法（TWM）與全譜法（FSA）計(jì)算結(jié)果對比Fig.4 Results comparison of L1891 calculated by TWM and FSA

5 結(jié) 論

1）利用Python 編程，用較少的代碼實(shí)現(xiàn)了航空伽馬全譜法數(shù)據(jù)處理和K、U、Th、Rn 等核素含量的計(jì)算；相比于傳統(tǒng)全譜法處理中調(diào)用大量數(shù)值計(jì)算軟件，Python 編程不但便捷、高效，且程序具有更好的兼容性和移植性。

2）云南昆明-昆陽磷礦試驗(yàn)區(qū)航放計(jì)算結(jié)果表明，基于Python 編程實(shí)現(xiàn)航空伽馬全能譜分析方法可行、結(jié)果可靠，處理結(jié)果局部特征更為突出，能有效地修正飛行超高引起的異常信息。

3）鑒于全譜法數(shù)據(jù)處理無需上測探測器數(shù)據(jù)進(jìn)行修正，即可進(jìn)行航空放射性數(shù)據(jù)處理；Python 編程在小型化航放測量數(shù)據(jù)處理以及其他全譜數(shù)據(jù)處理研究中值得推廣。