[发明专利]一种用于三维探地雷达的数据提取和展示方法

权利要求书

1.一种用于三维探地雷达的数据提取和展示方法：

步骤一，对测线漂移坐标进行过滤；

步骤二，对测线进行坐标插值及光顺平滑处理；

步骤三，对测线按照指定的位置进行雷达数据提取；

步骤四，对提取的雷达数据进行三维体绘制；

步骤一：对测线漂移坐标进行过滤，分为三个阶段，

阶段1：计算沿测线采集的两个相邻坐标点之间的直线距离D₁

坐标点P₁和P₂是通过实时动态差分设备(RTK)获得的，采用大地坐标系统，其坐标原点在地心，采用WGS-84参考椭球为基准面，由大地经度L、大地纬度B和大地高H表示，即WGS-84经纬度坐标；计算P₁和P₂两点之间的距离，需要把大地坐标系转换为高斯投影平面坐标系，通过高斯投影正算获得；

(1)高斯投影正算：

已知某点的(L,B)，求该点的(x,y)，即(L,B)→(x,y)的坐标变换；

(2)投影变换必须满足的条件：

中央子午线投影后为直线；

中央子午线投影后长度不变；

投影具有正形性质，即正形投影条件；

(3)投影过程：

在椭球面上有对称于中央子午线的两点P₁和P₂,它们的大地坐标分别为(L₁,B₁)或(l₁,B₁)及(L₂,B₂)或(l₂,B₂)式中l为椭球面上点的经度与中央子午线的经度差：l＝L-L₀,点在中央子午线之东,l为正，在西则为负，则投影后的平面坐标一定为P₁(x₁,y₁)和P₂(x₂,y₂)

(4)计算公式

转换精度精确至0.00lm时，用下式计算：

式中：

阶段2：计算沿测线采集的两个相邻测道之间的直线距离D₂

测道(trace)是通过硬件编码器设备(又称测距轮)获得的，对应阶段1中的坐标点P₁和P₂的测道分别为Tr₁和Tr₂；测距轮的精度Δd＝C/Nd，其中C为测距轮的周长，Nd为测距轮旋转一周的脉冲个数；探地雷达和RTK之间的数据采集通过测距轮同步触发；

经过标定的测距轮输出的测道信息与距离存在如下关系：

Distance＝Trace*DISTANCE_INTERVAL，其中DISTANCE_INTERVAL为道间距，Trace表示测道数；

所以：D₂＝(Tr_2–Tr₁)*DISTANCE_INTERVAL

阶段3：比较D₁和D₂

由于探地雷达和GPS之间的数据采集是通过测距轮同步触发设定过滤阈值threshold为5％，则过滤条件如下：

FABS(D₁-D₂)/D₂threshold，threshold＝0.05，FABS为浮点数求解绝对值；

顺序遍历测线上所有相邻点及对应的测道，计算D₁和D₂，根据过滤条件对发生漂移的坐标予以删除；由于在测线上通过RTK设备采集的相邻坐标点之间的距离通常在1米以内，且测线基于车载测距轮设备获取测道信息，视D₁、D₂之间的运动轨迹为直线，这是本算法成立的前提；

步骤二：对测线进行坐标插值及光顺平滑处理，对每条测线的位置文件内的所有坐标点进行三次贝塞尔曲线插值，使得每个测道都有对应的平面坐标，分为三个阶段：

阶段1：根据步骤一所述相同的计算方法，将每条测线的位置文件内的所有的大地坐标系经纬度坐标，通过高斯投影正算获得高斯投影平面坐标

阶段2：准备三次贝塞尔曲线函数

P(t)＝(1-t)3P1+3t(1-t)2P2+3t2(1-t)P3+t3P4

其中，P(t)为四个权重控制点的权重组合；

阶段3：对每条测线的投影平面坐标进行插值补齐

X(t)＝(+..i＝0,3)(Xi*C3i*pow(1-t,3-i)*pow(t,i)),t属于[0,1]

Y(t)＝(+..i＝0,3)(Yi*C3i*pow(1-t,3-i)*pow(t,i)),t属于[0,1]

其中，(+..i＝0,3)：累加符号，表示公式后面括号内将i分别用0,1,2,3代替，每个代替运算的结果全部加起来；Xi，Yi表示四个相邻的投影平面坐标；

C3i是二项式系数，其值为3！/(i！*(3-i)！)；

pow是求指数函数；

通过插值补齐使得每个测道都有对应的平面坐标；

步骤三：对测线按照指定的位置进行雷达数据提取，分为两个阶段：

阶段1：指定要提取的测线区域并查找测道号

(1)根据步骤一所述相同的计算方法，获得每条测线指定起始的P₁和P₂两个点的平面投影坐标(X₁,Y₁)、(X₂,Y₂)

(2)从步骤二插值补齐的坐标集合中，查找P₁和P₂的平面投影坐标(X₁,Y₁)、(X₂,Y₂)，得到对应的测道号(Tr₁)、(Tr₂)

测线1,P₁(X₁,Y₁),P₂(X₂,Y₂)→P₁(Tr₁),P₂(Tr₂)

测线2,P₁(X₁,Y₁),P₂(X₂,Y₂)→P₁(Tr₁),P₂(Tr₂)

......

测线N,P₁(X₁,Y₁),P₂(X₂,Y₂)→P₁(Tr₁),P₂(Tr₂)

阶段2：根据起始测道号提取每条测线指定区域的雷达振幅强度数据

1、雷达测线数据文件头总长度为3600字节，分两部分提取：

(1)文件头第一部分

长度为3200bytes，EBCDIC字符集，转换为ASCII码；

(2)文件头第二部分

长度为400bytes，32位、16位的整型，记录数据体信息；

2、雷达测线数据体由多个数据道组成，每道数据包含道头和采样数据：

(1)数据道的道头

长度为240bytes，32位、16位的整型，提取采样点数、采样间隔；

(2)采样数据

长度为采样数*采样点字节数，根据400字节头文件标示的数据体信息及阶段1获得的测道号，定位数据文件的位置，提取4字节IEEE浮点数，即雷达振幅强度数据；

数据体长度＝4*采样数*通道数*(Tr₂-Tr₁)

其中，三维阵列雷达内部通道数为16，采样数为256；

步骤四：对提取的雷达数据进行三维体绘制，分为两个阶段：

阶段1：为异常体准备原始数据

三维探地雷达采集的数据经过以上三个步骤的处理，得到异常体的原始数据，它们具有以X、Y和Z顺序排列的每个体素强度的8位强度值；

在OpenGL或DirectX中，为异常体制作的3D纹理中加载所有这些数据，但由于WebGL目前不支持存储或采用3D纹理，必须以2D纹理使用的方式存储它，由于这个原因，我们把强度值通过归一化生成并存储一个PNG图像文件，其中所有Z切片彼此相邻，形成2D切片的马赛克；我们利用开发的转换器工具，将提取的三维雷达数据转换为一个PNG图像镶嵌图，用于编码Alpha通道中每个体素的强度；一旦PNG文件作为2D纹理加载到内存中，我们就使用自定义的sampleAs3DTexture函数对其进行采样并生成形体渲染；

阶段2：对异常体进行体绘制

采用光线投射法将投射线从给定的角度上穿过异常体，对异常体内的像素信息作综合显示；结合深度、荫罩显示技术、旋转技术及信号强度切割技术，赋值影像以不同的伪彩与透明度，得到近似真实的三维结构的感受；

光线投射法基于图像序列的直接体进行绘制：从图像的每一个像素，沿固定方向发射一条光线，光线穿越整个图像序列，并在这个过程中，对图像序列进行采样获取颜色信息，同时依据光线吸收模型将颜色值进行累加，直至光线穿越整个图像序列，最后得到的颜色值就是渲染图像的颜色；

光线投射算法流程：

(1)正向面深度图渲染(Frontfacegeneration)

(2)背向面深度图渲染(Backfacegeneration)

(3)计算顶点位置、射线方向

(4)在片段着色程序中检索frontface，得到近距离(5)在片段着色程序中检索backface，得到远距离(6)计算当前射线方向上的最大穿越距离

(7)在片段着色程序中进行颜色合成和透明度累加(8)输出颜色值。