适用于全地形的航空γ能谱测量系统及其测量方法与流程

适用于全地形的航空
γ
能谱测量系统及其测量方法
技术领域
1.本发明属于航空γ能谱测量技术领域，具体为一种适用于全地形的航空γ能谱测量系统及其测量方法。


背景技术：

2.现行应用的航空γ能谱测量系统主要由配置大体积nai(tl)晶体探测器的航空γ能谱仪、导航定位系统、雷达高度计和数据收录仪等构成。航空γ能谱测量过程在记录每一测点γ能谱各道计数率数据的同时，还记录该测点的位置坐标和飞机离地高度信息。
3.在现行的常规测量中为便于能谱数据高度修正的程序自动化，采用负指数函数来模拟γ照射量率(计数率、计数)随测量高度的变化规律，采用如下公式进行高度衰减修正：
4.ns＝nm·
exp[μ
·
(h
0-h)]
[0005]
式中，ns是换算到规定测量高度h0上的计数率；
[0006]
nm是stp等等效高度上经剥离和本底修正的计数率；
[0007]
h为测量系统雷达高度计测量的离地高度值；
[0008]
μ为衰减系数。
[0009]
航空γ能谱各测点的能谱数据是来自该测点所对应整个测量周期作用带的γ辐射，而雷达高度计仅记录飞机所在测点位置的离地高度。因此进行高度衰减改正的前提是地面相对于待改正点为近似平坦且无限大的辐射体，即只适用于平原、戈壁和丘陵等地形比较平缓的测区。当在地形起伏变化的山区进行测量时，只通过高度改正必然产生误差，会导致凸地形上空改正不足，测量结果较真实值偏低；凹地形上空改正过大，测量结果偏高，往往造成异常的掩盖或者出现假异常，给后期的推断解释带来困难。因此在地形切割严重的测区只通过高度改正不能消除地形因素对航空γ能谱测量的影响。
[0010]
另一种航空γ能谱测量针对地形对航空γ能谱测量的影响，数据处理过程中地形被分割成若干个辐射体，根据辐射体γ场的公式计算探测点处所有辐射体的总γ场与半无限宽辐射体在该探测点同样测量条件下的比值，通过比值对该测点进行地形校正：其中方法一是高度改正过程中除了取测点本身的高度外,还沿测线方向前后各取若干点的高度值,并赋予它们由作用带所决定的不同系数，使高度改正考虑地形因素的影响获得的测量结果更准确；方法二是根据矩形辐射体航空γ辐射场理论，可以把空中任一测点的照射量率认为是地面上n个带状辐射体的累积贡献。如果这n个带状辐射体按照改正点对地面辐射体张角来划分，且改正点相对每个辐射体的张角皆为π/n，分别计算每个辐射体对改正点的影响，即按角度地形改正方法；方法三是作用带中垂直测线的方向将其视为处于不同地形高度的多个线状辐射源，线状辐射源的照射量率进行积分可得作用带在探测器处的照射量。
[0011]
上述针对地形较复杂区域的航空γ能谱测量，仍采用的是传统测量系统，对地形因素影响修正时，采用的是二维或准二维地形改正方法。当测区为复杂的地形时，二维地形改正有许多局限性，这种改正方法只考虑沿测线剖面方向的地形变化影响，未考虑垂直测
线方面以及与测线斜交方向的地形影响，而航空γ能谱仪探测器接收的信号来自整个360
°
方位作用带，故仅通过二维校正不能完全消除地形因素的影响。
[0012]
同时上述方法中沿测线剖面方向的地形数据一般采用系统测量过程中采集的gps高度值与雷达高度作差获得，而gps高度值与雷达高度的采集频率有限，因此获得的沿测线剖面方向的地形数据的精度也受到很大限制。
[0013]
随着gis和数字高程模型(dem)的普及，dem应用于航空γ能谱测量地形改正方法中，通过测区内的dem数据快速划分网格或采用有限元等方法进行三维地形影响校正。地形影响改正所采用的dem数据不是航空γ能谱测量过程中的实时地形数据，由于测量时间、测量环境等的差别，实时性较差，可能具有误差。同时通常获取的大图幅的dem数据是由不同扫描日期和不同数据源的数据通过拼接及滤波处理形成，该类型dem数据并不一定适合作为航空γ能谱单点地形改正所需较小范围作用带内的地形数据。


技术实现要素：

[0014]
本发明涉及适用于全地形测区的航空γ能谱测量系统，可实现对航空γ能谱测量数据三维地形影响改正，消除地形因素对测量结果的干扰，避免地形因素给航空γ能谱测量结果造成异常掩盖或出现假异常的影响，提高航空γ能谱测量的质量和精度，保证航空γ能谱测量技术在在地形平缓测区和地形复杂测区均可有效应用。
[0015]
具体的技术方案为：
[0016]
适用于全地形的航空γ能谱测量系统包括航空γ能谱仪、机载gps系统、imu、激光雷达扫描仪、数据采集存储计算机和飞机平台。
[0017]
机载gps系统分别向航空γ能谱仪、imu、激光雷达扫描仪提供pps信号和gps信息，以gps接收机输出的pps秒脉冲信号作为时间同步基准，对航空γ能谱仪、imu、激光雷达扫描仪三个子系统采样触发信号进行时间标定，实现各子系统输出结果的时间同步，通过gps输出信息中的utc时间，递推出不同测量频率下各子系统测量结果所对应的精确utc时间，为各子系统测量数据添加相应的时间戳，将不同测量子系统的结果统一到同一个时间轴上，实现在同一时间维度下数据精确采集。航空γ能谱仪、imu、激光雷达扫描仪将标有时间戳的测量数据包传输到数据采集存储计算机上。
[0018]
航空γ能谱测量单个测点的测量周期为1秒，机载gps系统输出的pps秒脉冲信号中的上升沿作为航空γ能谱仪上一测量周期结束控制信号，并作为下一测量周期的开始控制信号。航空测量为动态测量过程，在单个测量周期中航空γ能谱仪探测器位置是动态变化，将单个测量周期中航空γ能谱仪探测器移动范围的中间点位置作为该测量周期能谱数据对应的测点位置，航空γ能谱仪接收的机载gps系统提供更新率为2hz的gps信息，选用每秒的第2个gps信息中的坐标位置信息为该周期测量结果对应的测点坐标。
[0019]
将航空γ能谱测量作用带中所有激光雷达扫描仪获取的激光脚点视为点状放射源，将作用带表面细化为离散点源，航空γ能谱各测量周期的结果为作用带内多个激光脚点对应点源共同作用的结果。
[0020]
适用于全地形的航空γ能谱测量方法，包括以下步骤：
[0021]
s1系统获取数据处理
[0022]
s1.1pos数据处理
[0023]
(1)gps数据和imu数据的融合，经解算获得与imu相同采样频率的200hz gps天线坐标和imu姿态数据。
[0024]
(2)以激光雷达扫描仪数据采样频率为基准，对(1)中解算获得200hz的gps天线wgs84坐标数据进行内插，得到与激光点采集频率一致的20khz gps天线wgs-84坐标数据。
[0025]
(3)以激光雷达扫描仪数据采样频率为基准，对imu的200hz的姿态数据进行内插，得到与激光点采集频率一致的20khz imu平台姿态数据。
[0026]
s1.2激光扫描数据处理
[0027]
(1)通过激光雷达扫描仪获得的原始数据建立每个激光脚点在激光发射瞬时坐标系下坐标(0,0,ρ)。
[0028]
(2)通过激光雷达扫描仪瞬时扫描角θ，将每个激光脚点在激光发射瞬时坐标系下坐标(0,0,ρ)转化为激光脚点在扫描仪平台坐标系下坐标(x
ua
,y
ua
,z
ua
)。
[0029]
(3)通过偏心改正将激光脚点在扫描仪平台坐标系下坐标(x
ua
,y
ua
,z
ua
)转化为激光脚点在以gps天线为原点的imu平台坐标系下坐标(x
ga
,y
ga
,z
ga
)。
[0030]
(4)通过imu姿态角，将激光脚点在以gps天线为原点的imu平台坐标系下坐标(x
ga
,y
ga
,z
ga
)转化为激光脚点在以gps天线为原点在当地水平参考坐标系下的坐标(x
gia
,y
gia
,z
gia
)。
[0031]
(5)通过gps天线的wg-s84坐标，将激光脚点在以gps天线为原点的当地水平参考坐标系下的坐标(x
gia
,y
gia
,z
gia
)转化为激光脚点在wgs-84坐标系的坐标(x
84
,y
84
,z
84
)。
[0032]
(6)激光脚点在wgs-84坐标系的坐标(x
84
,y
84
,z
84
)转化为激光脚点的平面坐标(xn,ye)。
[0033]
(7)激光脚点以平面坐标(xn,ye)和海拔高度zh的形式来表示其三维坐标，形成测区范围内的激光点云三维坐标，即第j个激光脚点的三维坐标位置(x
fpj
,y
fpj
,z
fpj
)采用平面坐标系中坐标(x
fpj
,y
fpj
)和海拔高度z
fpj
表示。
[0034]
s1.3能谱窗口数据地形及高度综合改正公式的建立
[0035]
(1)建立航空γ能谱仪探测器在以gps天线为原点的平台坐标系下的坐标
[0036]
通过偏心测量建立与imu数据更新频率一致的200hz航空γ能谱仪探测器以gps天线为原点在平台坐标系下的坐标(x
gd
,y
gd
,z
gd
)。
[0037]
(2)建立航空γ能谱仪探测器以gps天线为原点的当地水平参考坐标系的坐标
[0038]
通过s1.1中pos数据解算获得的200hz imu姿态角将200hz的航空γ能谱仪探测器以gps天线为原点在平台坐标系下坐标(x
gd
,y
gd
,z
gd
)转化为以gps天线为原点的当地水平参考坐标系的坐标(x
gid
,y
gid
,z
gid
)。
[0039]
(3)建立航空γ能谱仪探测器在wgs-84坐标系的坐标
[0040]
通过s1.1中pos数据解算获得的200hz gps天线wgs-84坐标将200hz的航空γ能谱仪探测器以gps天线为原点的当地水平参考坐标系的坐标(x
gid
,y
gid
,z
gid
)转化为航空γ能谱仪探测器在wgs-84坐标系的坐标(x
84d
,y
84d
,z
84d
)。
[0041]
(4)建立航空γ能谱仪探测器的平面坐标
[0042]
航空γ能谱仪探测器在wgs-84坐标系下的坐标(x
84d
,y
84d
,z
84d
)转化为航空γ能谱仪探测器的平面坐标(x
nd
,y
ed
)。
[0043]
(5)建立航空γ能谱仪探测器以平面坐标和海拔高度表示的三维坐标
[0044]
采用航空γ能谱仪探测器平面坐标系下坐标(x
nd
,y
ed
)和海拔高度z
hd
来表示其三维坐标(x
nd
,y
ed
,z
hd
)。
[0045]
(6)构建航空γ能谱仪探测器γ射线信号接收边界曲面
[0046]
在平坦地区假设航空γ能谱仪探测器静止于某一飞行高度h处，经验总结航空γ能谱探测作用带(地表探测范围)为位于探测器正下方半径为h的圆面积,即航空γ能谱测量作用带的直径约为飞行高度的2倍，因此航空γ能谱仪探测器收γ射线信号的视场角为90
°
。
[0047]
根据航空γ能谱仪接收γ射线信号的视场角为90
°
，构建基于三维坐标的航空γ能谱仪探测器接收射线信号的边界曲面，是为以探测器为顶点、顶角为90
°
的圆锥曲面，表达式为：
[0048][0049]
其中(x,y,z)为曲面上任意一点的三维坐标，(x
nd
,y
ed
,z
hd
)为航空γ能谱仪探测器的三维坐标，采用其平面坐标系下坐标(x
nd
,y
ed
)和海拔高度z
hd
表示。
[0050]
(7)实际地形下航空γ能谱仪探测器作用带范围内激光脚点的确定
[0051]
根据辐射场的可叠加原理，任何形状的辐射源都可视为若干点源的叠加。航空γ能谱测量过程中空中探测器处的粒子注射率可视为作用带内无数个点源叠加作用的结果。将航空γ能谱测量作用带中所有激光雷达扫描仪获取得的激光脚点视为点状放射源，即可将作用带表面细化为离散点源，航空γ能谱各测量周期的结果为作用带内多个激光脚点对应点源共同作用的结果，即探测器处的粒子注射率为作用带内多个点源叠加作用的结果。
[0052]
单能不带电粒子的辐射场中，γ照射量率与粒子注量率成正比。单能γ辐射点源处于均匀介质中，距点源r处的γ射线粒子注量率φ点的计算公式为：
[0053][0054]
其中m为放射性物质的质量；μ为介质对γ射线的衰减系数；k是γ常数[c
·
m2/(kg
·s·
g)]；f为γ粒子注量率与γ照射量率的换算系数。
[0055]
在实际测区地形条件下探测器接收射线信号的边界曲面与地表曲面的交线为航空γ能谱测量作用带的边界，为不规则图形。激光雷达扫描仪获得的地表激光脚点视为地面点状辐射源，第j个激光脚点(点状辐射源)的三维坐标(x
fpj
,y
fpj
,z
fpj
)采用平面坐标系中坐标(x
fpj
,y
fpj
)和海拔高度z
fpj
表示。当激光脚点位置满足式1-3，该激光脚点在作用带范围内，将其视为作用带内的点源，参与航空γ能谱测量地形改正的计算。
[0056][0057]
其中(x
nd
,y
ed
,z
hd
)为采用航空γ能谱仪探测器的平面坐标系下坐标(x
nd
,y
ed
)和探测器的海拔高度z
hd
表示的三维坐标。
[0058]
(8)理想平坦地形条件下航空γ能谱仪探测器作用带范围内激光脚点的确定
[0059]
在理想平坦测区航空γ能谱仪探测器静止于某一飞行高度h处，根据探测器接收γ射线信号的90
°
视场角计算，航空γ能谱探测作用带为位于探测器正下方半径为h的圆面积，即航空γ能谱测量作用带的直径约为飞行高度的2倍；
[0060]
系统中的激光雷达扫描仪获得的地表激光脚点视为地面点状辐射源，第j个激光
脚点(点状辐射源)的三维坐标(x
fpj
,y
fpj
,z
fpj
)采用平面坐标系中坐标(x
fpj
,y
fpj
)和海拔高度z
fpj
表示。
[0061]
航空γ能谱仪探测器的三维坐标(x
nd
,y
ed
,z
hd
)采用探测器的平面坐标系下坐标(x
nd
,y
ed
)和探测器的海拔高度z
hd
来表示。
[0062]
航空γ能谱仪探测器与获得的第j个激光脚点的平面坐标分别为(x
nd
,y
ed
)和(x
fpj
,y
fpj
)，取(x
nd-x
fpj
)2+(y
ed-y
fpj
)2为最小值的激光脚点作为探测器在正下方地表上的投影，设该激光脚点的平面坐标为(x
′
fp y
′
fp
)，海拔高度为z
′
fp
。则假设理想平坦地形条件下各激光脚点海拔相同均为z
′
fp
，各激光脚点的三维坐标为(x
fpj
,y
fpj
,z
′
fp
)。探测器的离地高度hi的计算公式为式1-4：
[0063]
h＝z
hd-z
′
fp
ꢀꢀꢀ
1-4
[0064]
当理想平坦地形条件下激光脚点位置满足式1-5，视该激光脚点为作用带内的点源，参与航空γ能谱测量地形改正的计算。
[0065][0066]
(9)实际地形条件下航空γ能谱仪探测器静止于空中时的粒子注量率计算
[0067]
实际地形条件下航空γ能谱仪探测器静止于空中时的粒子注量率计算公式如下：
[0068][0069]
φc为实际地形条件下探测器静止于空中时的粒子注量率；
[0070]
φ
c点j
为实际测量地形条件下γ能谱测量作用带内的第j个点源在探测器处产生的γ射线粒子注量率；
[0071]
k为实际测量地形条件下航空γ能谱测量作用带内点源数量；
[0072]
m为放射性物质的质量；
[0073]
μ为介质对γ射线的衰减系数；
[0074]
k是γ常数[c
·
m2/(kg
·s·
g)]；f为γ粒子注量率与γ照射量率的换算系数。
[0075]rcj
为实际测量地形条件下航空γ能谱测量作用带内第j个点源距探测器的距离，计算公式如下:
[0076][0077]
式中(x
nd
,y
ed
)为探测器的平面坐标系下坐标，z
hd
为探测器的海拔高度；(x
fpj
,y
fpj
)为作用带内第j个激光脚点(点源)的平面坐标，z
fpj
为作用带内第j个激光脚点(点源)的海拔高度。
[0078]
(10)理想平坦地形条件下航空γ能谱仪探测器静止于空中时的粒子注量率计算
[0079]
理想平坦地形条件下航空γ能谱仪探测器静止于空中时的粒子注量率φ
pi
计算公式如下：
[0080][0081]
φ
p
为理想平坦地形条件下探测器静止于空中时的粒子注量率；
[0082]
φ
p点j
为理想平坦地形条件下γ能谱测量作用带内的第j个点源在探测器处产生的γ射线粒子注量率；
[0083]
n为理想平坦地形条件下航空γ能谱测量作用带内点源数量；
[0084]
m为放射性物质的质量；
[0085]
μ为介质对γ射线的衰减系数；
[0086]
k是γ常数[c
·
m2/(kg
·s·
g)]；f为γ粒子注量率与γ照射量率的换算系数。
[0087]rpj
为理想平坦地形条件下航空γ能谱测量作用带内第j个点源距探测器的距离，计算公式如下:
[0088][0089]
1-9式中(x
nd
,y
ed
)为探测器在平面坐标系中坐标，探测器海拔高度为z
hd
，理想平坦地形条件下作用带内第j个激光脚点(视为点源)的平面坐标为(x
fpj
,y
fpj
)，海拔高度统一为z
′
fp
。
[0090]
(11)地形改正公式建立
[0091]
航空γ能谱测量为动态测量过程，在1秒的测量周期内探测器处于运动状态，沿测线方向进行连续数据采集，每个测点处的计数率为相应测量周期中计数累加的结果。根据pos数据解算的输出频率，将每个测量周期内航空γ能谱仪探测器移动路径视为由探测器在空中200个不同位置组成，位置对应pos数据解算的探测器坐标；每个测量周期测量结果视为探测器处在空中200个不同位置测量共同作用的结果，即每个测量周期的测量对象为探测器静止在空中200个不同位置时所对应的作用带；
[0092]
每个测量周期为1秒的航空γ能谱测点的各能量窗口计数率地形及高度综合改正公式：
[0093][0094]
其中n为地形及高度改正后的计数率；n0为地形及高度改正前的计数率。
[0095]
pos数据解算的输出频率为200hz，即每个测量周期内探测器移动路径视为由探测器在空中200个不同坐标位置组成。
[0096]
φ
pi
φ
ci
分别为理想平坦地形条件下和实际地形条件下航空γ能谱仪探测器位于每个测量周期内200个不同坐标位置中的第i个位置时的粒子注量率。
[0097]hi
为航空γ能谱仪探测器位于每个测量周期内200个不同坐标位置中的第i个位置时的离地高度，h0为能谱地形和高度改正的统一目标高度。
[0098]
μ为γ射线高度衰减系数。
[0099]
将式1-6、1-8代入式1-10，每个测点航空γ能谱测量计数率地形及高度综合改正公式：
[0100]
[0101]
式中ki和ni分别为实际测量地形条件下和理想平坦地形条件下航空γ能谱仪探测器位于每个测量周期内200个不同坐标位置中的第i个位置时对应作用带内点源数量；
[0102]rcij
和r
pij
分别为实际测量地形条件下和理想平坦地形条件下航空γ能谱仪探测器位于每个测量周期内200个不同坐标位置中的第i个位置时对应作用带内第j个点源距探测器的距离，计算公式分别为公式1-7和1-9。
[0103]
s2航空γ能谱数据处理
[0104]
(1)分别对航空γ能谱测量获取的总窗计数率(tc)、钾窗计数率(k)、铀窗计数率(u)、钍窗计数率(th)进行飞机本底及宇宙射线本底改正；
[0105]
(2)总窗计数率(tc)、钾窗计数率(k)、铀窗计数率(u)、钍窗计数率(th)大气氡气本底改正；
[0106]
(3)钾道计数率(k)、铀道计数率(u)、钍道计数率(th)进行康普顿散射改正。
[0107]
(4)采用s1.3中建立的地形及高度综合改正公式，分别对总窗计数率(tc)、钾窗计数率(k)、铀窗计数率(u)、钍窗计数率(th)进行地形及高度综合改正。
[0108]
(5)将地形及高度综合改正后的钾道计数率(k)、铀道计数率(u)、钍道计数率(th)换算为对应的核素含量。
[0109]
本发明可消除地形因素对航空γ能谱测量结果的影响，避免地形因素给航空γ能谱测量造成异常掩盖或出现假异常的影响，提高航空γ能谱测量的质量和精度，提高后期地质推断解释的准确性，保证航空γ能谱测量技术在地形平缓测区和地形复杂测区均可有效应用。
附图说明
[0110]
图1为本发明的系统结构示意图；
[0111]
图2为本发明的pps信号；
[0112]
图3(1)为本发明空中测量数据处理流程之一；
[0113]
图3(2)为本发明空中测量数据处理流程之二；
[0114]
图4为本发明实际地形下γ射线探测器作用带范围示意图；
[0115]
图5为本发明理想平坦地形条件下γ射线探测器作用带范围示意图；
[0116]
图6为本发明航空γ能谱测量为动态测量过程探测器作用带范围示意图。
具体实施方式
[0117]
结合附图说明本发明的具体技术方案。
[0118]
适用于全地形的航空γ能谱测量系统由航空γ能谱仪、机载gps系统、imu、激光雷达扫描仪、数据采集存储计算机和飞机平台构成(如图1)。
[0119]
机载gps系统由机载gps天线及机载gps接收机构成，分别向航空γ能谱仪、imu、激光雷达扫描仪提供秒脉冲(pps)信号和gps信息(utc时间、gps天线位置、速度数据)；其中机载gps系统向航空γ能谱仪提供2hz的gps信息，向imu提供20hz的gps信息，向激光雷达扫描仪提供20hz的gps信息；以gps接收机输出的pps秒脉冲信号作为时间同步基准，对航空γ能谱仪、imu、激光雷达扫描仪三个子系统采样触发信号进行时间标定，通过gps输出信息中的utc时间，递推出不同测量频率下各子系统测量结果所对应的精确utc时间，为各子系统测
量数据添加相应的时间戳，实现不同的测量子系统在同一时间维度下进行数据精确采集，将结果数据统一到同一个时间轴上。
[0120]
航空γ能谱仪包括γ射线探测器和能谱仪主机，将航空γ能谱数据和相应的时间戳打包形成带时间戳的航空γ能谱数据包传输到数据采集存储计算机，更新频率为1hz；imu将imu数据和相应的时间戳打包形成带时间戳的imu数据包传输到数据采集存储计算机，更新频率为200hz；激光雷达扫描仪将所测激光测距数据和相应的时间戳打包形成带时间戳的激光测距数据包传输到数据采集存储计算机，更新频率为20khz。
[0121]
在航空γ能谱仪中gps输出的秒脉冲信号(pps)的上升沿作为每个测点测量周期的开始/停止测量控制信号。如图2所示，pps秒脉冲信号中的上升沿即作为上一测量周期结束控制信号，也作为下一测量周期的开始控制信号。
[0122]
航空γ能谱测量结果以每个测点的位置坐标和所对应能谱数据的形式表示，航空γ能谱测量是γ射线信号分析累加过程，单个测点能谱的测量周期为1秒。而航空测量为动态测量过程，飞机及γ射线探测器处于移动状态，位置是动态变化的，因此将单个测量周期中航空γ能谱仪探测器的移动范围的中间点位置作为该测量周期能谱数据对应的测点位置。航空γ能谱仪接收的gps系统提供的gps信息(gpgga格式数据)为2hz更新率，因此选用每秒的第2个gps信息中的坐标位置信息为该周期测量结果对应的测点坐标。
[0123]
适用于全地形的航空γ能谱测量方法，包括以下步骤：
[0124]
s1系统获取数据处理
[0125]
适用于全地形的航空γ能谱测量系统获取的空中测量数据经过如图3(1)和图3(2)所示数据处理流程，获得最终测量数据。具体的：
[0126]
s1.1 pos数据处理
[0127]
(1)gps数据(gps天线在wgs84坐标系下坐标)和imu数据(三轴加速度、三周角速度)的融合，经解算获得与imu相同采样频率的200hz gps天线坐标和imu姿态数据。
[0128]
(2)以激光雷达扫描仪数据采样频率为基准，对(1)中解算获得200hz的gps天线wgs84坐标数据进行内插，得到与激光点采集频率一致的20khz gps天线wgs84坐标数据。
[0129]
(3)以激光雷达扫描仪数据采样频率为基准，对imu的200hz的姿态数据进行内插，得到与激光点采集频率一致的20khz imu平台姿态数据。
[0130]
s1.2激光扫描数据处理
[0131]
(1)通过激光雷达扫描仪获得的原始数据建立每个激光脚点在激光发射瞬时坐标系下坐标(0,0,ρ)。
[0132]
(2)通过激光雷达扫描仪瞬时扫描角θ，将每个激光脚点在激光发射瞬时坐标系下坐标(0,0,ρ)转化为激光脚点在扫描仪平台坐标系下坐标(x
ua
,y
ua
,z
ua
)。
[0133]
(3)通过偏心改正将激光脚点在扫描仪平台坐标系下坐标(x
ua
,y
ua
,z
ua
)转化为激光脚点在以gps天线为原点的imu平台坐标系下坐标(x
ga
,y
ga
,z
ga
)。
[0134]
(4)通过imu姿态角，将激光脚点在以gps天线为原点的imu平台坐标系下坐标(x
ga
,y
ga
,z
ga
)转化为激光脚点在以gps天线为原点在当地水平参考坐标系下的坐标(x
gia
,y
gia
,z
gia
)。
[0135]
(5)通过gps天线的wgs84坐标，将激光脚点在以gps天线为原点的当地水平参考坐标系下的坐标(x
gia
,y
gia
,z
gia
)转化为激光脚点在wgs-84坐标系的坐标(x
84
,y
84
,z
84
)。
[0136]
(6)激光脚点在wgs-84坐标系的坐标(x
84
,y
84
,z
84
)转化为激光脚点的平面坐标(xn,ye)。
[0137]
(7)激光脚点以平面坐标(xn,ye)和海拔高度zh的形式来表示其三维坐标，形成测区范围内的激光点云三维坐标，即第j个激光脚点的三维坐标位置(x
fpj
,y
fpj
,z
fpj
)采用平面坐标系中坐标(x
fpj
,y
fpj
)和海拔高度z
fpj
表示。
[0138]
s1.3能谱窗口数据地形及高度综合改正公式的建立
[0139]
(1)建立航空γ能谱仪探测器在以gps天线为原点的平台坐标系下的坐标
[0140]
通过偏心测量建立与imu数据更新频率一致的200hz航空γ能谱仪探测器以gps天线为原点在平台坐标系下的坐标(x
gd
,y
gd
,z
gd
)。
[0141]
(2)建立航空γ能谱仪探测器以gps天线为原点的当地水平参考坐标系的坐标
[0142]
通过s1.1中pos数据解算获得的200hz imu姿态角将200hz的航空γ能谱仪探测器以gps天线为原点在平台坐标系下坐(x
gd
,y
gd
,z
gd
)转化为以gps天线为原点的当地水平参考坐标系的坐标(x
gid
,y
gid
,z
gid
)。
[0143]
(3)建立航空γ能谱仪探测器在wgs-84坐标系的坐标
[0144]
通过s1.1中pos数据解算获得的200hz gps天线wgs-84坐标将200hz的航空γ能谱仪探测器为以gps天线为原点的当地水平参考坐标系的坐标(x
gid
,y
gid
,z
gid
)转化为探测器在wgs-84坐标系的坐标(x
84d
,y
84d
,z
84d
)。
[0145]
(4)建立航空γ能谱仪探测器的平面坐标
[0146]
航空γ能谱仪探测器在wgs-84坐标系的坐标(x
84d
,y
84d
,z
84d
)转化为探测器的平面坐标(x
nd
,y
ed
)。
[0147]
(5)建立航空γ能谱仪探测器以平面坐标和海拔高度的形式表示的三维坐标
[0148]
γ射线探测器的三维坐标(x
nd
,y
ed
,z
hd
)，采用航空γ能谱仪探测器的平面坐标系下坐标(x
nd
,y
ed
)和探测器的海拔高度z
hd
来表示。
[0149]
(6)构建航空γ能谱仪探测器γ射线信号接收边界曲面
[0150]
在平坦地区假设航空γ能谱仪探测器静止于某一飞行高度h处，经验总结航空γ能谱探测作用带(地表探测范围)为位于探测器正下方半径为h的圆面积,即航空γ能谱测量作用带的直径约为飞行高度的2倍，因此航空γ能谱探测器收γ射线信号的视场角为90
°
。
[0151]
根据航空γ能谱仪接收γ射线信号的视场角为90
°
，构建基于三维坐标的航空γ能谱探测器接收射线信号的边界曲面，是为以探测器为顶点、顶角为90
°
的圆锥曲面，表达式为：
[0152][0153]
其中(x,y,z)为曲面上任意一点的三维坐标，(x
nd
,y
ed
,z
hd
)为航空γ能谱仪探测器的三维坐标，采用其平面坐标系坐标(x
nd
,y
ed
)和海拔高度z
hd
表示。
[0154]
(7)实际地形下航空γ能谱仪探测器作用带范围内激光脚点的确定(保护点)
[0155]
根据辐射场的可叠加原理，任何形状的辐射源都可视为若干点源的叠加。航空γ能谱测量过程中空中探测器处的粒子注射率可视为作用带内无数个点源叠加作用的结果。将航空γ能谱测量作用带中所有激光雷达扫描仪获取得的激光脚点视为点状放射源，即可
将作用带表面细化为离散点源，航空γ能谱各测量周期的结果为作用带内多个激光脚点对应点源共同作用的结果，即探测器处的粒子注射率为作用带内多个点源叠加作用的结果。
[0156]
单能不带电粒子的辐射场中，γ照射量率与粒子注量率成正比。单能γ辐射点源处于均匀介质中，距点源r处的γ射线粒子注量率φ点的计算公式为：
[0157][0158]
其中m为放射性物质的质量；μ为介质对γ射线的衰减系数；k是γ常数[c
·
m2/(kg
·s·
g)]；f为γ粒子注量率与γ照射量率的换算系数。
[0159]
√在实际测区地形条件下探测器接收射线信号的边界曲面与地表曲面的交线为航空γ能谱测量作用带的边界，为不规则图形。作用带内的激光脚点视为点源，探测器处的粒子注射率为作用带内多个点源叠加作用的结果，如图4。
[0160]
激光雷达扫描仪获得的地表激光脚点视为地面点状辐射源，第j个激光脚点(点状辐射源)的三维坐标位置(x
fpj
,y
fpj
,z
fpj
),采用平面坐标系中坐标(x
fpj
,y
fpj
)和海拔高度z
fpj
表示。当激光脚点位置满足式1-3，该激光脚点在作用带范围内，将其视为作用带内的点源，参与航空γ能谱测量地形改正的计算。
[0161][0162]
其中(x
nd
,y
ed
,z
hd
)为采用探测器的平面坐标系下坐标(x
nd
,y
ed
)和探测器的海拔高度z
hd
表示的三维坐标。
[0163]
(8)理想平坦地形条件下构建航空γ能谱仪探测器作用带范围内的激光脚点的确定(保护点)
[0164]
在理想平坦测区航空γ能谱仪探测器静止于某一飞行高度h处，根据探测器接收γ射线信号的90
°
视场角计算，航空γ能谱探测作用带(地表探测范围)为位于探测器正下方半径为h的圆面积，即航空γ能谱测量作用带的直径约为飞行高度的2倍，如图5。
[0165]
系统中的激光雷达扫描仪获得的地表激光脚点视为地面点状辐射源，第j个激光脚点(点状辐射源)的三维坐标(x
fpj
,y
fpj
,z
fpj
)采用平面坐标系中坐标(x
fpj
,y
fpj
)和海拔高度z
fpj
表示。
[0166]
航空γ能谱仪探测器的三维坐标(x
nd
,y
ed
,z
hd
)采用探测器的平面坐标系下坐标(x
nd
,y
ed
)和探测器的海拔高度z
hd
来表示。
[0167]
航空γ能谱仪探测器与获得的第j个激光脚点的平面坐标分别为(x
nd
,y
ed
)和(x
fpj
,y
fpj
)，取(x
nd-x
fpj
)2+(y
ed-y
fpj
)2为最小值的激光脚点作为探测器在正下方地表上的投影，设该激光脚点的平面坐标为(x
′
fp y
′
fp
)，海拔高度为z
′
fp
。则假设理想平坦地形条件下各激光脚点海拔相同均为z
′
fp
，各激光脚点的三维坐标为(x
fpj
,y
fpj
,z
′
fp
)。探测器的离地高度hi的计算公式为式1-4：
[0168]
h＝z
hd-z
′
fp 1-4
[0169]
当理想平坦地形条件下激光脚点位置满足式1-5，视该激光脚点为作用带内的点源，参与航空γ能谱测量地形改正的计算。
[0170]
[0171]
(9)实际地形条件下航空γ能谱仪探测器静止于空中时的粒子注量率计算
[0172]
实际地形条件下航空γ能谱仪探测器静止于空中时的粒子注量率φ
ci
计算公式如下：
[0173][0174]
φc为实际地形条件下探测器静止于空中时的粒子注量率；
[0175]
φ
c点j
为实际测量地形条件下γ能谱测量作用带内的第j个点源在探测器处产生的γ射线粒子注量率；
[0176]
k为实际测量地形条件下航空γ能谱测量作用带内点源数量；
[0177]
m为放射性物质的质量；
[0178]
μ为介质对γ射线的衰减系数；
[0179]
k是γ常数[c
·
m2/(kg
·s·
g)]；f为γ粒子注量率与γ照射量率的换算系数。
[0180]rcj
为实际测量地形条件下航空γ能谱测量作用带内第j个点源距探测器的距离，计算公式如下:
[0181][0182]
式中(x
nd
,y
ed
)为探测器的平面坐标系下坐标，z
hd
为探测器的海拔高度；(x
fpj
,y
fpj
)为作用带内第j个激光脚点(点源)的平面坐标，z
fpj
为第j个激光脚点(点源)的海拔高度。
[0183]
(10)理想平坦地形条件下航空γ能谱仪探测器静止于空中时的粒子注量率计算
[0184]
理想平坦地形条件下航空γ能谱仪探测器静止于空中时的粒子注量率计算公式如下：
[0185][0186]
φ
p
为理想平坦地形条件下探测器静止于空中时的粒子注量率；
[0187]
φ
p点j
为理想平坦地形条件下γ能谱测量作用带内的第j个点源在探测器处产生的γ射线粒子注量率；
[0188]
n为理想平坦地形条件下航空γ能谱测量作用带内点源数量；
[0189]
m为放射性物质的质量；
[0190]
μ为介质对γ射线的衰减系数；
[0191]
k是γ常数[c
·
m2/(kg
·s·
g)]；f为γ粒子注量率与γ照射量率的换算系数。
[0192]rpj
为理想平坦地形条件下航空γ能谱测量作用带内第j个点源距探测器的距离，计算公式如下:
[0193][0194]
1-9式中(x
nd
,y
ed
)为探测器在平面坐标系中坐标，探测器海拔高度为z
hd
，理想平坦地形条件下作用带内第j个激光脚点(视为点源)的平面坐标为(x
fpj
,y
fpj
)，海拔高度统一为z
′
fp
。
[0195]
(11)地形改正公式建立
[0196]
航空γ能谱测量为动态测量过程，在1秒的测量周期内探测器处于运动状态，沿测线方向进行连续数据采集，每个测点处的计数率为相应测量周期中计数累加的结果。根据pos数据解算的输出频率(200hz)，将每个测量周期内航空γ能谱仪探测器移动路径视为由探测器在空中200个不同位置组成，位置对应pos数据解算的探测器坐标；每个测量周期测量结果视为探测器处在空中200个不同位置测量共同作用的结果，即每个测量周期的测量对象为探测器静止在空中200个不同位置时所对应的作用带(如图6)。
[0197]
每个测量周期为1秒的航空γ能谱测点的各能量窗口计数率地形及高度综合改正公式：
[0198][0199]
其中n为地形及高度改正后的计数率；n0为地形及高度改正前的计数率。
[0200]
pos数据解算的输出频率为200hz，即每个测量周期内探测器移动路径视为由探测器在空中200个不同坐标位置组成。
[0201]
φ
pi
φ
ci
分别为理想平坦地形条件下和实际地形条件下航空γ能谱仪探测器位于每个测量周期内200个不同坐标位置中的第i个位置时的粒子注量率。
[0202]hi
为航空γ能谱仪探测器位于每个测量周期内200个不同坐标位置中的第i个位置时的离地高度，h0为能谱地形和高度改正的统一高度。
[0203]
μ为γ射线高度衰减系数。
[0204]
将式1-6、1-8代入式1-10，每个测点航空γ能谱测量计数率地形及高度综合改正公式：
[0205][0206]
式中ki和ni分别为实际测量地形条件下和理想平坦地形条件下航空γ能谱仪探测器位于每个测量周期内200个不同坐标位置中的第i个位置时对应作用带内点源数量；
[0207]rcij
和r
pij
分别为实际测量地形条件下和理想平坦地形条件下航空γ能谱仪探测器位于每个测量周期内200个不同坐标位置中的第i个位置时对应作用带内第j个点源距探测器的距离，计算公式分别为公式1-7和1-9。
[0208]
s2航空γ能谱数据处理
[0209]
(1)分别对航空γ能谱测量获取的总窗计数率(tc)、钾窗计数率(k)、铀窗计数率(u)、钍窗计数率(th)进行飞机本底及宇宙射线本底改正；
[0210]
(2)总窗计数率(tc)、钾窗计数率(k)、铀窗计数率(u)、钍窗计数率(th)进行大气氡气本底改正；
[0211]
(3)钾道计数率(k)、铀道计数率(u)、钍道计数率(th)进行康普顿散射改正。
[0212]
(4)采用s1.3中建立的地形及高度综合改正公式，分别对总窗计数率(tc)、钾窗计数率(k)、铀窗计数率(u)、钍窗计数率(th)进行地形及高度综合改正。
[0213]
(5)将地形及高度综合改正后的钾道计数率(k)、铀道计数率(u)、钍道计数率(th)换算为对应的核素含量。