劉吉波，王志紅

(貴州工程應用技術學院 礦業(yè)工程學院，貴州 畢節(jié) 551700)

等值線圖是一種非常重要的專題圖型，常用的等值線圖有等高線圖、等深線圖、等溫線圖、等壓線圖等。由離散點生成等值線，可采用Delaunay三角剖分算法根據離散點構建三角網，然后插值生成等值線，這種算法應用較廣[1-5]。對于規(guī)則的格網狀數據，由于其數據結構的特殊性，可不用建立三角網而直接生成等值線。這種方法直觀，易編程實現(xiàn)，也有很多學者研究并實現(xiàn)[6-11]。本文結合開采沉陷數據處理特殊性，充分利用對角線上的數據，提高等值線精度，實現(xiàn)規(guī)則格網數據等值線的快速生成。

1 數據結構設計

設預計點有N行M列，分別位于網格點上，其位置由所在的行號和列號確定，如(i，j)。通過行列號的引用，每個矩形網格可自動分割成兩個三角形，實際并未實施剖分，故稱之為虛擬剖分，可大幅提高算法效率，如圖1所示。

圖1 規(guī)則網格點的虛擬三角剖分

設等值距為V，則可根據網格角點數據內插生成每條邊(包括對角線)位于等值線上的點。當網格角點屬性值恰好為等值距V的整數倍時，亦即網格點也是等值線上的點，將使等值線構造變復雜。為了降低算法復雜性，將該角點屬性值加上一個不影響等值線生成精度的極小值ε，則可避免等值線過網格頂點時需多次進行路徑判斷問題，且不影響等值線精度。

對每個網格內的數據進行以下約定：對于網格(i，j)，i、j為其左下角點的行號和列號，上部三角形編號為Ⅰ，下部三角形編號為Ⅱ，四條邊及對角線的編號分別為①、②、③、④、⑤，每個網格編碼方式相同，如圖2所示。

圖2 網格(i，j)內三角形和邊編號規(guī)則

2 等值線生成算法

2.1 起始三角形的確定

從編號為(0，0)的網格開始遍歷，找到第一個需要繪制等值線的三角形所在網格。首先計算三角形三條邊上當前等值線內插值的個數n，n可能為1，2或3。(Bb、Be)分別用來標記一個三角形內等值線段起點和終點所在的邊的編號。當n=1時，Bb=Be=插值點所在的邊的編號；當n=2時，第一段等值線的可能編號方式為(①，②)、(②，⑤)、(⑤，①)(對于△Ⅰ)或(③，④)、(④，⑤)、(⑤，③)(對于△Ⅱ)，每個三角形有3種情形，若為逆時針，則將Bb和Be互換；當n=3時，則直接將這三個點連接構成閉合等值線。

2.2 下一個三角形的確定

如何尋找下一個要繪制的三角形是生成等值線算法的關鍵。當Be=①時，下一個要繪制的三角形位于網格(i，j-1)，編號為Ⅱ；當Be=②時，下一個要繪制的三角形位于網格(i+1，j)，編號為Ⅱ；當Be=③時，下一個要繪制的三角形位于網格(i，j+1)，編號為Ⅰ；當Be=④時，下一個要繪制的三角形位于網格(i-1，j)，編號為Ⅰ；當Be=⑤，位于△Ⅰ時，下一個要繪制的三角形位于網格(i，j)，編號為Ⅱ；當Be=⑤，位于△Ⅱ時，下一個要繪制的三角形位于網格(i，j)，編號為Ⅰ。

確定了下一個三角形位置后則可計算出新三角形中另一條等值線插值點所在的邊，存儲該點位置信息，并將其所在的邊的編號賦值給Be。

重復上述過程，即可實現(xiàn)等值線的自動追蹤，從而完成等值線繪制。

2.3 等值線追蹤停止條件判斷

一條等值線追蹤停止或繪制完成的條件分以下2種情況：一種是當前等值線段的終點所在邊的編號Be等于整條等值線的起點所在邊的編號Bb，則停止追蹤，此時等值線閉合，繪制完成，如圖3中等值線A所示；另一種情況是當前等值線段的終點Be所在的邊為網格數據的邊界線時，停止追蹤。對于后一種情況，多數情況下只是繪制了半條等值線。此時，需將已生成的等值線反向存儲，并令Be=Bb，重復過程(2.2)和(2.3)，生成等值線的剩余部分，從而完成整條等值線的繪制，此時等值線不閉合，如圖3中等值線B所示。

圖3 等值線跟蹤

當生成的等值線為逆時針存儲時，需將等值線端點逆序重新存儲，便于統(tǒng)一管理。

2.4 特殊等值線處理

在開采沉陷數據處理中，下沉邊界等值線是由下沉值為10 mm的點構成，而不是下沉值為0 mm的點。同時，除下沉，傾斜、水平移動、曲率等變形都有正負之分(極少數情況下下沉也會出現(xiàn))，需要對正負零變形邊界分別處理。

對于下沉邊界，可以生成下沉值為10 mm的等值線，將其屬性賦值為0 mm即可；而對于其他變形的正負0邊界，可以用±ε(ε是個極小值)這兩條特殊等值線代替，其屬性值分別賦值為±0。

2.5 程序流程圖

根據以上算法設計，由預計數據生成等值線的程序流程圖如圖4所示。

圖4 程序流程圖

3 等值線平滑

由于預計效率要求，預計網格點密度總是有限的，導致直接生成的等值線會呈現(xiàn)鋸齒狀，故需進行平滑處理。在要求不高的情況下，可采用質心平滑算法。即對于相鄰三點構成的三角形，中間點的坐標用三角形質心坐標代替。采用平滑算法處理等值線簡單且效果較好。

4 算法應用

在進行煤礦開采沉陷預計時，通常將預計點布設成網格狀，根據預計網格點的變形值繪制相應的變形曲線。根據本文闡述的等值線生成算法，采用VC++ 6.0平臺，即可實現(xiàn)規(guī)則網格數據等值線自動生成，等值線生成的函數代碼如下，程序界面如圖5所示。

圖5 等值線繪制程序界面

void DrawContours(CFile ConFile)

{

int i,j,k,Flag,oneconSize;

double Value;

SP ConP;

CArraytemparray;

CString Str1,ConLable,Formats;

GetBeginRowCol(0,0);

for(i=0;i

for(j=0;j<2*(nx-1);j++)

{

GetBeginRowCol(beginrow,begincol);

i=beginrow;

j=begincol;

if(j==2*(nx-1))

break;

ConP=GetFirstP(i,j);

Formats.Format("%%.%df ",3);

ConLable.Format(Formats,ConP.v);

Value=ConP.v;

onecon.RemoveAll();

flagone=0;

flagclose=0;

if(DrawFirstTri(i,j,Value))

{

Flag=DrawNextTri(i,j,Value);

while(Flag)

Flag=DrawNextTri(nr,nc,Value);

if((IsEqual(onecon.GetAt(onecon.GetSize()-1).x,onecon.GetAt(0).x))&&(IsEqual(onecon.GetAt(onecon.GetSize()-1).y,onecon.GetAt(0).y)))

flagclose=1;

if((flagclose!=1)&&(flagone!=1))

{

be=bb;

temparray.RemoveAll();

for(k=onecon.GetSize()-1;k>=0;k--)

temparray.Add(onecon.GetAt(k));

onecon.RemoveAll();

if(be==1)

{

GetPtAndSet0(grids[i][GridCol(j)].pny,grids[i][GridCol(j)].pointsy,Value,1,&ConP);

SetGridFlag(i,j);

if(j-1>=0)

SetGridFlag(i,j-1);

}

if(be==2)

{

GetPtAndSet0(grids[i][GridCol(j)].pnx,grids[i][GridCol(j)].pointsx,Value,1,&ConP);

SetGridFlag(i,j);

if(i-1>=0)

SetGridFlag(i-1,j-1);

}

if(be==3)

{

GetPtAndSet0(grids[i][GridCol(j)+1].pny,grids[i][GridCol(j)+1].pointsy,Value,1,&ConP);

SetGridFlag(i,j);

if(j+1<2*(nx-1))

SetGridFlag(i,j+1);

}

if(be==4)

{

GetPtAndSet0(grids[i+1][GridCol(j)].pnx,grids[i+1][GridCol(j)].pointsx,Value,1,&ConP);

SetGridFlag(i,j);

if(i+1

SetGridFlag(i+1,j+1);

}

if(be==5)

{

GetPtAndSet0(grids[i][GridCol(j)].pnxy,grids[i][GridCol(j)].pointsxy,Value,1,&ConP);

SetGridFlag(i,j);

if(Odd(j))

SetGridFlag(i,j-1);

else

SetGridFlag(i,j+1);

}

Flag=DrawNextTri(i,j,Value);

while(Flag)

Flag=DrawNextTri(nr,nc,Value);

for(k=0;k

temparray.Add(onecon.GetAt(k));

onecon.RemoveAll();

for(k=0;k

onecon.Add(temparray.GetAt(k));

}

if(flagclose==1)

oneconSize=onecon.GetSize()-1;

else

oneconSize=onecon.GetSize();

if(oneconSize>=2)

{

ConFile.Write(ConLable,ConLable.GetLength());

if(flagclose==1)

{

for(k=0;k

{

Str1.Format("%.3f,%.3f ",onecon.GetAt(k).x,onecon.GetAt(k).y); ConFile.Write(Str1,Str1.GetLength());

}

Str1="C ";

ConFile.Write(Str1,Str1.GetLength());

}

else

for(k=0;k

{

Str1.Format("%.3f,%.3f ",onecon.GetAt(k).x,onecon.GetAt(k).y);

ConFile.Write(Str1,Str1.GetLength());

}

}

}

else

{

if(onecon.GetSize()>=2){

ConFile.Write(ConLable,ConLable.GetLength());

for(k=0;k

{

Str1.Format("%.3f,%.3f ",onecon.GetAt(k).x,onecon.GetAt(k).y);

ConFile.Write(Str1,Str1.GetLength());

}

Str1="C ";

ConFile.Write(Str1,Str1.GetLength());

}

}

}

}

以某煤礦開采沉陷預計為例，單一工作面開采，預計網格數為40×50個，可繪制下沉等值線如圖6所示。

圖6 某煤礦開采沉陷預計下沉等值線圖

5 結 論

本文闡述了根據規(guī)則格網數據快速生成開采沉陷預計移動變形等值線的生成算法原理，并利用VC++ 6.0編程實現(xiàn)。研究表明，根據文中算法，程序編寫較其他算法更為簡單, 生成的等值線圖效果較好，可用于其他規(guī)則形式數據處理。該算法主要處理對象為規(guī)則數據，具有專用性，適用范圍受到限制。對于任意離散數據，可以由線性內插法生成規(guī)則數據后再應用本算法。