基于ArcGIS Engine的沉積物粒徑趨勢(shì)分析模型開(kāi)發(fā)與應(yīng)用

2016-12-27 09:08:15閆吉順方海超孫家文于永海

地理空間信息 2016年7期

關(guān)鍵詞：趨勢(shì)分析

閆吉順，方海超，孫家文*，于永海，李家

（1. 遼寧師范大學(xué) 城市與環(huán)境學(xué)院，遼寧大連 116029；2. 國(guó)家海洋環(huán)境監(jiān)測(cè)中心，遼寧大連 116023）

基于ArcGIS Engine的沉積物粒徑趨勢(shì)分析模型開(kāi)發(fā)與應(yīng)用

閆吉順1，方海超2，孫家文2*，于永海2，李家1

（1. 遼寧師范大學(xué) 城市與環(huán)境學(xué)院，遼寧大連 116029；2. 國(guó)家海洋環(huán)境監(jiān)測(cè)中心，遼寧大連 116023）

以Microsoft Visual Studio 2010為開(kāi)發(fā)平臺(tái)，嵌入ArcGIS Engine集成組件，以Gao-Collins方法為基礎(chǔ)，利用C#程序語(yǔ)言開(kāi)發(fā)沉積物粒度分析模型、半方差函數(shù)模型和沉積物粒徑趨勢(shì)分析模型等分析模型工具；并以大連市金石灘84個(gè)站位采樣數(shù)據(jù)為實(shí)例，利用該工具進(jìn)行了沉積物粒徑趨勢(shì)分析。

Gao-Collins；沉積物粒度；半方差函數(shù)；粒徑趨勢(shì)分析；ArcGIS Engine

在海洋沉積學(xué)中，通過(guò)沉積物粒度參數(shù)的空間變化（粒徑趨勢(shì)），反推沉積物凈輸運(yùn)方向的思想和相應(yīng)的一套計(jì)算方法稱為粒徑趨勢(shì)分析[1]。眾多研究者進(jìn)行了深入研究[2]，例如劉濤[1]等基于“動(dòng)力組分”思想建立了沉積物粒徑趨勢(shì)模型；Mclaren和Bowles把平均粒徑、分選系數(shù)和偏態(tài)3個(gè)粒度參數(shù)結(jié)合起來(lái)，提出了一維粒徑輸運(yùn)趨勢(shì)模型[3]；Gao和Collins提出了一種粒徑趨勢(shì)分析方法，即Gao-Collins方法[4]；馬菲[5]等采用地統(tǒng)計(jì)法對(duì)北部灣東部海域沉積物粒徑趨勢(shì)進(jìn)行分析；徐志偉[6]等提出利用多元統(tǒng)計(jì)及物源分析對(duì)北部灣東部海域沉積物粒徑趨勢(shì)進(jìn)行分析。

粒徑趨勢(shì)分析方法已經(jīng)日益成熟，但過(guò)程繁雜，計(jì)算規(guī)模較大，若實(shí)現(xiàn)模型計(jì)算自動(dòng)化與可視化將能提高研究工作的效率與準(zhǔn)確性。基于地理數(shù)據(jù)具有空間相關(guān)性的特征，本文選擇Microsoft Visual Studio 2008為開(kāi)發(fā)平臺(tái)，嵌入ArcGIS Engine集成組件，不僅可保持地理數(shù)據(jù)的空間相關(guān)性，而且能滿足繁瑣與較大規(guī)模的計(jì)算需求。

1 相關(guān)模型與理論

1.1 沉積物粒度分析模型

粒度分析的目的是獲得沉積物粒徑的概率分布曲線。習(xí)慣上將沉積物樣品分為粗顆粒和細(xì)顆粒（以0.063 mm粒徑為界），細(xì)顆粒物質(zhì)用移液管法分析，粗顆粒物質(zhì)用篩法分析，最后將兩部分綜合起來(lái)，獲得完整的粒徑分布曲線。要注意的是，沉積學(xué)中的粒徑通常不以mm表示，而是以DΦ來(lái)表示：

式中，D為以mm計(jì)的粒徑；DΦ為Φ變換后的無(wú)量綱粒徑。

對(duì)于一定的粒度分布曲線，可以計(jì)算出一系列粒度參數(shù)。最常用的粒度參數(shù)包括平均粒徑、分選系數(shù)和偏態(tài)系數(shù)[7]。

粒度參數(shù)的計(jì)算方法主要有圖解法和矩值法，包括Folk和Ward的圖解法公式、Collias等的矩法公式和McManus的矩法公式。對(duì)3種公式綜合對(duì)比后認(rèn)為Collias等的矩法公式反映沉積物粒度分布特征更加靈敏、可靠[8]，本文采用該公式。

式中，μ、δ、Sk和Ku分別為平均粒徑、分選系數(shù)、偏態(tài)和峰態(tài)；Xi為粒級(jí)組中值（Ф）；fi為各粒級(jí)范圍的百分含量。

1.2 半方差函數(shù)模型

對(duì)表層沉積物采樣間距的分析方法，使用了地學(xué)統(tǒng)計(jì)法中的半方差分析，即對(duì)于平均粒徑、分選系數(shù)和偏態(tài)等參數(shù)，構(gòu)造以下函數(shù)[9]：

式中，h為采樣間距；Nh為樣品的個(gè)數(shù)；Zi為原點(diǎn)處的參數(shù)值；Zi+h為h之外站位的參數(shù)值；γh為半方差值。以半方差γh為縱坐標(biāo)、采樣間隔h為橫坐標(biāo)作圖得到半方差圖。

1.3 Gao-Collins方法沉積物粒徑趨勢(shì)分析模型[10]

粒徑趨勢(shì)指沉積物粒度參數(shù)平面分布的變化趨勢(shì)。對(duì)于一個(gè)海區(qū)，可以布設(shè)一定的采樣網(wǎng)格進(jìn)行底質(zhì)取樣。在采樣網(wǎng)格中，如果考慮任意兩個(gè)相鄰的采樣點(diǎn)A和B，則這兩個(gè)采樣點(diǎn)的粒度參數(shù)之間有多種可能的空間變化，如從采樣點(diǎn)A至采樣點(diǎn)B，分選系數(shù)可能減小，平均粒徑可能變大，這些都代表粒徑趨勢(shì)的不同類型。用多個(gè)粒度參數(shù)可以形成組合的粒徑趨勢(shì)，一般而言，用n個(gè)粒度參數(shù)可構(gòu)成2n種粒徑趨勢(shì)，如用平均粒徑、分選系數(shù)和偏態(tài)系數(shù)，從采樣點(diǎn)A到采樣點(diǎn)B可構(gòu)成如表1所列的8種基本類型；還可用基本類型進(jìn)一步形成復(fù)合的類型，如表1中的類型1與類型2之和就是一種復(fù)合類型。

表1 用3種粒度參數(shù)所構(gòu)成的粒徑趨勢(shì)類型

上述任何一種類型的粒徑趨勢(shì)都可用一個(gè)矢量來(lái)表示，稱為粒徑趨勢(shì)矢量。該矢量的方向是從采樣點(diǎn)A指向B，其大小定義為一個(gè)單位長(zhǎng)度。對(duì)于一個(gè)采樣網(wǎng)格，每一個(gè)采樣點(diǎn)的粒度參數(shù)都可以與相鄰采樣點(diǎn)的參數(shù)進(jìn)行比較，從而找出對(duì)應(yīng)于這個(gè)采樣點(diǎn)的所有粒徑趨勢(shì)矢量。粒徑趨勢(shì)矢量還有不同的定義方法，但其目的都是為了將粒徑趨勢(shì)定量化或?qū)⒘节厔?shì)表示為平面圖形。

進(jìn)一步考察粒徑趨勢(shì)矢量，可以發(fā)現(xiàn)在一個(gè)采樣網(wǎng)格中各個(gè)方向上出現(xiàn)某一類型的粒徑趨勢(shì)的概率是不同的。如果將某一采樣點(diǎn)的所有粒徑趨勢(shì)矢量相加，求出這些矢量的合矢量，則各個(gè)采樣點(diǎn)的合矢量往往構(gòu)成一種有序的分布。在各向同性的情況下，合矢量的長(zhǎng)度應(yīng)該接近于零。因此，合矢量的有序分布代表一種各向異性的特征。在天然的海洋環(huán)境中，粒徑趨勢(shì)矢量的各向異性的程度隨地點(diǎn)而異。

Gao-Collins方法，把粒徑趨勢(shì)矢量的平面分布圖看成一幅同時(shí)包含信息和噪聲的圖像，從而用圖像處理技術(shù)來(lái)提取平面二維粒徑趨勢(shì)矢量圖像中所含的沉積物輸運(yùn)信息。對(duì)于所考慮的海域可通過(guò)各采樣點(diǎn)底質(zhì)的粒度分析獲得粒徑參數(shù)的平面分布圖式，然后經(jīng)粒徑趨勢(shì)分析獲得沉積物輸運(yùn)圖式。粒徑趨勢(shì)分析首先是在采樣點(diǎn)網(wǎng)格上對(duì)每?jī)蓚€(gè)相鄰的采樣點(diǎn)進(jìn)行比較，找出所有的粒徑趨勢(shì)矢量。兩個(gè)采樣點(diǎn)是否“相鄰”，可用特征距離DΦ來(lái)衡量（DΦ通常為最大采樣間距）。如果兩個(gè)采樣點(diǎn)的實(shí)際間距小于DΦ，則判定為“相鄰”，否則判定為“不相鄰”。

然后用式（7）求出每個(gè)采樣點(diǎn)的趨勢(shì)矢量的和：

式中，n為所考慮的采樣點(diǎn)的趨勢(shì)矢量總數(shù)；r（x，y）為趨勢(shì)矢量。

2 模型開(kāi)發(fā)

2.1 環(huán)境配置

開(kāi)發(fā)平臺(tái)為Microsoft Visual Studio 2010，開(kāi)發(fā)語(yǔ)言為C#，組件配置為嵌入ArcGIS Engine組件，運(yùn)行環(huán)境為Framework 3.5和ArcGIS Engine Runtime。

2.2 功能分析

通過(guò)對(duì)ArcGIS軟件的二次開(kāi)發(fā)，實(shí)現(xiàn)數(shù)據(jù)的自動(dòng)化分析和可視化顯示。具體功能如圖1所示。

圖1 功能分析流程圖

2.3 代碼實(shí)現(xiàn)

模型運(yùn)行代碼：

1）初始數(shù)據(jù)轉(zhuǎn)化為Φ無(wú)量綱粒徑。

public void lijing_zhuanhua(string []id_str,int zw_count)

{

int n = zw_count;//站位個(gè)數(shù)

int m = 93;//粒徑數(shù)量

double[,] lj = new double[n, m];//各個(gè)站位粒徑

double[,] tj = new double[n, m];//各個(gè)站位微分體積

double[,] lj_fai = new double[n, m];//各個(gè)站位粒徑轉(zhuǎn)化成Φ值

double[,] tj_pinlv = new double[n, m];//各個(gè)站位微分體積轉(zhuǎn)成頻率

string []name_str=new string[n];

for (int i = 0; i ＜ name_str.Length; i++)

{

name_str[i] = id_str[i];

}

string sql = string.Empty;

DataTable dt=null;

string strConn = "Provider=Microsoft.Jet.Oledb.4.0;Data Source=" + Environment.CurrentDirectory + "\lidu.mdb";

OleDbConnection conn;

conn = new OleDbConnection(strConn);

conn.Open();

for (int j = 0; j ＜ name_str.Length; j++)

{

sql = "select * from "+name_str[j];

DataSet ds = new DataSet();

OleDbDataAdapter da=new OleDbDataAdapter(sql, conn);

da.Fill(ds);

dt = ds.Tables[0];

for (int k = 0; k ＜ 93; k++)

{

lj[j,k] = double.Parse(dt.Rows[k]["lj"].ToString());

lj_fai[j, k] =- Math.Log((lj[j,k]/1000),2);

tj[j,k] = double.Parse(dt.Rows[k]["tj"].ToString());

tj_pinlv[j, k] = tj[j, k] / 100;

}

2）平均粒徑、分選系數(shù)、偏態(tài)和峰態(tài)計(jì)算。

public void sizhi_jisuan(string []id_str,int zw_count)

{

double avg_lj=0 ;//平均粒徑

double fxxs = 0;//分選系數(shù)

double pt = 0;//偏態(tài)

double ft = 0;//峰態(tài)

double[] pingjunlijing = new double[name_str.Length];

double[] fenxuanxisu = new double[name_str.Length];

double[] piantai = new double[name_str.Length];

double[] fengtai = new double[name_str.Length];

for (int o = 0; o ＜ name_str.Length; o++)

{

for (int p = 0; p ＜ 93; p++)

{

avg_lj += ((lj_fai[o, p] * tj_pinlv[o, p])/100);

}

for (int f1 = 0; f1 ＜ 93; f1++)

{

fxxs += (Math.Pow((lj_fai[o, f1] - avg_lj), 2) * tj_pinlv[o, f1])/100;

pt += (Math.Pow((lj_fai[o, f1] - avg_lj), 3) * tj_pinlv[o, f1])/100;

ft += (Math.Pow((lj_fai[o, f1] - avg_lj), 4) * tj_pinlv[o, f1])/100;

}

pingjunlijing[o] = avg_lj;

fenxuanxisu[o] = Math.Sqrt(fxxs);

if (pt ＜ 0)

{

piantai[o] = (-pt)/Math.Pow(fxxs, 3);

}

else

{

piantai[o] = pt/Math.Pow(fxxs, 3);;

}

fengtai[o] = ft/Math.Pow(fxxs,4);

avg_lj = 0;

fxxs = 0;

pt = 0;

ft = 0;

}

add_tezhengzhi(pingjunlijing, "pjlj");

add_tezhengzhi(fenxuanxisu, "fxxs");

add_tezhengzhi(piantai, "pt");

add_tezhengzhi(fengtai, "ft");//將值添加到屬性表中

}

3）半方差函數(shù)模型。

public void bianfangcha()

{

IFeatureLayer pfeaturelayer;

IFeatureClass pfeatureclass;

ITable ptable;

IQueryFilter pfilter = new QueryFilterClass();

int u2 = 0;

if (axMapControl1.LayerCount == 0)

{

MessageBox.Show("您沒(méi)有添加數(shù)據(jù)，請(qǐng)?zhí)砑觭hp文件");

return;

}

pfeaturelayer = axMapControl1.Map.get_Layer(0) as IFeatureLayer;

pfeatureclass = pfeaturelayer.FeatureClass;

ptable = pfeatureclass as ITable;

u2 = ptable.RowCount(pfilter);

double[] x_1 = new double[u2];

double[] y_1 = new double[u2];

double[,] distance_xy = new double[u2,u2];

double[] pjlj_1 = new double[u2];

double[] pt_1 = new double[u2];

double[] fxxs_1 = new double[u2];

double[] ft_1 = new double[u2];

string[] str_zw = new string[u2];

double[] bianfangcha1=new double [10];

double[] bianfangcha2=new double [10];

double[] bianfangcha3=new double [10];

ICursor pcursor = ptable.Search(null, false);

IRow prow = pcursor.NextRow();

for (int qq = 0; qq ＜ u2; qq++)

{

x_1[qq] = double.Parse(prow.get_Value(ptable. FindField("pmx")).ToString());

y_1[qq] = double.Parse(prow.get_Value(ptable. FindField("pmy")).ToString());

pjlj_1[qq] = double.Parse(prow.get_Value(ptable. FindField("pjlj")).ToString());

pt_1[qq] = double.Parse(prow.get_Value(ptable. FindField("pt")).ToString());

fxxs_1[qq] = double.Parse(prow.get_Value(ptable. FindField("fxxs")).ToString());

ft_1[qq] = double.Parse(prow.get_Value(ptable. FindField("ft")).ToString());

str_zw[qq] = prow.get_Value(ptable.FindField("zw_id")). ToString();

prow = pcursor.NextRow();

}

for (int fa = 0; fa ＜ u2; fa++)

{

for (int fb = 0; fb ＜ u2 ; fb++)

{

distance_xy[fa,fb] = Math.Sqrt((x_1[fa] - x_1[fb]) * (x_1[fa] -x_1[fb]) + (y_1[fa] - y_1[fb]) * (y_1[fa] - y_1[fb]));

}

bianfangcha1=bianfangcha_jisuan(distance_xy,pjlj_1);

bianfangcha2=bianfangcha_jisuan(distance_xy,pt_1);

bianfangcha3=bianfangcha_jisuan(distance_xy,fxxs_1);

//調(diào)用添加字段函數(shù)

}

public double[] bianfangcha_jisuan(double [,]distance_xy, double []fenlei)

{

double[] di = { 300, 500, 800, 1000, 1200, 1500, 2000, 2500, 2800, 3000 };

double count=0;

double bianfangcha=0;

double []bianfangcha_shuzu=new double [10];

for (int i = 0; i ＜ di.length;i++)

{

for (int zu = 0; zu ＜ fenlei.length; zu++)

{

for (int zo = 0; zo ＜ fenlei.length; zo++)

{

if (distance_xy[zu, zo] ＞ 0 && distance_xy[zu, zo] ＜ di[i])

{

count++;

bianfangcha+=Math.Pow((fenlei[zu]-fenlei[zo]),2);

}

bianfangcha_shuzu[i]=bianfangcha/(2*count);

}

return bianfangcha;

}

4）Gao-Collins方法沉積物粒徑趨勢(shì)分析模型。

①每個(gè)采樣點(diǎn)的趨勢(shì)矢量和。

public void qushi_tezhengjuli_fenxi(double [,]distance_ xy,double []pjlj_1,double []fxxs_1,double []pt_1,double tezhengjuli,int u2,double[,]distance_x_yipie_1,double[,]distance_ y_yipie_1)

{

for (int zu = 0; zu ＜ u2; zu++)

{

for (int zo = 0; zo ＜ u2; zo++)

{

if (distance_xy[zu, zo] ＞ 0 && distance_xy[zu, zo] ＜= tezhengjuli)

{

if (fxxs_1[zu] ＜ fxxs_1[zo] && pjlj_1[zu] ＜ pjlj_1[zo] && pt_1[zu] ＞ pt_1[zo])

{