地球化学参数统计

Answer 1

（一）统计参数的理论基础

在进行地球化学统计时，一般采用2类统计参数：一类反映属性数据集中特性，包括频数、平均数、中数、众数；另一类反映数据离散程度，包括极差、方差、标准差、变异系数等。其相应的数学模型如下所述（张文彤，2002；余建英等，2003）。

（1）频数

将变量x_i（i＝1，2，…，n）按大小顺序排列，并按一定的间距分组。变量在各组出现或发生的次数称为频数，一般用ƒ_i 表示。各组频数与总频数之比叫做频率，按如下公式计算：

浙江省农业地质环境GIS设计与实现

根据大数定理，当n相当大时，频率可近似地表示事件的概率。

计算出各组变量的频率后，若以纵轴表示频率，横轴表示分组，就可做出频率直方图，用以表示事件发生的频率和分布状况。变量的频数分析可以使人们非常清楚地了解变量值的分布情况，而在进行数理统计时，通常要求数据呈正态分布。因此，在针对地球化学元素的基准值统计之前，进行频数分析可以展现区域内各个元素的数据分布状态，符合正态分布的可以进行后续统计分析，不符合正态分布的要剔除异常值，使之符合正态分布。

（2）平均数

平均数反映了数据取值的集中位置，常以

表示。对地球化学元素的基准值求平均数，可以从整体上大致掌握元素的含量水平。对于数据X_i（i＝1，2，…，n）通常有简单算术平均数和加权算术平均数。其中简单算术平均数的计算公式为

浙江省农业地质环境GIS设计与实现

加权算术平均数的计算公式为

浙江省农业地质环境GIS设计与实现

式中：P_i 为数据X_i 的权值。

（3）众数

众数是具有最大可能出现的数值。如果数据X是离散的，则称X中出现最大可能性的值x为众数；如果X是连续的，则以X分布的概率密度P（x）取最大值的x为X的众数。显然，众数可能不是唯一的。

地球化学元素基准值统计中针对某元素求得的众数表明该元素在此含量水平上的取值最多，也反映了区域内元素含量的基本水平。

（4）中数

对于有序数据集X，如果有一个数x，能同时满足以下2式：

浙江省农业地质环境GIS设计与实现

则称x为数据集X的中数，记为x_Me。

若X的总项数为奇数，则中数为

浙江省农业地质环境GIS设计与实现

若X的总项数为偶数，则中数为

浙江省农业地质环境GIS设计与实现

可见，频数（率）、平均数、众数、中数都反映了数据的总体取值水平，当数据分布状态最佳时，平均数、众数、中数在理论上应该是相等的，但是这种理想状态在实际中很少出现。在地球化学统计工作中，采用以上4个参数就足以反映区域内元素的集中含量水平。

（5）极差

极差是一组数据中最大值与最小值之差，即

R＝max｛x₁，x₂，…，x_n｝ -min｛x₁，x₂，…，x_n｝

（6）方差与标准差

方差是均方差的简称，是以离差平方和除以变量个数求得的，记为s²，即

浙江省农业地质环境GIS设计与实现

标准差是方差的平方根，记为

浙江省农业地质环境GIS设计与实现

极差、方差与标准差是反映数据离散程度的参数，这几个参数的统计结果表明区域内地球化学元素在不同采样点的含量水平差距，其值越大，说明样点间的含量水平差异越大；反之说明各个样点的元素含量较为均匀。

由于元素在土壤中的地球化学分布与分配受多种环境条件的限制，如地质作用、成矿作用、表生作用及人为活动等，因此寻找地球化学异常是较为重要的一项工作。从统计学角度来看，通常认为地球化学异常是在一定置信度下相对于背景的偏离。浙江省农业地质环境调查工作中，以元素背景值的平均数加2倍标准差作为元素含量的异常上限，以元素背景值的平均数减2倍标准差作为元素含量的异常下限，经过参数统计找到元素的地球化学异常，就可以从地质因素、人为影响等方面入手，进一步分析其导致异常的原因。

（7）变异系数

在衡量区域内地球化学元素的含量水平差距时，如果在元素间进行比较，经常会遇到不同元素含量相差比较大的情况，有时甚至相差几个数量级，这时使用方差或是标准差来衡量很不方便。因此采用另一个参数——变异系数。

变异系数用来衡量数据在时间和空间上的相对变化的程度，记为C：

浙江省农业地质环境GIS设计与实现

式中：s为标准差；

为平均数。

可以看出变异系数是无量纲的量，它消除了不同化学元素含量在数量级上的差异，反映的也是数据取值的离散程度，值越大，含量差异越大；值越小，含量取值越均匀。

（二）基于GIS的参数统计

地球化学参数统计常规方法是先根据获取的数据源，使用电子表格的计算功能进行半自动化统计。但在实际生产过程中，往往要求按专业分类、行政区划进行参数统计，如统计某县市范围内按土壤类型的变异系数时会变得很困难。不仅统计的速度慢、效率低，而且统计的结果容易出错。因此，采用基于GIS的参数统计十分适合复杂的生产要求。

图7-1 基于GIS的参数统计流程图

1.实现思路

要实现县（市）范围内按土壤类型的地球化学参数统计，首先要整理测试分析数据，通过测试分析数据的经纬度投影成区域测试数据点位图，要求与行政区划图层的格式、坐标系一致。实现属性表（测试分析表数据）向空间数据（区域地球化学分析测试点位图层）的转换。

为获取统计范围内的数据，可采用行政区划作为裁剪框与地球化学分析测试点位图层进行图形裁剪，并以土壤分布类型作为统计单元，再将裁剪后的结果点数据图层与土壤类型分布区图层再进行点与区相交分析，同时测试点位图层与土壤类型分布图层的属性也会相应地分类合并，在此基础上再选择需统计的字段（元素）进行地球化学的参数统计，形成地球化学参数统计报表。具体流程如图7-1所示。

2.关键代码

（1）点与区相交分析

Private Sub分析_点_区（）

Dim lyr1 As MapObjects2.maplayer，lyr2 As MapObjects2.maplayer，lyr3 As MapOb⁃jects2.maplayer

Setlyr1＝m_map.Layers（Comb_1）＇点图层

Setlyr2＝m_map.Layers（Comb_2）＇区图层

Dim rd1 As MapObjects2.RecordSet，rd2 As MapObjects2.RecordSet

Dim tagrd As MapObjects2.RecordSet

Dim shp1 As MapObjects2.point，newshp As MapObjects2.point

Dim i As Long

Dim nums As Long

Screen.MousePointer＝vbHourglass

If m_判别方法＝相减Then

CreateNewLyr lyr1，lyr2，点，False

Else

CreateNewLyr lyr1，lyr2，点，True

End If

Set tagrd＝m_tagerlyr.Records

tagrd.AutoFlush＝False

Set rd1＝lyr1.Records

nums＝ReturnNums（rd1）

i＝0

Do Until rd1.EOF

Setshp1＝rd1.Fields（＂shape＂）.Value

Select Case m_判别方法

Case判别

Add点区相交 tagrd，rd1，lyr2，shp1

Add点区相减 tagrd，rd1，lyr2，shp1

Case相交

Add点区相交 tagrd，rd1，lyr2，shp1

Case相减

Add点区相减 tagrd，rd1，lyr2，shp1

End Select

i＝i + 1

ProgressBar1.Value＝i/nums*100＃

ProgressBar1.Refresh

rd1.MoveNext

Loop

tagrd.AutoFlush＝True

ProgressBar1.Value＝0

ProgressBar1.Refresh

Screen.MousePointer＝vbDefault

End Sub

＇*********************点对区相交分析****************************

PrivateSubAdd点区相交（TagerRdAsMapObjects2.RecordSet，srdAsMapObjects2.RecordSet，lyr2 AsMapObjects2.maplayer，tshp AsMapObjects2.point）

Dim rd2 As MapObjects2.RecordSet

Dim flds As MapObjects2.Fields，fld As MapObjects2.Field

Dim tfldstr As String

Setrd2＝lyr2.SearchShape（tshp，moAreaIntersect，＂＂）

Dim t_Id As Long

If returnRdEof（rd2）＝ False Then

With TagerRd

.AddNew＇相应属性合并

t_Id＝t_Id + 1

.Fields（＂objectid＂）.Value＝t_Id

Set flds＝srd.Fields

For Each fld In GetSortFld（srd）

tfldstr＝fld.name

If Not .Fields（tfldstr）Is Nothing Then

.Fields（tfldstr）.Value＝fld.Value

End If

Next

Set flds＝rd2.Fields

For Each fld In GetSortFld（rd2）

tfldstr＝newfldname（fld.name，TagerRd.Fields）

If Not .Fields（tfldstr）Is Nothing Then

.Fields（tfldstr）.Value＝fld.Value

End If

Next

.Fields（＂shape＂）.Value＝tshp

.Update

End With

End If

End Sub

（2）参数统计

Public Function dataEDIT（t数据类别As String，t数据表单As String，t元素字段As String，trd As ADODB.RecordSet，mrd As MapObjects2.RecordSet，sField As String，Is去处离散数据As Boolean）As Boolean

dataEDIT＝False

Dim lstrSQL As String

Dim Dcount As Integer

Dim i As Long，j As Long

Dim TempData（）As Double

Dim dbname As String

Dim str As String

Dim FV（13）As String

Dim tt As Integer

Dim Avg，s，ss，cv，Dsum As Double

Dim k，n，m As Integer

On Error GoTo ErrCl

Dim TNums As Long

Dcount＝ReturnNums（mrd）

ReDim TempData（2，Dcount）

i＝0

Do Until mrd.EOF

If mrd.Fields（sField）.Value＞0 Then

TempData（0，i）＝ mrd.Fields（sField）.Value

i＝i + 1

End If

mrd.MoveNext

Loop

If i＜2 Then

Exit Function

End If

＇开始数据统计

dbname＝＂test＂

FV（0）＝＂深层数据＂

FV（2）＝＂＂

FV（3）＝ sField

FV（8）＝ TempData（0，Dcount 1）

FV（9）＝ TempData（0，0）

FV（12）＝＂＂

Dsum＝0＃

For i＝0 To Dcount 1

Dsum＝Dsum + TempData（0，i）

Next

If TempData（k，i）＜（Avg+2 s）And TempData（k，i）＞（Avg 2 s）Then

Avg＝Dsum/Dcount ＇**第一次平均值μ**

Dsum＝0＃

For i＝0 To Dcount 1

Dsum＝Dsum +（TempData（0，i）Avg）＾2

Next

s＝Sqr（Dsum/（Dcount 1））＇**第一次离差＇=**

k＝0

j＝1

m＝Dcount 1

Do While True

＇**去除离散数据**

If Is去处离散数据 Then

n＝0

For i＝0 To m

TempData（j，n）＝ TempData（k，i）

n＝n + 1

End If

Next

TNums＝n

Else

n＝0

For i＝0 To m

TempData（j，n）＝ TempData（k，i）

n＝n + 1

Next

TNums＝n

End If

m＝n - 1

＇**求样品数、平均值、离差

Dsum＝0＃

For i＝0 To m

Dsum＝Dsum + TempData（j，i）

Next

Avg＝Dsum /（m + 1）

Dsum＝0＃

For i＝0 To m

Dsum＝Dsum +（TempData（j，i）Avg）＾2

Next

If Dsum＝0 Then

GoTo allok

Else

ss＝Sqr（Dsum/m）

End If

If Abs（ss s）＜＝0.0001 Then GoTo allok

s＝ss

tt＝k

k＝j

j＝tt

Loop

allok：

m＝m + 1

FV（4）＝ m

FV（5）＝ Avg

FV（6）＝ ss

FV（7）＝ ss/Avg

If（m 2*（m/2））＝0Then＇求中位数

i＝m/2 1

FV（11）＝ TempData（j，i）

i＝i + 1

FV（11）＝（FV（11）+ TempData（j，i））/2＃

Else

i ＝（m + 1）/2 1

FV（11）＝ TempData（j，i）

End If

m＝m 1

Dimzzmax（2，5000）AsDouble＇求众数

k＝0

zzmax（0，k）＝ TempData（j，0）

zzmax（1，k）＝ 0

For i＝1 To m

If TempData（j，i）＝ zzmax（0，k）Then

zzmax（1，k）＝ zzmax（1，k）+ 1＃

Else

k＝k + 1

zzmax（0，k）＝ TempData（j，i）

End If

Next

tt＝zzmax（1，0）

s＝zzmax（0，0）

For i＝1 To k

If zzmax（1，i）＞tt Then

tt＝zzmax（1，i）

s＝zzmax（0，i）

End If

Next

If tt＞0 Then

FV（10）＝ s

Else

＇*****************************存入数据库******************************

trd！平均值＝Format（FV（5），＂＃.＃＃＃＂）

trd！统计样品数＝TNums

trd！标准离差＝Format（FV（6），＂＃.＃＃＃＂）

trd！变异系数＝Format（FV（7），＂＃.＃＃＃＂）

trd！极大值＝Format（FV（8），＂＃.＃＃＃＂）

trd！极小值＝Format（FV（9），＂＃.＃＃＃＂）

trd！众值＝Format（FV（10），＂＃.＃＃＃＂）

FV（10）＝＂＂

End If

Dim InsertSQL As String

trd.AddNew

trd！数据类别＝t数据类别

trd！数据表单＝t数据表单

trd！元素字段＝t元素字段

trd！原始样品数＝FV（4）

trd！中位数＝Format（FV（11），＂＃.＃＃＃＂）

dataEDIT＝True

trd.Update

ErrCl：

End Function

（三）实例

农业地质环境调查工作获取了大量的样品分析数据针对不同调查目的、不同样品，所测试的项目各有不同，但无论是何种测试数据，其分析统计方法是相同的。对于地球化学元素基准值，通常都是统计以上介绍的几个常用参数，用以分析其数据特征。下面以“萧山区农业地质环境调查”子项目中，表层土壤的环境地球化学背景值统计分析为例，展示AGEIS系统的地球化学参数统计功能。

选萧山区不同的土壤类型为统计单元，调取表层土壤单点样分析数据，共计227件样品，分析测试项目为8项元素全量，对其背景值进行剔除离差后的参数统计。统计结果见表7-1。