本发明属于地理信息技术应用领域，具体涉及一种基于地质体面图层和等高线图层，运用霍夫变换算法，实现直立岩层自动化识别和提取的方法。
背景技术：
根据构造产状类型，岩层可分为水平岩层、倾斜岩层与直立岩层。其中，直立岩层是指倾角大于85°的岩层，一般出现在构造强烈的地区。其地质界线沿其走向做直线延伸，不受地形影响。在常规的地质应用中，主要是人工通过读取地质图，寻找边界表现为平行直线形态的要素来识别直立岩层。而现有的计算机自动化读图技术主要应用于建筑图纸领域(杨华飞,杨若瑜,路通,等.人工读图机理分析及其在计算机读图中的应用[J].计算机科学,2008,35(2))，对于基于地图知识自动化读图的研究较少。地质图的读图仍然是通过专业人员结合知识规则，按需提取地质要素的方式实现。但人工识别方法，效率低，准确度因人而异，质量难以保证。
技术实现要素：
本发明所要解决的技术问题在于，克服现有技术存在的缺陷，提出了一种直立岩层自动化识别和提取的方法，利用直立岩层在地质图中出露为直线边界这一特征，基于霍夫变换的直线识别算法，提取出地质图中平直的岩层界线，并通过断层过滤、平缓界线过滤与非平行成对界线过滤来筛选组成直立岩层的界线，最终实现直立岩层的自动化识别和提取，实现直立岩层在各类地质图中的计算机自动绘制。
判别直立岩层的主要依据是平直地质界线的存在。霍夫变换作为图像处理中识别几何形状的基本方法之一，最初设计就是用于实现图像直线的识别提取。为此，本专利拟基于霍夫变换方法，通过对地质界线中平直部分的分段识别和提取，实现直立岩层的自动识别。
本专利的主要思路是：1)将地质图中的地质体面要素处理为线状的岩层界线要素，并消除多边形内岛(洞)及图幅边界的影响；2)对岩层界线要素上的点进行霍夫变换可将岩层界线要素变换为一组曲线，通过分析这组曲线，可判断当前线要素是否为满足要求的直线；3)对提取出的平直岩层界线进行断层界线过滤、平缓界线过滤和非平行成对界线的过滤，筛选出组成直立岩层的界线，并绘制直立岩层。
本发明直立岩层自动化识别和提取的方法，具体步骤如下：
步骤一、岩层边界读取及预处理
步骤1.1加载shp格式的岩层面要素图层StratumLayer，得到所有面要素集合Stratums＝{stai|i＝1，2，3，...，n}，n为岩层面要素的数量，stai是各个岩层面要素。
步骤1.2读取集合Stratums中每个岩层面要素stai的边界，存为线要素集合OriL＝{oli|i＝1，2，3，...，n}，每一个面要素stai的边界存储为一个线要素oli；每个线要素oli对应一个点集OLPi＝{opij(xij，yij)|j＝1，2，3，...，mi}，点opij为点集OLPi中的第j个点，j为点opij在点集OLPi中的下标，其坐标表示为(xij，yij)。mi为线要素oli对应的点集OLPi中的点的数量；oli存储“id”属性来记录其归属的岩层面要素stai的编号。
步骤1.3线要素集合OriL中可能存在部分线要素oli为多部件要素，对多部件要素oli进行拆分，得到新的线要素集合SinL＝{slk|k＝1，2，3，...，l}，l为拆分后岩层面要素边界的数量总和，线要素slk属性“id”继承拆分前oli线要素的属性“id”。
步骤1.4线要素集合SinL中必然存在线要素slk，其部分位于图幅边界而并非真实的岩层边界，需要对位于岩层边界的部分进行删除。
步骤1.5：遍历岩层边界删除后的线要素集合SinL＝{slk|k＝1，2，3，...，l}，若存在属性“id”相同，且空间邻接的线要素sli与slj(i,j＝1,2,3,…,l，且j≠i)，对线要素sli与slj进行合并，得到新的线要素集合FinL＝{flu|u＝1，2，3，...，p}，p为预处理完成后的岩层界线的数量总和。
步骤二、岩层平直界线提取
建立直线集合AllSL用于存放所有线要素提取出的平直界线，对线要素集合FinL中的每个线要素flu进行平直界线提取。针对单个线要素flu，创建一个空的直线集合StaightL存放flu中提取出的直线，线要素flu的平直界线提取方法如下：
步骤2.1给定一个长度标准Llimit来限制提取的直线的最短长度。
步骤2.2每个线要素flu表示均为点集FLPu＝{flpuv|v＝1，2，3，...，qu}，qu为线flu上的顶点个数；对于每个点flpuv(xuv，yuv)进行霍夫变换。霍夫变换的原理为笛卡尔坐标系上的线通过霍夫变换映射为一点，点通过霍夫变换映射为一曲线。设定角度容限Δθ作为提取的平直界线的走向限制，即Δθ越大，提取的平直界线可能越弯曲，反之，提取的平直界线则越接近直线。在地质图中纯粹的直线较少，因此，需要调整Δθ以提取出最合适的平直界线。对一个点进行霍夫变换，以角度容限Δθ为间隔，取根据公式(1)可得到对应的ρ0，ρ1，ρ2，ρ3,…,ρj，点(θ0,ρ0),(θ1,ρ1),(θ2,ρ2),(θ3,ρ3),…,(θd,ρd)可连成一条霍夫变换后的曲线HoughL{(θ0,ρ0),(θ1,ρ1),(θ2,ρ2),…,(θd,ρd)}；
ρ＝x cosθ+y sinθ (1)
则对于点flpuv，以Δθ为间隔，将其坐标(xuv，yuv)带入公式(1)，可得到点flpuv霍夫变换后的曲线HoughLuv＝{(0，xuv),(Δθ，xuvcosΔθ+yuvsinΔθ),(Δθ*2，xuvcos(Δθ*2)+yuvsin(Δθ*2)),(π，xuv)}。
步骤2.3创建一个直线点集StaightLineP，用于存放组成一条直线的点；首先记录开始下标s为1，将flpus加入集合StaightLineP。
步骤2.4将FLPu中的点依次添加进StaightLineP，每次往StaightLineP中添加新点时，即对StaightLineP进行如下操作：
1)若StaightLineP包含的点数少于3，则继续往StaightLineP中添加新点；反之，求取StaightLineP中所有点flput对应的HoughLut(t＝1,2,…,c，c为当前StaightLineP中点的数量)；
2)取ρmax为所有HoughLut中最大的ρ值，取ρmin为最小的ρ值由容限Δθ，和Δρ＝(ρmax-ρmin)/d，得到容限窗口(Δθ，Δρ)；
3)依次计算HoughLut与HoughLu(t+1)的交点pt(θt(t+1)，ρt(t+1))(t＝1,2,3,…,c-1)，得到点集InterP＝{p1(θ12，ρ12),p2(θ23，ρ23),p3(θ34，ρ34),…,pc-1(θ(c-1)c，ρ(c-1)c)}，每个交点pt代表笛卡尔坐标系中的一条直线，其坐标(θt(t+1)，ρt(t+1))为该直线霍夫变换后的坐标。判断所有交点的分布是否超过设定的容限窗口(Δθ，Δρ)，只有一个交点则不进行该判断：
若所有交点的分布不超出容限，且新加的点flput不是线flu上最后一点，则继续添加点，重复步骤2.4的1)-3)步骤；
若所有交点的分布不超出设定的容限，且新加的点flput即为线flu上的最后一点，则将StaightLineP＝{flpus，flpu(s+1)，...，flput}存成平直界线li，根据公式(2)计算其所在直线的霍夫坐标(lθx，lρx)，存储至平直界线li的属性“θ”与“ρ”；将平直界线li加入集合StaightL，完成线要素flu中直线要素的提取；
其中，lθx与lρx为平直界线li霍夫变换映射的点坐标，θt(t+1)，ρt(t+1)为交点pt的坐标，c为当前StaightLineP中点的数量。
若所有交点的分布超出设定的容限，且将点flput从StaightLineP中移除后，线SL<flpus,flpu(t-1)>的长度大等于Llimit，则将StaightLineP＝{flpus，flpu(s+1)，...，flpu(t-1)}存成平直界线li，根据公式(2)计算其所在直线的霍夫坐标(lθx，lρx)，存储至平直界线li的属性“θ”与“ρ”。将平直界线li加入集合StaightL。然后清空StaightLineP，重新记录开始下标s为t-1，将点flpu(t-1)加入StaightLineP，重复步骤2.4的1)-3)步骤。
若所有交点的分布超出设定的容限，且将点flput从StaightLineP中移除后，线SL<flpus,flpu(t-1)>的长度小于Llimit，则将点flpus从StaightLineP中移除，重新记录开始下标s为s+1，将点flpu(t-1)加入StaightLineP，重复步骤2.4的1)-3)步骤。
步骤2.5对于StaightL中的平直界线li与lj(j≠i)，若满足：顶点重合且属性“θ”差值小于容限Δθ，则将两个平直界线li与lj进行合并。
步骤2.6StaightL中的所有平直界线li的属性“id”继承线要素flu的同名属性。并将StaightL中的元素加入所有直线集合AllSL。
步骤三、直立岩层识别绘制
步骤3.1对所有线要素提取平直界线得到直线集合AllSL＝{li|i＝1，2，3，...，w}，w是对所有岩层的平直界线的数量总和。
步骤3.2对直线集合AllSL进行倾向断层过滤；对于一组相互有重叠的直线要素分别记录这组直线要素两侧的岩层“地层符号”序列，若两侧的序列一致，则认为这组平直界线组成一条断层界线，对这组平直界线进行过滤操作。
步骤3.3加载等高线图层，得到等高线集合Contour＝{clj|j＝1，2，3，...，r}，基于等高线对AllSL进行平缓界线过滤；设定夹角阈值α，对于平直界线li，顺序记录其与等高线的交点，并计算交点处等高线与界线li的夹角{lci1，lci2，...，lcik}；求取平均夹角lci，若其小于设定的夹角阈值α，AllSL＝SSL1∪SSL2∪SSL3，...，SSLg，g为过滤后仍保留平直界线的岩层数量；有成对平行的平直界线的岩层视作直立岩层，对非平行成对的平直界线进行过滤；本发明取容限Δθ为筛选平行平直界线的阈值，判断岩层SSLt中是否满足，存在成对的平直界线lx与ly(y≠x)属性“θ”差值小于容限Δθ，对不满足上述条件的平直界线进行过滤操作；
步骤3.5对于满足直立岩层条件的岩层，其平直界线lx与ly围起来的部分即为直立的，自动连接平直界线lx与ly的首尾点，绘制面要素sltai，即为直立岩层；最终得到直立岩层集合LStatum＝{lstai|i＝1，2，3，...，f}，f为识别得到的直立岩层的数量。
本发明方法，结合地质图读图的知识，利用几何计算将人工读图的过程转变为计算机读图。结合地址读图知识，利用几何方式实现直立岩层的自动化识别相对人工识别更加快捷，准确度更高。为普通人员读取地质图提供了帮助，也可满足地质建模、构造解析等学术研究以及工程选址等生产建设的需求。
附图说明
图1本发明方法的流程图。
图2实施例实验数据(某地的地形地质图)。
图3实施例岩层边界线提取效果。
图4示例点霍夫变换折线图。
图5示例直线提取效果图。
图6实施例直线提取结果图。
图7示例断层过滤示意图。
图8实施例断层过滤效果图。
图9示例等高线过滤示意图。
图10实施例平缓界线过滤效果图。
图11实施例岩层非平行界线过滤效果图。
图12实施例直立岩层提取结果。
具体实施方式
下面结合附图并通过描述一个自动识别并标注直立岩层的实例，来进一步说明本发明的效果。本实施例选择某地的地形地质图为示例数据(主要包括等高线图层与地质分区图层)，如图2所示。
(一)岩层边界读取和预处理
步骤1.1加载shp格式的地质体面图层数据“地质分区”，读取岩层面要素读入集合Stratums＝{stai|i＝1，2，3，...，25}，本实施例中共包含25个岩层面要素。
步骤1.2加载的地形地质图图层基本参数为：宽2316.33m，高1705.86m，图幅边界点集BoundaryP＝{(461.44,-1999.47),(507.27,-1999.54),(716.37,-1999.88),(1,646.99,-2001.37),(1749.90,-2001.53),(1,891.30,-2001.76),(2378.58,-2,002.53),(2459.91,-2002.66),(2777.77,-2003.17),(2769.89,-1095.25),(2769.89,-1095.25),(2763.88,-402.55),(2762.97,297.38),(1941.99,-298.70),(1686.37,-299.11),(1302.08,-299.73),(1058.81,-300.12),(817.65,-300.50),(700.11,-300.69),(457.74,-301.08),(459.83,1259.61),(460.39,-1514.49)}。实施例直线标准长度Llimit取值230m；霍夫变换的窗口15*15，容限Δθ取值为π/15，即0.209rad，同时此值设置为平行线筛选的阈值。平缓界线过滤的夹角阈值设置为20°，即0.349rad。
步骤1.3读取stai的边界，存为线要素oli，得到集合OriL＝{oli|i＝1，2，3，...，25}。每一个面要素stai存储为一个oli。oli继承stai的“id”属性。
步骤1.4本实施例中不存在岛(洞)的情况，所有边界线都为单个环线。得到SinL＝{sli|i＝1，2，3，...，25}，sli继承oli的“id”属性；再利用位于图幅边界的点集分割线要素，并按条件进行筛选和连接，最终得到集合FinL＝{flu|u＝1，2，3，...，28}，以ol3为例，其界线被边界分割而得到界线fl3和fl4，且fl3和fl4的“Id”属性为3。
步骤1.5全图对岩层取边界，得到预处理后的地质界线图层如图3所示。
(二)岩层平直界线提取
直线提取针对每个线要素分别进行提取，以线要素fl3为例，其点集FLP3＝{(507.27，-1999.54),(509.12，-1992.93),(518.21，-1985.494),…,(1730.64，-1981.40),(1749.90，-2001.53)}，包括38个点。
步骤2.1对点集FLP3进行霍夫变换，以π/15为间隔计算得每个点的霍夫变换曲线。以点flp31(507.27，-1999.54)为例，霍夫变换后得到曲线HoughL31<(0.507.27),(0.21，80.46),(0.42，-349.87),(0.63，-764.91),(0.84，-1146.52),…,(2.93，-911.91),(3.14，507.27)>，绘制曲线如图4所示。
步骤2.2创建一个直线点集StaightLineP，将FLP3中的第一个点flp31加入直线点集，此时StaightLineP中只包含一个要素。
步骤2.3将FLP3中的点依次添加进StaightLineP，每次往StaightLineP中添加新点时，即对StaightLineP进行操作。以线要素fl3为例，将点flp32加入集合StaightLineP，对其进行如下操作：
1)此时StaightLineP中仅包含两个点，则往集合中添加新点flp33。
2)获取点flp31、点flp32与flp33霍夫变换后的曲线HoughL31、HoughL32与HoughL33。得到ρmax＝518.21，ρmin＝-2058.43，计算得到Δρ＝171.78，确定容限窗口为(0.21，171.78)。
3)计算HoughL31与HoughL32的交点得到点p1(2.87，-1021.44)，计算HoughL2与HoughL3的交点得到点p2(2.27，-1857.81)，两交点的差值(0.60，836.37)，超出容限窗口，则将点flp33从StaightLineP中移除，计算线SL<flp31,flp32>的长度，其值为6.86m，小于Llimit，则将flp31从StaightLineP中移除，将点flp33重新加入StaightLineP。
4)此时StaightLineP中仅包含两个点，则往集合中添加新点flp34。
5)获取点flp32、点flp33与点flp34霍夫变换后的集合HoughL32、HoughL33与HoughL34。得到ρmax＝542.38，ρmin＝-2058.43，计算得到Δρ＝173.39，确定容限窗口为(0.21，171.78)。
6)计算HoughL32与HoughL33的交点得到点p1(2.27，-1857.81)，计算HoughL33与HoughL34的交点得到点p2(2.20，-1901.60)，两交点的差值(0.07，43.79)，不超过容限窗口，继续添加新点flp35，重复对StaightLineP进行判断。
7)添加点fpl36、fpl37、flp38,…,flp311，集合StaightLineP内霍夫变换后曲线的交点皆在容限窗口内，直到添加点flp312(875.36，-1751.91)，此时对于霍夫曲线{HoughL32,HoughL33,…,HoughL312}，容限窗口为(0.21，293.88)，交点集InterP＝{p1(2.27，-1857.81)，p2(2.20，-1901.60)，p3(2.27，-1887.48)，…，p9(2.11，-1944.95),p10(2.04，-1950.72)}交点分布已超出容限窗口，则将flp312从StaightLineP移除，计算线SL<flp32,flp311>的长度，其值为385.56m，大于Llimit，则将StaightLineP＝{flp32，flp33，...，flp311}另存为直线要素l5，将其加入集合StaightL，通过计算，其属性“θ”赋值为2.09，属性“ρ”赋值为-1918.73。
8)清空StaightLineP，将点flp311加入集合，继续判断直到最后一个点flp338。
步骤2.4对直线完成平直界线提取之后，对集合StaightL中邻接且近似的线进行合并。以线要素fl3为例，其可提取出2段平直界线{l5，l6}，其中l5与l6属性“θ”与“ρ”差距小于容限，且邻接，则对l5与l6进行合并，对其进行平直界线提取效果及合并效果如图5所示。
步骤2.5对集合StaightL中的元素进行属性“id”的赋值，以线要素fl3为例，其对应集合StaightL中{l5}的“id”属性赋值为3。
(三)直立岩层识别绘制
步骤3.1完成平直界线提取后，得到所有岩层的平直界线如图6所示，全图平直界线集合AllSL＝{li|i＝1，2，3，...，48}。
步骤3.2对AllSL进行倾向断层过滤，对一组相互重叠的直线要素两侧的岩层符号序列进行判断。如图7所示，线l3、l8与线l43、l10相互重叠，且l3、l8归属的岩层其地层符号序列为{“C_2”,“C_3”}，l43、l10归属的岩层其地层符号序列为{“C_2”,“C_3”}，两个序列相同，则判断线l3、l8、l43和l10组成一组断层线，将线l3、l8、l43和l10从集合AllSL中移除。有时会遇到两个序列长度不相等的情况，则以一者包含另一者为标准进行判断。经断层过滤后得到的集合AllSL＝{li|i＝1，2，3，...，39}，图层绘制如图8所示。
步骤3.3加载等高线图层，得到等高线集合Contour＝{clj|j＝1，2，3，...，78}。基于等高线对AllSL进行平缓界线过滤。以直线要素l3、l4、l5、l6、l7和l8为例(如图9所示)，其属性与判断情况如下表所示：
对AllSL中的直线要素逐一进行平缓直线过滤，最终得到集合AllSL＝{li|i＝1，2，3，...，25}，图层绘制如图10所示。
步骤3.4将属性“Id”相等的直线lx加入同一集合SSLt，最后得到集合AllSL＝SSL1∪SSL2∪SSL3，...，SSL15，分别对应属性“id”为2、4、6、8、10、11、12、13、14、15、16、19、20、24和25的岩层面要素经过断层过滤和平缓界线过滤后的平直界线。对其进行非平行成对的界线进行过滤，依次判断SSLt中是否存在成对的直线要素lx与ly，属性“θ”差值小于容限Δθ，对不满足条件的平直界进行过滤，最后得到满足条件的界线对(如图11所示)信息见下表：
步骤3.5对于满足直立岩层条件的要素，将其成对的直线要素lx与ly首尾连接绘制新的面要素，即为岩层直立的部分。分别连接5个直立岩层部分，得到直立岩层如图12所示。