krigingmethod.pro

在IDL环境下

;kriging interpolation Method
;ysh 20090313 start


;---------------------------------------------------------------------------------------------
;导入已知点的坐标
;格式：xi,yi,zi(X坐标，Y坐标，Z高程)
;返回所得点坐标的一个 N * 3 的数组以及N的值

Function KrgInitPointInput

	FileName = DIALOG_PICKFILE(FILTER = "*.txt", /MUST_EXIST,   $
	      						PATH = "F:\IDL", TITLE ="Open", $
	      						FILE = "line.txt")

	if not FILE_TEST(FileName) then Begin
		FailMes = DIALOG_MESSAGE("Fail To Open File")
		return,0
	End


	;将数据文件写入矩阵
	OPENR, fp, FileName, /Get_lun
	pp = ''
	while not EOF(fp) do Begin
		READF, fp, pp
		cc = STRSPLIT(pp, ",", /EXTRACT)
		;store the data follow by column
		if N_ELEMENTS(cc) eq 1 then Begin
	    	cc = STRSPLIT(pp," ", /EXTRACT)
	    	if N_ELEMENTS(cc) eq 3 then Begin
				if N_ELEMENTS(Data) eq 0 then Begin
					Data = Transpose([double(cc[0]), double(cc[1]), double(cc[2])])
				Endif else Begin
					Data = [Data,Transpose([double(cc[0]), double(cc[1]), double(cc[2])])]
				Endelse
	   		 Endif
		Endif
	EndWhile
	Free_lun, fp

	Data = Transpose(Data)

;	绘制面x，y，z的坐标数组
	dXCoodinates = data[0, *]
	dYCoodinates = data[1, *]
	dZCoodinates = data[2, *]

	;存储点阵中的最大最小值，用于网格边界的确定
	Common gKrgCommon, GridBorder
	GridBorder = [Min(dXCoodinates), Max(dXCoodinates), Min(dYCoodinates), Max(dYCoodinates)]

	return, Data
End

;---------------------------------------------------------------------------------------------
; Kriging 变差函数的求取
;此处先得到距离矩阵

Function KrgVariogram, nPointNum, DataMatrix

	;创建点间距离数组
	dPointDistArr = INDGEN(nPointNum, nPointNum, /DOUBLE)

	;计算点数组中各个点之间的距离，并存入 dPointDistArr
	for i = 0, nPointNum - 1 do Begin
		for j = 0, nPointNum - 1 do Begin
			dPointDistArr[i, j] = Sqrt((DataMatrix[1,j] - DataMatrix[1,i]) * (DataMatrix[1,j] - DataMatrix[1,i]) $
								  + (DataMatrix[0,j] - DataMatrix[0,i]) * (DataMatrix[0,j] - DataMatrix[0,i]))
		End
	End
	dPointDistArr = Transpose(dPointDistArr)

	return, dPointDistArr
End

;---------------------------------------------------------------------------------------------
;根据事先设定的值做出网格矩阵

Function KrgGrid, nXDiv, nYDiv

	;------------------------------------------------------------------------------------------
	;网格矩阵的生成，亦是一个 N * 3 的矩阵，其中 [N,0],[N,1] 为 x，y 坐标,[N,2]则为所求的估算值
	Common gKrgCommon, GridBorder

	nXCoodinate = (GridBorder[1] - GridBorder[0]) / nXDiv
	nYCoodinate = (GridBorder[3] - GridBorder[2]) / nYDiv

	GridMatrix = INDGEN(3, (nXDiv + 1) * (nYDiv + 1), /DOUBLE)
	for i = 0, nXDiv do Begin
		for j = 0, nYDiv do Begin
				GridMatrix[0,j + (nYDiv + 1) * i] = GridBorder[0] + nXCoodinate * j
				GridMatrix[1,j + (nYDiv + 1) * i] = GridBorder[2] + nYCoodinate * i
				GridMatrix[2,j + (nYDiv + 1) * i] = 0
		Endfor
	Endfor

	return, GridMatrix
End

;---------------------------------------------------------------------------------------------
;Kriging插值算法系数的求取，也即计算K，M矩阵
;nPointNum   : 数据点个数
;dStepLength : 变差函数步长值
;dPointDisArr: 各个已知点间的距离数组
;GridMatrix  : 划分的网格数组
;测试所用的是线形变差函数

;K,M矩阵
;    |r(v1,v1)    r(v1,v2)  ...  r(v1,vn)    1  |        | r(v0,v1) |
;    |r(v2,v1)    r(v2,v2)  ...  r(v2,vn)    1  |        | r(v0,v2) |
;K = |                      ...                 |    M = |   ...    |
;    |r(vn,v1)    r(vn,v1)  ...  r(vn,vn)    1  |        | r(vn,v0) |
;    |    1        1                1        0  |        |    1     |


Function CalCoefficeient,ithPoint , nPointNum, dStepLength, dPointDisArr, GridMatrix, DataMatrix

	;求K
	KMatrix = INDGEN(nPointNum + 1, nPointNum + 1, /DOUBLE)
	for i = 0, nPointNum - 1 do Begin
		for j = 0, nPointNum - 1 do Begin
			if (i NE j) then Begin
				KMatrix[i,j] = 3 * Ceil(dPointDisArr[i,j] / dStepLength)	;r(h) = 3 * h--假设的线形变差函数
			Endif else Begin
				KMatrix[i,j] = 0
			Endelse
		Endfor
		KMatrix[i,nPointNum] = 1
	Endfor
	for n = 0, nPointNum - 1 do KMatrix[nPointNum,n] = 1
	KMatrix[nPointNum,nPointNum] = 0

	;求M
	MMatrix = INDGEN(nPointNum + 1, /DOUBLE)

	MMatrix[nPointNum] = 1
	for m = 0, nPointNum - 1 do Begin
		dTemp = Sqrt((DataMatrix[1,m] - GridMatrix[1,ithPoint]) * (DataMatrix[1,m] -GridMatrix[1,ithPoint]) $
				+ (DataMatrix[0,m] - GridMatrix[0,ithPoint]) * (DataMatrix[0,m] - GridMatrix[0,ithPoint]))
		MMatrix[m] = 3 * Ceil(dTemp / dStepLength)
	Endfor

	;求得系数Kriging系数矩阵
	dCoef = Invert(KMatrix) # MMatrix
	dValueofZ = DataMatrix[2,*]
	;所得到的系数矩阵中含有u，故在所给点值的对应最后位置家一个零，使其消去
	dValueofZ = [Transpose(dValueofZ),0.0]
	dValueofZ = Transpose(dValueofZ)

	dResult = dValueofZ # dCoef
	return, dResult
End

;---------------------------------------------------------------------------------------------
;Kriging算法
;

Pro KrigingMethod

	;全局变量，边界数组的存储
	Common gKrgCommon, GridBorder

	;得到输入的数据
	DataMatrix = KrgInitPointInput()

	; N * 3 的数组，由此得数据个数 nPointNum给
	nPointNum = N_ELEMENTS(DataMatrix) / 3

	;各个点间距离数组
    dPointDisArr = KrgVariogram(nPointNum, DataMatrix)


	;根据所给定的网格划分方式确定步长
	nXDiv = 5	;待定，选择给出
	nYDiv = 5	;待定，选择给出

    ;网格矩阵的生成
    GridMatrix = KrgGrid(nXDiv, nYDiv)

	;确定变差函数步长，其中 nGroupNum >= 4
    nGroupNum = 5	;待定项，可选择输入
	dStepLength = (Max(dPointDisArr) - Min(dPointDisArr))/nGroupNum


	;kriging插值算法求得所有网格点的值

	for i = 0, ((nXDiv + 1) * (nYDiv + 1) - 1) do Begin
		GridMatrix[2,i] = CalCoefficeient(i, nPointNum, dStepLength, dPointDisArr, GridMatrix, DataMatrix)
	End
End