GIS空间分析51.docx
- 文档编号:7544876
- 上传时间:2023-01-24
- 格式:DOCX
- 页数:16
- 大小:79.49KB
GIS空间分析51.docx
《GIS空间分析51.docx》由会员分享,可在线阅读,更多相关《GIS空间分析51.docx(16页珍藏版)》请在冰豆网上搜索。
GIS空间分析51
空间形态
曲面样条函数插值方法原理
第一节曲面样条函数插值计算程序(SSPL2A.FOR)
一、程序功能
本程序使用曲面样条函数插值计算二元空间曲面。
在一个指定的矩形区域内,若已有一组空间曲面观测数据(xi,yi,zi)(i=1,2,...,n),如钻孔数据,利用本程序可推断空间曲面(如地层面)的曲面样条函数表达式z=f(x,y)并得到网格化数据。
曲面样条函数插值计算的特点是:
已知点坐标不必在平面上呈规则矩形排列,所得曲面函数任意阶可微;故尤为适合推断地质曲面。
除计算待估区域的曲面样条函数值(Z值)外,本程序还可计算该空间曲面在X、Y方向的一价偏导数值、二价偏导数值和梯度值。
本程序将计算得到的网格化数据按GRAPHER系列绘图软件中.GRD文件格式存于一数据文件,所以本程序的结果可直接供TOPO、SURF模块使用,产生空间曲面的等值线图和三维立体图。
二、数学模型及实现方法
1.曲面样条函数表达式的求解
曲面样条函数可看成是一块无限大平板纯弯曲时的变形。
根据弹性力学理论,联系挠度(变形)W(x,y)和作用于该板上负载q(x,y)之间的微分方程是
D4W=q(1-1)
其中D为板的抗弯刚度。
若将该板的变形看作是由n个给定的独立点(xi,yi,i=1,2,...,n)上的点负载pi(i=1,2,...,n)决定的,引入极坐标,并采用边界条件:
在远离负载的地方,曲面样条函数是"平坦"的或"平面状的",则可推导出曲面样条函数的表达式为:
(1-2)
式中:
a0,a1,a2,Fi(i=1,2,...,n)为待定系数;
ri2=(x-xi)2+(y-yi)2;
调节曲面曲率大小的经验参数。
当曲面曲率变化较大时,要取得小些,反之取大些。
一般对平坦曲面取=1~10-2,对有奇性的曲面取=10-5~10-6。
上式中的n+3个未知数a0,a1,a2,Fi(i=1,2,...,n)可以通过下列方程组求得:
(1-3)
式中,cj=16D/kj。
kj是关于点j的弹性常数,若kj,则cj0。
一般在地质曲面插值中都取cj=0,以使得求出的曲面样条函数在已知点吻合原始数据。
该方程组的矩阵表示式为:
(1-4)
式中
这是一个对称方程组。
在地质应用中,一般cj(j=1,2,...,n)都取为零,此时矩阵A的主对角线元素全为零。
为保证解的稳定性,可用Householder变换法解该方程组。
此法是先用镜象变换将A化为三对角对称矩阵,再用追赶法三对角对称方程组。
曲面方程确定之后,还可以根据需要确定曲面方向。
曲面方向可用曲面上任意一点的法向量描述。
法向量可用一阶偏导数表示,其计算公式为:
(1-5)
(1-6)
同理可以计算二价偏导数和梯度
(1-7)
(1-8)
(1-9)
2.网格化数据的计算
利用上面得到的曲面样条函数表达式,代入每一点的x、y值,就可得到该点的z值。
3.网格化数据文件(.GRD)的格式‘
将网格化数据写成TOPO、SURF能接受的格式,就能使用TOPO、SURF绘制等值线图和三维立体图。
文件的头五行为:
form四字节字符常数,文本文件为DSAA,二进制文件为DSBB;
nx,ny网格列数、排数;
xlo,xhix坐标的最小值、最大值;
ylo,yhiy坐标的最小值、最大值;
zlo,zhiz坐标的最小值、最大值。
此后为网格中各排(固定y坐标)z的数据,各排中x坐标从小到大排列;第一排的y坐标最小。
若一排数据一行写不下,可换行写,但每一排均需从新行开始。
如是文本文件,采用自由格式。
三、程序说明
1.主要标识符说明
LDATA原始数据数组最大行数;
JDATA原始数据数组最大列数;
JGRD网格最大列数;
N观测值个数;
DATA二维数组,存原始数据观测值;
X二维数组,存从原始数据中挑选出的X、Y、Z值;
ZKM字符型数组,存钻孔名;
EPS即;
XMAX、XMIN、YMAX、YMIN网格化区域(插值区)x的上、下界,y的上、下界;
DX、DYx方向、y方向的网格间距;
F一维数组,存Fi(i=1,2,...,n)、a0、a1、a2;
A一维大数组,存放矩阵A,;如果X和ZKM、F数组扩大,本数组也应相应扩大,其最小设置为(n+3)(n+4)/2。
2.子程序功能说明
SSPL1生成矩阵A,并调用HOUSE子程序解之;
HOUSE即用HOUSEHOLDER法求解逆矩阵A-1;
SSPLD0利用曲面样条函数表达式计算待估点的z值;
SSPLD1利用曲面样条函数表达式计算待估点x方向的一价偏导数值;
SSPLD2利用曲面样条函数表达式计算待估点y方向的一价偏导数值;
SSPLD3利用曲面样条函数表达式计算待估点x方向的二价偏导数值;
SSPLD4利用曲面样条函数表达式计算待估点y方向的二价偏导数值;
SSPLD5利用曲面样条函数表达式计算待估点的梯度值。
3.输入输出
1)文件
输入文件即观测点数据文件,每个观测点一行,自由格式。
若有钻孔名,则放在第一列,后面放x、y、z等数据,这些数据列排列无限制。
输出文件即网格化数据文件,文本格式。
2)屏幕输出
输出曲面样条函数系数即Fi(i=1,2,...,n)、a0、a1、a2。
四、操作过程
本程序采用交互式方法控制运行。
(1)程序运行后,将显示:
InputYourDataFilename
此时输入数据文件名。
如文件不能打开,要求重新输入。
(2)显示:
InputTotalNumberofColumns
此时输入需读入的数据文件的列数。
(3)显示:
DoYouHaveDrill-Name(y/n)?
如有钻孔名,则键入Y,否则键入N。
(4)显示:
ToSelectDataCoordinate(IX,IY,IZ)
此时输入X、Y、Z在数据文件中的位置(列号)。
(5)显示:
InputDataLimit(X01,X02,Y01,Y02)
此时输入参加运算的数据限制条件,即只有X01XX02和Y01YY02的数据才参加运算。
(6)显示:
Printoriginaldatamatrix(Y/N)?
若要显示原始数据矩阵,键入Y,否则键入N。
(7)显示:
InputEPS
此时输入小量ε。
(8)然后要求输入插值区边界:
InputXMIN,XMAX,YMINandYMAX
(9)接下来选择网格化数据结果类型:
InputMapType:
0-Contour1-1stOrderDerivativeinX-Direction
2-1stOrderDerivativeinY-Direction3-2ndOrderDerivativeinX-Direction
4-2ndOrderDerivativeinY-Direction5-Gradient
即选择0-待估点的Z值(可用于作等值线图),1-X方向的一价偏导数,2-Y方向的一价偏导数,3-X方向的二价偏导数,4-Y方向的二价偏导数,5-梯度值。
(10)接下来计算网格化数据,先输入结果文件名:
InputGRIDFile-Name
(11)再输入网格间距:
InputDXandDY
(12)一个网格化文件完成后,将接连显示:
AnotherGrd-file(Y/N)?
ToChangeXMIN,XMAX,YMINandYMAX(Y/N)?
ToChangeEPS(Y/N)?
ToChangeDataCoordinateorDataLimit(Y/N)?
第一个响应可以改变DX、DY,即返回到(10)继续执行;第二个响应可以改变插值区域,返回到(8)继续执行;第三个响应则可以改变EPS,返回到(7)继续执行;最后一个响应可以改变X、Y、Z在数据文件中的位置(列号)和参加运算的数据限制条件,返回到(4)继续执行。
若回答都是N,则往下执行直到结束。
第二节曲面样条函数分区插值计算程序(SSPL7A.FOR)
一、程序功能
本程序使用曲面样条函数进行分区插值计算二元空间曲面。
在一个矩形区域内,可将其分为不同的插值区,并使用不同的原始数据对这些区分别插值计算推断空间曲面,最后将这些区合并成一个网格化文件。
二、数学模型及实现方法
1.曲面样条函数插值
曲面样条函数的插值原理同SSPL2A.FOR(见第一节)。
2.分区边界坐标数据文件
各分区(可以是任意不规则单连通区域)的边界坐标组成一个闭合曲线(折线)。
分区边界坐标数据文件的格式为:
n1,0
x11,y11
x12,y12
............
x1n,y1n
n2,0
x21,y21
x22,y2
............
x2n,y2n
............
其中n1为第一个区域的边界点数,后接第一个区域边界的第一个坐标、第二个坐标、...、第n1个坐标;n2为第二个区域的边界点数,后接第二个区域边界的第一个坐标、第二个坐标、...、第n2个坐标;以下为第三个区域,...。
每个区域的第一个坐标与最后一个坐标相同,以保证闭合。
各区域边界可以是逆时针连接,也可以是顺时针连接。
3.点坐标区域归属判定原理
显然,本程序的关键是对矩形区内的任一点判定其归属区,即判定其是在各分区的里边还是外边。
×
。
图1图2
如图1所示,"。
"在区域里边,"×"在区域外边,人工判定很容易,但计算机判定并非易事。
本程序采用的方法为,对每一点,作过该点的水平线,求出水平线与区域边界的交点和两条交线;由待判定点与两条交线的端点组成两个矢量,两条交线也作为两个矢量,再计算待判定点与交线的端点组成的矢量与交线矢量之间的夹角,便可判断该点是位于区域内还是区域外。
对于被多个区域包围的点,根据包围区域数是奇数还是偶数决定其是否插值,具体来讲,包围区域数是奇数就插值,偶数不插值。
如图2,画水平线的区域插值,中间圆形区不插值。
所以,本程序可适用于各种复杂的情况。
对于不插值的点,程序对其进行白化处理,即其Z值为1.70141×1038。
三、程序说明
1.主要标识符说明
LDATA原始数据数组最大行数;
JDATA原始数据数组最大列数;
JGRD网格最大列数;
KDATA0所有区域边界坐标最多个数;
KDATA单个区域边界坐标最多个数;
N观测值个数;
KNUM所有区域边界坐标个数;
KEND2分区区域数;
ZBLANK白化值,等于1.70141×1038
DATA二维数组,存原始数据观测值;
X二维数组,存从原始数据中挑选出的X、Y、Z值;
ZKM字符型数组,存钻孔名;
EPS即;
XMAX、XMIN、YMAX、YMIN网格化区域(插值区)x的上、下界,y的上、下界;
DX、DYx方向、y方向的网格间距;
F一维数组,存Fi(i=1,2,...,n)、a0、a1、a2;
A一维大数组,存放矩阵A,;如果X和ZKM、F数组扩大,本数组也应相应扩大,其最小设置为(n+3)(n+4)/2;
XY0二维数组,存放所有边界坐标数据。
2.子程序功能说明
SSPL1生成矩阵A,并调用HOUSE子程序解之;
HOUSE即用HOUSEHOLDER法求解逆矩阵A-1;
SSPLD0利用曲面样条函数表达式计算待估点的z值;
SSPLD1利用曲面样条函数表达式计算待估点x方向的一价偏导数值;
SSPLD2利用曲面样条函数表达式计算待估点y方向的一价偏导数值;
SSPLD3利用曲面样条函数表达式计算待估点x方向的二价偏导数值;
SSPLD4利用曲面样条函数表达式计算待估点y方向的二价偏导数值;
SSPLD5利用曲面样条函数表达式计算待估点的梯度值;
PJUDG0读入并整理分区区域边界坐标数据;
PJUDG1通过调用PJUDG2判断某一点与所有分区区域的关系;
PJUDG2判断某一点与某一分区区域的关系。
3.输入输出
1)文件
输入文件有两个,一为观测点数据文件,数据格式同第一节;二为分区区域边界数据文件,数据格式见上。
输出文件即网格化数据文件,文本格式。
2)屏幕输出
输出曲面样条函数系数即Fi(i=1,2,...,n)、a0、a1、a2。
四、操作过程
本程序采用交互式方法控制运行。
(1)程序运行后,将显示:
InputYourDataFilename
此时输入数据文件名。
如文件不能打开,要求重新输入。
(2)显示:
InputTotalNumberofColumns
此时输入需读入的数据文件的列数。
(3)显示:
DoYouHaveDrill-Name(y/n)?
如有钻孔名,则键入Y,否则键入N。
(4)显示:
ToSelectDataCoordinate(IX,IY,IZ)
此时输入X、Y、Z在数据文件中的位置(列号)。
(5)显示:
InputDataLimit(X01,X02,Y01,Y02)
此时输入参加运算的数据限制条件,即只有X01XX02和Y01YY02的数据才参加运算。
(6)显示:
Printoriginaldatamatrix(Y/N)?
若要显示原始数据矩阵,键入Y,否则键入N。
(7)显示:
InputEPS
此时输入小量ε。
(8)显示:
PleaseSelectDestination-FileType:
0--NEWGRIDFile,1--OLDGRID
FILE,2--None
选择计算类型,0-建立一个新的网格化文件,1-将计算结果叠加到一个已存在的网格化文件上,2-不进行插值计算,转至(14)。
(9)显示
InputGRIDFile-Name
输入结果文件名。
当(8)中选择1时,以下3步不执行。
(10)接下来选择网格化数据结果类型:
InputMapType:
0-Contour1-1stOrderDerivativeinX-Direction
2-1stOrderDerivativeinY-Direction3-2ndOrderDerivativeinX-Direction
4-2ndOrderDerivativeinY-Direction5-Gradient
即选择0-待估点的Z值(可用于作等值线图),1-X方向的一价偏导数,2-Y方向的一价偏导数,3-X方向的二价偏导数,4-Y方向的二价偏导数,5-梯度值。
(11)然后要求输入插值区边界:
InputXMIN,XMAX,YMINandYMAX
(12)再输入网格间距:
InputDXandDY
(13)输入分区区域边界坐标数据文件名:
InputBoundarydatafilename
(14)一个网格化文件完成后,将接连显示:
AnotherDestination-FileType(Y/N)?
ToChangeEPS(Y/N)?
ToChangedatacoordinateorDataLimit(Y/N)?
第一个响应可以选择计算类型重新计算,即返回到(8)继续执行;第二个响应可以改变EPS,返回到(7)继续执行;最后一个响应可以改变X、Y、Z在数据文件中的位置(列号)和参加运算的数据限制条件,返回到(4)继续执行。
若回答都是N,则往下执行直到结束。
第三节地质曲面边界追踪程序(BLNGRDA.FOR)
一、程序功能
本程序通过对第二节生成的GRD文件进行操作,对地质曲面边界进行追踪,生成边界坐标数据文件。
所谓地质曲面边界,指某一地层缺失区的边界。
在第二节的分区插值处理中,对未插值的点进行了白化处理,故这些未插值的点组成的区域可看作地层缺失区。
若将某些特殊区域(如断层区)进行白化处理,利用本程序可追踪这些区域的边界。
二、实现方法
本程序的具体操作比较复杂,主要思路是通过对每一点与其各个方向邻近点的关系,勾勒出各个地层缺失区的边界。
三、程序说明
1.主要标识符说明
JGRD网格最大列数;
JDATA2所有边界点之和最大数;
JDATA3缺失区最多个数。
ZBLANK白化值,等于1.70141×1038
2.输入输出
1)输入
输入文件为待处理的网格化数据文件,文本格式。
2)输出
输出文件为边界坐标数据文件,文本格式。
三、操作过程
本程序采用交互式方法控制运行。
(1)程序运行后,将显示:
InputGrd-fileName
此时输入网格化文件名。
(2)然后将显示:
InputDestinationFileName
此时输入结果文件名。
第四节网格化文件Z值限制程序(GRD1A.FOR)
一、程序功能
本程序对GRD文件进行Z值限制处理操作,以得到符合实际的地质曲面Z值数据。
由于在插值处理中,资料空白区的Z值是推断产生的,有可能偏高或偏低,与实际情况严重不符。
利用本程序可以将Z值限制在一定范围内,将超出Z值范围的区域赋白化值,即进行白化处理。
二、程序说明
1.主要标识符说明
JGRD网格最大列数;
ZBLANK白化值,等于1.70141×1038
2.输入输出
1)输入
输入文件为待处理的网格化数据文件,文本格式。
2)输出
输出文件为处理后的网格化数据文件,文本格式。
三、操作过程
本程序采用交互式方法控制运行。
(1)程序运行后,将显示:
InputGrd-fileName
此时输入网格化文件名。
(2)显示:
InputDestinationFileName
此时输入结果文件名。
(3)然后将显示:
InputMinimumDepthandMaximumDepth
此时输入该地质曲面的最小深度值和最大深度值。
第五节网格化文件叠加处理程序(GRD21A.FOR)
一、程序功能
本程序对两个GRD文件进行叠加处理,生成一个新的网格化数据文件。
如若第一个网格化数据文件代表A地层的标高,另一个网格化文件代表B地层与A地层之间的厚度,新产生的网格化文件则代表B地层的标高。
本程序可以将第二个网格化数据的值限制在一定范围内,将超出该值范围的区域进行白化处理。
二、程序说明
1.主要标识符说明
JGRD网格最大列数;
ZBLANK白化值,等于1.70141×1038
2.输入输出
1)输入
输入文件为两个待处理的网格化数据文件,文本格式。
2)输出
输出文件为处理后的网格化数据文件,文本格式。
三、操作过程
本程序采用交互式方法控制运行。
(1)程序运行后,将显示:
Input1stGrd-fileName
此时输入第一个网格化文件名。
(2)显示:
Input2ndGrd-fileName
此时输入第二个网格化文件名。
(3)显示:
InputDestinationFileName
此时输入结果文件名。
(4)显示:
InputMinimumDepthandMaximumDepthof2ndGrd-file
此时输入第二个网格化文件Z的最小值和最大值,用以限制第二个地层的Z值(如地层厚度值)的范围。
(5)然后将显示:
InputSequence(1=Grd1:
Upper,-1=Grd2:
Upper)
此时决定两个地层的层序关系,若第一个地层在上,输入1,反之输入-1。
第六节网格化文件协调处理程序(GRD22A.FOR)
一、程序功能
本程序对两个GRD文件进行协调处理,生成一个新的网格化数据文件来代替第二个网格化文件。
假设第一个网格化数据文件代表A地层的标高,第二个网格化文件代表B地层的标高,且A地层在B地层之上,如果在某一点B地层的标高大于A地层的标高,则两个地层在这一点上不协调(不符合地质实际情况)。
本程序将第一个GRD文件作为基准文件,将第二个GRD文件中的不协调点进行白化处理,即看作地层缺失点。
二、程序说明
1.主要标识符说明
JGRD网格最大列数;
ZBLANK白化值,等于1.70141×1038
2.输入输出
1)输入
输入文件为两个待处理的网格化数据文件,文本格式。
2)输出
输出文件为处理后的网格化数据文件,文本格式。
三、操作过程
本程序采用交互式方法控制运行。
(1)程序运行后,将显示:
Input1stGrd-fileName
此时输入第一个网格化文件名。
(2)显示:
Input2ndGrd-fileName
此时输入第二个网格化文件名。
(3)显示:
InputDestinationFileName
此时输入结果文件名。
(4)然后将显示:
InputSequence(1=Grd1:
Upper,-1=Grd2:
Upper)
此时决定两个地层的层序关系,若第一个地层在上,输入1,反之输入-1。
第七节剖面图数据生成程序(GRDSECA.FOR)
一、程序功能
本程序从GRD文件中提取剖面图数据。
剖面可以沿任意联线方向前进,即可以是多条折线组成的剖面。
二、程序说明
1.主要标识符说明
JGRD网格最大列数;
KDATA控制剖面线的最多点数;
KDATA1剖面线截取的最多网格数;
ZBLANK白化值,等于1.70141×1038。
2.输入输出
1)输入
输入文件有两个,一为网格化数据文件,文本格式;二为控制剖面线的点坐标数据,文件格式同上述的边界坐标数据文件,但本程序一次只处理一条剖面线。
2)输出
输出文件为剖面线坐标和获取的地层曲面(标高)数据,每一行代表一个点,共有6列:
第1列--距剖面线起始点的距离;
第2列--该点的地质曲面(标高)值;
第3列--该点的X坐标;
第4列--该点的Y坐标;
第5列--该点在网格化文件中的行号;
第6列--该点在网格化文件中的列号。
利用第1、第2列就可作出剖面图。
三、操作过程
本程序采用交互式方法控制运行。
(1)程序运行后,将显示:
Input1stGrd-fileName
此时输入第一个网格化文件名。
(2)显示
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- GIS 空间 分析 51