CN109472072A

CN109472072A - Interaction prediction method between ephemeral stream and underground water based on simulating river

Info

Publication number: CN109472072A
Application number: CN201811277574.8A
Authority: CN
Inventors: 陆垂裕; 李天辰; 曹国亮; 孙青言; 严聆嘉; 赵勇; 王建华; 张博; 李泽鹏; 秦韬
Original assignee: China Institute of Water Resources and Hydropower Research
Current assignee: China Institute of Water Resources and Hydropower Research
Priority date: 2018-10-30
Filing date: 2018-10-30
Publication date: 2019-03-15
Anticipated expiration: 2038-10-30
Also published as: CN109472072B

Abstract

The invention discloses a kind of interaction prediction method between ephemeral stream and underground water based on simulating river, devises a kind of algorithm for automatically determining simulaed path according to the characteristic of the network of waterways, is applied to large complicated natural/ditch network.When carrying out river number be in the present invention it is random, without considering the confluence relationship and divided relation of the network of waterways, facilitate the work of user；The function that the perfect more rivers of the present invention shunt simultaneously, can simulate the quantitative shunting in a plurality of river and shunt situation in proportion, and can divide situation of the water more than or less than gross water requirement；In addition, using the present invention to the volume forecasting to interact between ephemeral stream and underground water accuracy with higher.

Description

Interaction prediction method between ephemeral stream and underground water based on simulating river

Technical field

The invention belongs to hydrographic water resource data information automatic processing technical fields, and in particular to one kind is based on river mould The design of interaction prediction method between quasi- ephemeral stream and underground water.

Background technique

Mutually converting for surface and ground water is a kind of generally existing hydrology phenomenon of nature, it is the base of water resource One of this attribute.In carrying out water resources assessment and planning management, it is necessary to very focus on the close waterpower of surface and ground water The two is coupled together and considers that solid foundation could be provided for the evaluation and planning management of water resource by connection.Surface water with It is most important with the meaning that stream -aquiferinteraction is studied in underground water system with interaction, i.e. ephemeral stream and underground The relationship of water.Ephemeral stream generally refer to dry season river water stop riverbed it is exposed and in abundance of water season river water Pentium Ephermeral stream.

For the relationship for portraying subsurface flow and ephemeral stream, needs to carry out Underground river system sunykatuib analysis, most start, ground What lower Fluid Dynamics and prediction mostly used is simple hydrologic(al) budget method and hydrogeologic comparison, for boundary condition understanding usually River is handled as head boundary or given head boundary is determined, and this in many cases and does not meet practical situation. The water level in river, the flow exchanged with water-bearing layer are all dynamic changes in practical situation, this complex mechanism tradition side Method is difficult to portray.With the development of computer technology, numerical simulation is applied in the analysis of underground water meter is calculated, computational theory It has all made significant headway with actual measurement underground water analysis ability, the dynamic relationship of underground water and surface water has also obtained preferably It portrays.Underground river system method for numerical simulation mainly has finite difference calculus, Finite Element, boundary element method and limited bulk at present Method etc., based on these analogy methods, various groundwater Numerical Simulation softwares are continued to bring out in decades, including Visual The popular softwares such as MODFLOW, FEFlow, GMS.At present using MODFLOW as the finite difference calculus based on unit center of representative Numerical simulator solves the partial differential equation Definite problem of the description various states of underground water with the method for iterative solution.

However, research area has been carried out time and spatially by numerical simulation software based on center cell finite difference method It is discrete, replace difference quotient with difference coefficient approximation on discrete grid, the differential equation and its definite condition turned into algebraic equation, repeatedly In generation, solves and obtains head value.And item is output and input since each discrete grid cell contains, outflow of upstream etc. Inbound traffics in downstream, it is therefore necessary to be ranked up grid cell according to the confluence sequence in river.Most simulation softward is at this Aspect requires user voluntarily to specify the confluence relationship in river: serialization is carried out to all rivers, by each river in river system by confluence Sequential arrangement the river in sequence is named one by one at sequence, and with serial number, in order to carry out simulation trial later.But This intricate operation and fallibility in practical criss-cross Complex River, if network of waterways relationship or confluence sequence misarrangement will be led The analog result of entire program is caused to occur calculating mistake, when program is without inspection, the hiding error of generation, it is difficult to check.

In addition to this, existing simulation softward works as river when generating the network of waterways to the artificial processing for shunting river and not flexible Stream needs when being shunted according to split ratio, and user is needed to be input to each shunting from being about to split ratio and be converted into firm discharge value , and can not handle can point water be less than all rivers Fen Shui in downstream gross water requirement situation.

Summary of the invention

The purpose of the present invention is to solve the existing simulation softwards of use to phase interaction between ephemeral stream and underground water The above problem existing for numerical prediction proposes a kind of based on the ephemeral stream of simulating river and phase interaction between underground water Use prediction technique.

The technical solution of the present invention is as follows: Interaction Predicting side between ephemeral stream and underground water based on simulating river Method, comprising the following steps:

S1, the river network system of ephemeral stream is spatially divided into the network of waterways, river system, river and the level of section four.

S2, river and section are numbered at random.

S3, the data information that every river includes and the data information that each section includes are obtained, and River Data is believed Breath is with section data information memory into the different field of database.

S4, stream structure body is constructed according to River Data information and section data information.

S5, concentration of river network structure is constructed according to all stream structure bodies.

S6, the simulation order in each river in concentration of river network structure is identified, generates simulaed path.

Flow between S7, the calculating simulation path river Zhong Meitiao and underground water.

S8, the flow for judging whether between the river and underground water of all stress phases have all been calculated and have been finished, if then entering Step S9, otherwise return step S4.

S9, calculated result is saved into database, is completed pre- to what is interacted between ephemeral stream and underground water It surveys.

Further, the River Data information in step S3 includes:

The number in river is stored in field SEGMID since No. 1.

Downstream number indicates that the downstream river number in the river, the value are -1 no downstream river of expression, is stored in field In NEXTUNIT.

Head end inbound traffics are stored in field INFLOW.

Shunting conditions, the value are not -1 river number for indicating that this river to be shunted, for -1 indicate this river not from Other river shunts, and is stored in field DIVSEGM.

Shunting mode, if this river is shunted from other river, which is that target flow shunting is pressed in 0 expression, indicates to press for 1 and divide Water is stored in field DIVOPT than shunting.

Divide water ratio, indicate that this river shunts shared ratio, codomain is 0~1, is stored in field DIVRATIO.

Water level is calculated, which is 1 this river of expression according to Manning formula calculating water level, it is specified for 0 expression water level by user, It is stored in field BCALSTAGE.

Section data information in step S3 includes:

Stress phase number, is stored in field IPERIOD.

Grid number indicates floor, the row, column number of cell where this section, be stored respectively in field ILYR, IROW and In ICOL.

Affiliated river number, is stored in field SEGMID.

It is numbered in river, indicates that this section in the number in river, is stored in field CELLID.

Water level is stored in field STRSTAGE.

Hydraulic conductivity indicates the hydraulic conductivity between this section and underground reservoir, is stored in field COND.

Riverbed cope level is stored in field STRTOP.

Riverbed Bottom Altitude is stored in field STRBTM.

Riverbed width is stored in field STRWDT.

Gradient of river bed is stored in field STRSLP.

Manning roughness coefficient is stored in field STRNDC.

Further, when the value of field DIVOPT is 1, the shunting mode in river be according to divide water ratio from be split river into Row shunts, and shunts formula are as follows:

Wherein n is the quantity for shunting river, and i indicates i-th shunting river and i=1, and 2 ..., n, QRIV expression are split River End flow amount, QDIV_iIndicate the flow that i-th distributary river flow point obtains, α_iThe system setting point in river is shunted for i-th Water ratio, β_iFor i-th shunt river actually divide water ratio, QNRIV indicates the remaining stream for flowing into downstream river after shunting Amount.

If each item shunt river divide the sum of water ratio to be less than or equal to 1, stream flow will be gone out according to dividing water ratio to carry out It shunts, residual flow imports downstream river；If the sum of water ratio is divided to be greater than 1, lower reaches of river does not have flow, and each river that shunts is pressed According to dividing water proportional diverting.

When the value of field DIVOPT is 0, the shunting mode in river is to be shunted according to target flow from river is split, Shunt formula are as follows:

Wherein Q_iIt indicates the target flow that i-th river is got, is located in field INFLOW；If be split river goes out stream Amount is greater than or equal to each target flow for shunting river, then is shunted according to target flow, and residual flow imports downstream river； If the outflow for being split river is less than each target flow for shunting river, lower reaches of river does not have flow, and each river that shunts is pressed It is shunted according to the ratio of target flow.

Further, the stream structure body constructed in step S4 includes:

Upstream units array of pointers, for save oriented this river unit confluence upstream river unit pointer.

Dividing cell array of pointers, for saving the pointer for the river unit that this river unit to be shunted.

Downstream units pointer, the pointer of the downstream river unit for saving this river unit.

River cell data array, including river element number data, river cell attribute data and river unit point Flow data.

Section array of data, for saving the data information of all sections in this river unit.

Further, step S5 include it is following step by step:

S51, first stream structure body is selected according to the number in river.

S52, judge whether the stream structure body has downstream river structural body, if then entering step S53, otherwise enter step Rapid S54.

S53, the upstream units array of pointers that the downstream units pointer of the stream structure body is directed toward to downstream river structural body, Enter step S55.

S54, it marks the stream structure body for top river outlet, and is added into top river outlet array, into step Rapid S55.

S55, judge whether that all stream structure bodies are complete pointer relationship setting, if then entering step S56, otherwise Select next stream structure body, return step S52.

S56, first stream structure body is selected according to the number in river.

S57, judge whether the stream structure body is to shunt stream structure body, if then entering step S58, otherwise enter step Rapid S59.

S58, the pointer in the upstream units array of pointers of the shunting stream structure body is directed toward upstream stream structure body Dividing cell array of pointers, enters step S59.

S59, judge whether that all stream structure bodies are complete divided relation setting, if then completing concentration of river network structure Building, enters step S6, otherwise selects next stream structure body, return step S57.

Further, the process identified in step S6 to the simulation order in each river includes river system scanning and river system Traversal.

River system scanning the following steps are included:

A1, Initialize installation scanning times iSumScan=0.

A2, Initialize installation scanning terminate mark iFinish=0.

River system ordinal number i=1 in A3, setting concentration of river network structure.

A4, it calls recursive function RoutTest (i) to traverse i-th of river system, and obtains function RoutTest's (i) Return value；If the return value of function RoutTest (i) is 0, indicate that the simulation order in every river in this river system is all determined, If the return value of function RoutTest (i) is -1, indicate there is river not can determine that simulation in the scanning of this river system in this river system Sequentially.

A5, the return value of this river system function RoutTest (i) is added in scanning termination mark iFinish.

A6, river system ordinal number i is enabled to add 1.

A7, judge the river system number whether i is greater than in concentration of river network structure, if then entering step A8, otherwise return to step Rapid A4.

A8, judge that scanning terminates whether mark iFinish is equal to 0, if then indicating the simulation in all rivers in this river system Sequence is determined, and is terminated river system scanning process, is otherwise entered step A9.

A9, scanning times iSumScan is enabled to add 1.

A10, judge whether scanning times iSumScan is greater than unit total number in river in concentration of river network structure, if then table Show that the physical structure of the network of waterways is unreasonable, abnormal, end river system scanning process occurs in simulation, otherwise return step A2.

River system traversal the following steps are included:

The stream structure body ordinal number that B1, setting currently traverse is j.

B2, judge whether j-th of stream structure body has determined simulation order in river system scanning process before, if The return value that recursive function RoutTest (j) is then arranged is 0, enters step B13, otherwise enters step B3.

B3, judge whether the upstream units array of pointers of j-th of stream structure body has value, if then indicating j-th of river There is tributary stream structure body in structural body upstream, enters step B4, otherwise indicates that j-th of stream structure body upstream does not have tributary river Structural body enters step B11.

B4, Initialize installation end value iResult=0.

B5, Initialize installation upstream branching flow stream structure body ordinal number k=1.

B6, recursive function RoutTest (k) is called to traverse k-th of upstream branching flow stream structure body.

B7, the return value for obtaining function RoutTest (k) simultaneously accumulate it in end value iResult；If function The return value of RoutTest (k) is 0, indicates that the simulation order of the upstream branching flow stream structure body is all determined, if function The return value of RoutTest (k) is -1, indicates that the upstream branching flow stream structure body not can determine that simulation order.

B8, upstream branching flow stream structure body ordinal number k is enabled to add 1.

B9, judge whether k is greater than upstream branching flow stream structure body number, if then entering step B10, otherwise return step B6。

Whether B10, judging result value iResult are equal to 0, if then expression has determined the mould of j-th of stream structure body Quasi- sequence, the return value of setting RoutTest (j) are 0, enter step B13, otherwise indicate not determining j-th of stream structure body Simulation order, setting RoutTest (j) return value be -1, enter step B13.

B11, judge whether the dividing cell array of pointers of j-th of stream structure body has value, if then indicating j-th of river Structural body is shunted from other stream structure bodies, is entered step B12, is otherwise indicated j-th of stream structure body not from other river knots Structure body shunts, and is the upstream initial stage river of this river system, and then determine the simulation order of j-th of stream structure body, setting The return value of RoutTest (j) is 0, enters step B13.

B12, judge j-th of stream structure body shunt stream structure body in river system ergodic process before whether Determine simulation order, if then determining the simulation order of j-th of stream structure body, the return value of setting RoutTest (j) is 0, Enter step B13；Otherwise the simulation order of j-th of stream structure body can not be determined, the return value of setting RoutTest (j) is- 1, enter step B13.

B13, the return value for exporting RoutTest (j) terminate the traversal of j-th of stream structure body.

Further, the flow rate calculation formula in step S7 between river and underground water are as follows:

The wherein flow between QRIV ' expression river and underground water is stored in field STREAM, and value is timing expression The leakage of river direction underground reservoir indicates that underground reservoir flows to the excretion in river when value is negative；CRIV is indicated Hydraulic conductivity between river-underground reservoir is stored in field COND；h_i,j,kIndicate the unit meter where this section Head is calculated, RBOT indicates the elevation of certain point of riverbed bases, is stored in field STRBTM；HRIV indicates the water level of section, It is stored in field STRSTAGE.

Further, the determination method of the water level HRIV of section are as follows:

If the value of field BCALSTAGE is 0, the water level HRIV of section is specified by user；If the value of field BCALSTAGE It is 1, then section water level, calculation formula is calculated according to Manning formula are as follows:

Wherein Q is the inbound traffics of each section in river, and n is Manning roughness coefficient, is stored in field STRNDC, and c is river Hydraulic conductivity between section and underground reservoir, is stored in field COND；W is the riverbed width of section, is stored in field In STRWDT；S is the gradient of river bed of section, is stored in field STRSLP.

The beneficial effects of the present invention are:

(1) when carrying out river number be in the present invention it is random, without considering the confluence relationship and divided relation of the network of waterways, The work of user is facilitated, user is after the completion of all river numbers to research area, it is only necessary to further point out every river stream The downstream river course number entered, it is subsequent can be by searching in the network of waterways each simulating river sequence algorithm automatically the river in the entire network of waterways Stream is finished according to confluence order division and is sequentially named one by one, to meet the needs of quantitative analysis.

(2) function that the perfect more rivers of the present invention shunt, can simulate the quantitative shunting in a plurality of river and divide in proportion Situation is flowed, and the situation that water can be divided to be more than or less than gross water requirement.

(3) with higher to the volume forecasting to interact between ephemeral stream and underground water accurate using the present invention Property, compared with actual conditions, relative error is only 0.005%.

Detailed description of the invention

Fig. 1 interacts between showing the ephemeral stream and underground water provided in an embodiment of the present invention based on simulating river Prediction technique flow chart.

Fig. 2 show the relation schematic diagram of underground water difference unit provided in an embodiment of the present invention and river, section.

Fig. 3 show stream structure body schematic internal view provided in an embodiment of the present invention.

Fig. 4 show concentration of river network structural schematic diagram provided in an embodiment of the present invention.

Fig. 5 show building concentration of river network structure flow chart provided in an embodiment of the present invention.

Fig. 6 show river system scanning process figure provided in an embodiment of the present invention.

Fig. 7 show river system traversal flow chart provided in an embodiment of the present invention.

Fig. 8 show level of ground water provided in an embodiment of the present invention and drops to schematic diagram under riverbed bottomset.

Fig. 9 show the schematic bottom view that level of ground water provided in an embodiment of the present invention is higher than riverbed bottomset.

Figure 10 show river networks structure and tapping condition schematic diagram provided in an embodiment of the present invention.

Figure 11 show concentration of river network relation schematic diagram provided in an embodiment of the present invention.

Figure 12 show unit outflow schematic diagram in the network of waterways provided in an embodiment of the present invention.

Figure 13 show unit leakage schematic diagram in the network of waterways provided in an embodiment of the present invention.

Specific embodiment

Carry out detailed description of the present invention illustrative embodiments with reference to the drawings.It should be appreciated that shown in attached drawing and The embodiment of description is only exemplary, it is intended that is illustrated the principle and spirit of the invention, and is not limited model of the invention It encloses.

The embodiment of the invention provides a kind of Interaction Predicting between ephemeral stream and underground water based on simulating river Method, as shown in Figure 1, including the following steps S1~S9:

The network of waterways refers to the surface water channel system that all rivers staggeredly in length and breadth are constituted within the scope of research area.In the network of waterways River may export and flow out through multiple basins, therefore the network of waterways is divided into multiple river systems, the corresponding basin outlet of each river system. River system is made of a plurality of river, and river therein is a certain channel of insulated stream warp, the concept with upstream and downstream, water Amount flow to downstream via the upstream in every river, and imports its next stage river, becomes the inbound traffics of next stage upper river, Water passes through the confluence step by step in river, finally flows out from the lower exit (also referred to as basin outlet) in the most junior river of river system.

In the embodiment of the present invention, each river is independent from each other surface water channel unit, in order to simulate river surface water The effect of intercoupling between underground water is needed river river by the spatial relationship between its underground water difference gridding Stream is divided into each section.Section refers to one section that certain river is distributed in certain underground water grid cell, is simulating riverway earth's surface The base unit of water and underground water interaction.Water is flowed from upstream toward downstream in river, therefore river is from upstream to down Trip is divided paragraph by paragraph by underground water difference gridding unit and section is numbered in order, and from small to large relatively suitable is numbered in section Sequence represents the flow direction of water flow.Fig. 2 shows river, underground water difference gridding unit, the schematic relationships between section, figure There are 7 rivers in 2, wherein the 1st article of river is divided into 2 sections, the 2nd article of river is divided into 4 sections.

S2, river and section are numbered at random.

All there is unique number to be used to distinguish different rivers in every river in research area, and the number in river is from 1 Start serial number, therefore maximum river number is to study the quantity in area river.In the embodiment of the present invention, the number in river There is no any other meaning other than identifying the river, therefore the number for studying the river area Zhong Getiao can be random disorder , there is very high number freedom degree.Existing other technologies method usually requires that the number in river while also indicating that river exists Tandem in simulation process, such as MODFLOW, it means that user must first identify the confluence in the network of waterways between river Serial number is carried out to river again after relationship and divided relation, and suitable before and after representing river in simulation process with number order Sequence, this is not only very heavy work load, and error-prone for the more complex network of waterways.In the present invention user into Without considering the confluence relationship and divided relation of the network of waterways when row river number, this facilitate that the work of user.User is to research After the completion of all rivers number in area, it is only necessary to the downstream river course number that every river flows into is further pointed out, it can in subsequent step To establish concentration of river network structure automatically by respective algorithms, and identify the simulation order in the network of waterways river Zhong Getiao.It is illustrated in Fig. 2 Study the number in 7 rivers in area.

Every river is divided into different sections by its distribution in underground water difference gridding, and Fig. 2 illustrates section Division and coding rule.Cell (1,3) indicates the cell for being located at the 1st row the 3rd column.There are two sections in river 1, along river The direction of flowing starts from unit (1,3), terminates to unit (2,3).Some of flows are introduced in river 2, remaining Flow enters river 3 from river 1, and river 3 has 4 sections.River 2 and river 4 cross in unit (5,3), to be formed River 5 comprising 2 sections, river 3, river 5 and river 6 cross in unit (5,4), to form river 7.River 3 and river Two fraction rivers in stream 6 are not included in numbering plan, because they are just through cell (2,4) and (3,5) Corner.In addition, multiple sections in different rivers can be assigned to identical model unit lattice, such as second of the 5th article of river Section；4th section in the 3rd article of river；5th section in the 6th article of river；1st section in the 7th article of river is all located at unit (5,4) in.

Wherein, the data information that river includes is as shown in table 1, and different data information is stored in the different field of database In.

Table 1

In the embodiment of the present invention, when the value of field DIVOPT is 1, the shunting mode in river be according to divide water ratio from being divided Stream river is shunted, and formula is shunted are as follows:

Wherein n is the quantity (n >=1) for shunting river, and i indicates that i-th shunting river and i=1,2 ..., n, QRIV are indicated It is split river End flow amount, QDIV_iIndicate the flow that i-th distributary river flow point obtains, α_iThe system for shunting river for i-th Setting divides water ratio, β_iFor i-th shunt river actually divide water ratio, QNRIV indicates remaining after shunting and flows into downstream river The flow of stream.In the embodiment of the present invention, since each shunting river of system setting divides water ratio α_iThe sum of be likely larger than 1, because This introducing actually divides water ratio beta_iThe actual point of water ratio in i-th shunting river when as simulation.

When the value of field DIVOPT is 0, the shunting mode in river is to be divided according to target flow from river is split Stream shunts formula are as follows:

The data information that section includes is as shown in table 2, these information datas are simulation river surface and ground water waterpower The basis of connection.

Table 2

On the basis of river, section data information are complete, the data structure of river unit can be designed, thus in mould Facilitate tissue and the use of information data during quasi-, the internal structure of stream structure body is as shown in figure 3, each stream structure body In include:

On the basis of stream structure body, the remittance of research area's entirety network of waterways can be constructed using stream structure body pointer to object Flow structure.Stream structure body object will be generated for the data information in each river and its section, recycles the upper and lower of river Trip relationship, river divided relation and stream structure body pointer to object connect different stream structure bodies, are formed in simulation Concentration of river network structure, as shown in Figure 4.

In concentration of river network structure building process, need to consider natural and artificial both sides factor, first is that in the network of waterways not With the upstream and downstream relationship in river, i.e., the topological structure of the natural network of waterways；Second is that the artificial shunting connection between different rivers, i.e. people Water transfer relationship under class intervention.Under the natural situation of no manual intervention, river is always downstream flowed by the primary river of upstream Domain outlet confluence.As shown in figure 5, the building process of concentration of river network structure includes the following steps S51~S59:

S51, first stream structure body is selected according to the number in river.

Before step S51, it is necessary first to according to all words in parsing database about ephemeral stream data boundary Section therefrom extracts river information and attribute information, and determines river number.Check that number continuity uses Hash table structure later The number information for storing river, is inquired when subsequent analysis being facilitated to calculate.

S56, first stream structure body is selected according to the number in river.

Judge whether some stream structure body is to shunt stream structure body method particularly includes: when certain stream structure body unit Head end in data becomes a mandarin magnitude when not being 0 and shunting river number for -1, which is to shunt stream structure Body.Wherein, head end inbound traffics and shunting river number parse to obtain from INFLOW the and DIVSEGM field in database respectively.

The automatic identification of each simulating river sequence is core of the invention in concentration of river network structure.The embodiment of the present invention In, a simulation order array is set first, as soon as the simulation order in river is often identified, by the object of the stream structure body Pointer is placed in the end of simulation order array.In this way after the simulation order in all rivers identifies, the simulation order number The sequencing of each element just represents the simulation order in each river in group.

The algorithm flow of identification simulating river sequence includes the steps that two big: river system scanning is traversed with river system.If river Nobody's work point water contacts between each river in net, and the simulation order in each river can be directly by carrying out each river system The river traversal for being from upstream to downstream determines；If existing between each river in the network of waterways manually divides water to contact, need to lead to It crosses river system traversal and river system scans the method combined and determines simulaed path.

Due to dividing water relationship in the presence of complicated in concentration of river network structure, if the river shunted exists to Mr. Yu river It is simulated before the simulating river not yet, then this river cannot be determined in simulation is got from the river of shunting The process of flow, shunting cannot be realized.It, can without rule since divided relation is totally independent of relationship of naturally converging step by step Speech, it is therefore desirable to just can determine that final reasonable simulating river is suitable by being constantly scanned test to river system each in the network of waterways Sequence.

Recursive function RoutTest is called one by one to river system all in concentration of river network structure, to complete all river systems certainly The traversal of body, referred to as single pass.There is no when divided relation, the confluence sequence of each river system is all between river in the network of waterways It is independent, therefore 1 scanning is only needed, the simulation order in each river can determine in river network system.However it is deposited when between river When dividing water relationship, certain river in some river system may need to wait certain simulating river in other river system it is complete after It can determine its upland water amount, cause the single pass process of the river system that can't complete all simulating rivers in the river system suitable The sequence of sequence.In the embodiment of the present invention, the mould for being repeatedly scanned with the algorithm of the network of waterways to ensure every river in each river system is proposed Quasi- sequence can be determined reasonably, as shown in fig. 6, river system scanning includes the following steps A1~A10:

A1, Initialize installation scanning times iSumScan=0.

A2, Initialize installation scanning terminate mark iFinish=0.

A6, river system ordinal number i is enabled to add 1.

A9, scanning times iSumScan is enabled to add 1.

Every time when updating scanning times iSumScan, need to check it.Since the river unit in the network of waterways is total Number is fixed, and scans the simulation order that at least can more determine a river every time, therefore scanning times iSumScan exists One threshold value.If the value of scanning times iSumScan is greater than the total number of river unit in the network of waterways, it is abnormal to show that simulation occurs, The physical structure of the network of waterways is unreasonable, exits simulation and terminates river system scanning process.

Each river system has a top river outlet (the river unit of most downstream) in the network of waterways, and the process of river system traversal is adopted With the method for recursive call, upstream stream structure body is traversed since the outlet of the basin of each river system, until finding upstream most Primary river (not from the river that other stream structure body shunts), then downstream sorts since primary river.Such as Fig. 7 institute Show, river system traversal includes the following steps B1~B13:

B4, Initialize installation end value iResult=0.

Flow rate calculation formula between river and underground water are as follows:

It calculates flow QRIV ' and needs to use three parameters: the water between the water level HRIV of section, river-underground reservoir The elevation location RBOT of power coefficient of conductivity CRIV and riverbed bases, and need to tear open the parameter in QRIV ' expression formula (3) Divide and substitutes into the difference equation based on water quantity exchange.

Work as h_i,j,kWhen≤RBOT, as shown in figure 8, h_i,j,kIt has fallen under riverbed bottomset bottom surface, in its lower formation One unsaturation band.If it is assumed that riverbed bottomset itself keeps being saturated, the head of certain point of riverbed bases can be simply It is considered as the elevation of the point, i.e. RBOT, then the flow between river and underground water is a constant value.Calculate when need by- CRIV (HRIV-RBOT) is added in Difference Terms RHS；Work as h_i,j,kWhen > RBOT, as shown in figure 9, h_i,j,kHigher than riverbed bottomset Bottom surface, the flow between river and underground water is directlyed proportional to the head difference of river and underground reservoir, with head h_i,j,kAs Variable.It needs for-CRIV to be added in Difference Terms HCOF when calculating, and-CRIV*HRIV is added in Difference Terms RHS, from And QRIV ' fractionation is substituted into the difference equation based on water quantity exchange.

The difference equation simultaneous that all computing units are obtained, and solved using the method for iteration.In the mistake of iteration Cheng Zhong, the result of each iteration will all be used for calculating next time after treatment.Different algorithms has different processing methods, Under normal circumstances, the change of water level after each iteration is gradually reduced, and is finally reached convergence.This completes a periods Head calculate.Whether restrain, is usually determined by a convergence index predetermined, when the most flood iterated to calculate twice When head difference is less than the convergence index, referred to as restrain.Since initial head, at the end of each step finds out each period Head value, and the value is used and such process to be repeated, until required time knot as the initial value of next period Beam can thus obtain the flow in simulaed path between every river and underground water.

The determination method of the water level HRIV of section are as follows:

In the embodiment of the present invention, either specified water level, which is calculated, carries out automatic water level calculating according further to Manning formula, and one It, will be as this after the upstream inbound traffics of a section subtract the leakage (or seepage discharge plus underground water to this section) of this section The downstream outflow of section flows to next section, becomes the upstream inbound traffics of next section.It so gradually can be by river stream It measures from the section of the most upstream in river and calculates most downstream section to river, the flow and the upper of next river will be become Inbound traffics are swum, which will be continued until that river water leaves river network system.

Effect of the invention is described further with a specific embodiment below:

As shown in Figure 10, there are 9 rivers in the network of waterways, the number in river is not carried out according to the sequence of concentration of river network, each river Manning roughness coefficient is 0.013.

Wherein, there are 6 sections in river 1, length be respectively 2333.21m, 5203.9m, 3093.2m, 5303.9m, 5589.5m,2748.8m；Riverbed cope level is -0.5m；Riverbed Bottom Altitude is -1.1m；The width in riverbed is 12m；The gradient in riverbed is 0.0012.

There are 7 sections in river 2, length be respectively 4302.4m, 5003.7m, 5003.7m, 5003.7m, 5003.7m, 5003.7m,3279.8m；Riverbed cope level is -2m；Riverbed Bottom Altitude is -3m；The width in riverbed is 20m；River The gradient of bed is 0.001.

There are 10 sections in river 3, length be respectively 44615.4m, 6349.3m, 6350.4m, 6350.4m, 6350.4m, 6350.4m,6350.4m,6350.4m,6350.4m,3657.2m；Riverbed cope level is -2.5m；Riverbed Bottom Altitude is equal For -3.3m；The width in riverbed is 20m；The gradient in riverbed is 0.001.

There are 4 sections in river 4, and length is respectively 5584.0m, 2750.6m, 4756.1m, 6507.3m；Riverbed cope level It is -2.6m；Riverbed Bottom Altitude is -3.4m；The width in riverbed is 20m；The gradient in riverbed is 0.001.

There are 9 sections in river 5, length be respectively 2872.5m, 4636.0m, 3969.2m, 5190.0m, 5190.0m, 3572.2m,4983.6m,5190.0m,5190.0m；Riverbed cope level is -3m；Riverbed Bottom Altitude is -4m；Riverbed Width be 20m；The gradient in riverbed is 0.001.

There are 4 sections in river 6, and length is respectively 2647.4m, 6812.4m, 8235.4m, 4735.2m；Riverbed cope level It is -2.3m；Riverbed Bottom Altitude is -3.1m；The width in riverbed is 20m；The gradient in riverbed is 0.001.

There are 8 sections in river 7, length be respectively 5463.0m, 8445.2m, 8445.2m, 8445.2m, 8445.2m, 3686.0m,8445.2m,4759.2m；Riverbed cope level is -3.2m；Riverbed Bottom Altitude is -3.8m；The width in riverbed Degree is 20m；The gradient in riverbed is 0.001.

There are 4 sections in river 8, and length is respectively 7161.3m, 13369.2m, 16279.3m, 6149.1m；It is high at the top of riverbed Cheng Junwei -5m；Riverbed Bottom Altitude is -9m；The width in riverbed is 20m；The gradient in riverbed is 0.001.

There are 10 sections in river 9, length be respectively 2781.2m, 3001.0m, 2321.0m, 2543.0m, 2442.0m, 2797.0m,2611.0m,2351.0m,2124.0m,2100.0m；Riverbed cope level is -1m；Riverbed Bottom Altitude is- 2m；The width in riverbed is 12m；The gradient in riverbed is 0.0012.

Unit grid is 8m/s along the infiltration coefficient of line direction；Vertical hydraulic conductivity is 0.00002；At the beginning of grid cell Beginning head is -1m.The initial flow in river 1 is 12000m³/ min, the initial flow in river 8 are 10000m³/ min, river 3 from River 7 divides water to be 4000m³/min.River 5, river 9,2 three, river river divide the ratio of water to be respectively 3:6 from river 1: 1.The mikey of simulation is m；Chronomere is day.

Concentration of river network structure is constructed using step S2~S5 first, the command keyword in searching database, including The fields such as SEGMID, NEXTUNIT, DIVSEGM.As shown in table 3, the downstream river number in river 1 is river 2；Under river 2 Swimming river number is river 4；The downstream river number in river 4 is -1, i.e., river 4 does not have downstream river, is top river outlet. So each stream structure body is analyzed, establishes concentration of river network structure.

Table 3

Water transfer relationship is arranged on the basis of the natural network of waterways: river 3 is shunted from river 7；River 5 divides with river 9 from river 1 Stream；DIVSEGM field is that -1 corresponding river does not shunt.It finds all associated river units and pointer relationship is set, shape At concentration of river network relationship it is as shown in figure 11.

Then pass through the river system traversal and the automatic search simulaed path of river system scanning in step S6.In concentration of river network structure altogether There are two river systems, are left river system (including river 1,2,3,4) and right river system (river 5,6,7,8,9).Gone out by the top basin of river system Upstream river traverses mouth.

The scanning of first time river system is carried out, first left river system is traversed.Since river system outlet, river 4 is traversed, discovery has Upstream branching flow river 2,3 first carries out recursive call to river 2.River 2 is traversed, discovery has upstream river 1, passs river 1 Return calling.River 1 is traversed, river 1 is left river system upstream initial stage river, does not have upstream river, is added into simulated series number Group and return value 0 terminate the recursive call in river 1, back to the recursive call function in river 2.River 2 obtains return value 0, Upstream has determined that simulation order, and simulated series array and return value 0 is added in river 2.Terminate the recursive call in river 2, returns To the recursive call function in river 4.River 4 obtains first upstream branching flow return value 0, then traverses Article 2 upstream river 3. River 3 is shunted from river 7 and river 7 does not determine simulation order, is stopped left river system and is traversed and return to -1.On 4 second, river Swimming return value is -1, acquires the sum of two return values and is -1 and returns to left river system.Then right river system is traversed.From river system Outlet starts, and traverses river 6, and discovery has upstream branching flow river 9,7, first carries out recursive call to river 9.Traverse river 9, river 9 from the shunting of river 1 and river 1 has determined that simulation order, and simulated series array and return value 0 is added in river 9.Terminate river 9 Recursive call, back to the recursive call function in river 6.River 6 obtains first upstream return value 0, then traverses second Upstream river 7.River 7 is traversed, discovery has upstream branching flow river 5,8, first carries out recursive call to river 5.River 5 is traversed, River 5 is shunted from river 1 and river 1 has determined that simulation order, and simulated series array and return value 0 is added in river 5.Terminate river The recursive call of stream 5, back to the recursive call function in river 7.River 7 obtains first upstream return value 0, then in traversal Swim river 8.River 8 is traversed, river 8 is the top river in upstream, is added into simulated series array and return value 0.River 7 obtains Second upstream return value 0, the sum of all upstream return values are that 0 i.e. all upstream rivers are determined simulation order, by river Simulated series array and return value 0 is added in stream 7.The recursive call for terminating river 7, back to the recursive call function in river 6.River Stream 6 obtains second upstream return value 0, and all upstream rivers are determined simulation order, and simulated series number is added in river 6 Group simultaneously returns to right river system for 0, and all rivers of right river system have determined that simulation order at this time, completes the scanning of first time river system.

The sum of two river system return values are -1, carry out second of river system scanning.Left river system is traversed, river 4 is traversed, Its upstream river 2 has determined that simulation order return value 0, terminates the recursive call in river 2, back to the recursive call letter in river 4 Number.Then Article 2 upstream river 3 is traversed, river 3 is shunted from river 7 and river 7 has determined that simulation order, and river 3 is added Simulated series array and return value 0.River 4 obtains second upstream return value 0, and the sum of all upstream return values are 0, will River 4 is added simulated series array and returns to left river system for 0.Right river system return value is 0, then does not traverse to right river system, Complete second of river system scanning.

The sum of two river system return values are 0, it has been determined that the simulation order in all rivers in the network of waterways stops scanning.

The flow between river and underground water is predicted by step S7~S9, and prediction result is analyzed, Available technical effect of the invention are as follows:

(1) ability of automatic search simulaed path is embodied.By the analysis to the network of waterways, river in the entire network of waterways according to Confluence sequence is numbered one by one, and the results are shown in Table 4, and wherein IORDER field is simulated series, and UNITID field is that river is former Begin number.Therefore when carrying out river original number be it is random, without considering the confluence relationship and divided relation of the network of waterways, facilitate The work of user.

Table 4

IPERIOD	IORDER	UNITID
			0	1	1
0	2	2
			0	3	9
0	4	5
			0	5	8
0	6	7
			0	7	3
0	8	6
			0	9	4

(2) function that a plurality of river is shunted by quantitative shunting and in proportion is realized.River 5, river 9,2 three, river river Stream divides the flow of water to be respectively as follows: 1979.95m from river 1³/min、3959.89m³/min、659.99m³/min.Ratio is respectively 3:6:1；7 outflow of river is 12916.03m³/ min divides to walk 4000m through river 3³8915.73m after the flow of/min³/ min's Water imports downstream river 6.

(3) accuracy of simulation and forecast result is sufficiently illustrated.The outflow and leakage of network of waterways each unit grid are as schemed 12, shown in 13, the absolute water balance error of the whole district is -0.822m³；Relative water content balance error is 0.005%.Wherein, season The leakage of section property covalency ground water regime is 14303.2m³；Ground water regime is to the excretion of ephemeral stream 27.9m³。

Those of ordinary skill in the art will understand that the embodiments described herein, which is to help reader, understands this hair Bright principle, it should be understood that protection scope of the present invention is not limited to such specific embodiments and embodiments.This field Those of ordinary skill disclosed the technical disclosures can make according to the present invention and various not depart from the other each of essence of the invention The specific variations and combinations of kind, these variations and combinations are still within the scope of the present invention.

Claims

1. interaction prediction method between ephemeral stream and underground water based on simulating river, which is characterized in that including following Step:

S1, the river network system of ephemeral stream is spatially divided into the network of waterways, river system, river and the level of section four；

S2, river and section are numbered at random；

S3, obtain every river data information for including and data information that each section includes, and by River Data information with Section data information memory is into the different field of database；

S4, stream structure body is constructed according to River Data information and section data information；

S5, concentration of river network structure is constructed according to all stream structure bodies；

S6, the simulation order in each river in concentration of river network structure is identified, generates simulaed path；

Flow between S7, the calculating simulation path river Zhong Meitiao and underground water；

S8, the flow for judging whether between the river and underground water of all stress phases have all been calculated and have been finished, if then entering step S9, otherwise return step S4；

S9, calculated result is saved into database, completes the prediction to interacting between ephemeral stream and underground water.

2. interaction prediction method between ephemeral stream according to claim 1 and underground water, which is characterized in that described River Data information in step S3 includes:

The number in river is stored in field SEGMID since No. 1；

Downstream number indicates that the downstream river number in the river, the value are -1 no downstream river of expression, is stored in field In NEXTUNIT；

Head end inbound traffics are stored in field INFLOW；

Shunting conditions, which, which is not -1, indicates this river river to be shunted number, indicates this river not from other for -1 River shunts, and is stored in field DIVSEGM；

Shunting mode, if this river is shunted from other river, which is 0 to indicate to shunt by target flow, is indicated for 1 by dividing water ratio It shunts, is stored in field DIVOPT；

Divide water ratio, indicate that this river shunts shared ratio, codomain is 0~1, is stored in field DIVRATIO；

Water level is calculated, which is 1 this river of expression according to Manning formula calculating water level, is specified, is stored by user for 0 expression water level In field BCALSTAGE；

Section data information in the step S3 includes:

Stress phase number, is stored in field IPERIOD；

Grid number, the floor of cell, row, column number, are stored respectively in field ILYR, IROW and ICOL where indicating this section In；

Affiliated river number, is stored in field SEGMID；

It is numbered in river, indicates that this section in the number in river, is stored in field CELLID；

Water level is stored in field STRSTAGE；

Hydraulic conductivity indicates the hydraulic conductivity between this section and underground reservoir, is stored in field COND；

Riverbed cope level is stored in field STRTOP；

Riverbed Bottom Altitude is stored in field STRBTM；

Riverbed width is stored in field STRWDT；

Gradient of river bed is stored in field STRSLP；

Manning roughness coefficient is stored in field STRNDC.

3. interaction prediction method between ephemeral stream according to claim 2 and underground water, which is characterized in that described When the value of field DIVOPT is 1, the shunting mode in river is to shunt formula according to dividing water ratio to be shunted from being split river Are as follows:

Wherein n is the quantity for shunting river, and i indicates i-th shunting river and i=1, and 2 ..., n, QRIV expression are split river End flow amount, QDIV_iIndicate the flow that i-th distributary river flow point obtains, α_iThe system setting for shunting river for i-th divides water ratio Example, β_iFor i-th shunt river actually divide water ratio, QNRIV indicates the remaining flow for flowing into downstream river after shunting；

If each item shunt river divide the sum of water ratio to be less than or equal to 1, stream flow will be gone out according to dividing water ratio to shunt, Residual flow imports downstream river；If the sum of water ratio is divided to be greater than 1, lower reaches of river does not have flow, it is each shunt river according to point Water proportional diverting；

When the value of the field DIVOPT is 0, the shunting mode in river is to be shunted according to target flow from river is split, Shunt formula are as follows:

Wherein Q_iIt indicates the target flow that i-th river is got, is located in field INFLOW；If the outflow for being split river is big In or equal to each target flow for shunting river, then shunted according to target flow, residual flow imports downstream river；If by The outflow for shunting river is less than each target flow for shunting river, then lower reaches of river does not have flow, and each river that shunts is according to finger The ratio of constant flow is shunted.

4. interaction prediction method between ephemeral stream according to claim 2 and underground water, which is characterized in that described The stream structure body constructed in step S4 includes:

Upstream units array of pointers, for save oriented this river unit confluence upstream river unit pointer；

Dividing cell array of pointers, for saving the pointer for the river unit that this river unit to be shunted；

Downstream units pointer, the pointer of the downstream river unit for saving this river unit；

River cell data array, including river element number data, river cell attribute data and river unit porthole bridge According to；

5. interaction prediction method between ephemeral stream according to claim 4 and underground water, which is characterized in that described Step S5 include it is following step by step:

S51, first stream structure body is selected according to the number in river；

S52, judge whether the stream structure body has downstream river structural body, if then entering step S53, otherwise enter step S54；

S53, the upstream units array of pointers that the downstream units pointer of the stream structure body is directed toward to downstream river structural body, enter Step S55；

S54, it marks the stream structure body for top river outlet, and is added into top river outlet array, enter step S55；

S55, judge whether that all stream structure bodies are complete pointer relationship setting, if then entering step S56, otherwise select Next stream structure body, return step S52；

S56, first stream structure body is selected according to the number in river；

S57, judge whether the stream structure body is to shunt stream structure body otherwise to enter step if then entering step S58 S59；

S58, the shunting that the pointer in the upstream units array of pointers of the shunting stream structure body is directed toward to upstream stream structure body Unit pointer array, enters step S59；

S59, judge whether that all stream structure bodies are complete divided relation setting, if then completing concentration of river network structure structure It builds, enters step S6, otherwise select next stream structure body, return step S57.

6. interaction prediction method between ephemeral stream according to claim 5 and underground water, which is characterized in that described The process identified in step S6 to the simulation order in each river includes river system scanning and river system traversal；

River system scanning the following steps are included:

A1, Initialize installation scanning times iSumScan=0；

A2, Initialize installation scanning terminate mark iFinish=0；

River system ordinal number i=1 in A3, setting concentration of river network structure；

A4, it calls recursive function RoutTest (i) to traverse i-th of river system, and obtains the return of function RoutTest (i) Value；If the return value of function RoutTest (i) is 0, indicate that the simulation order in every river in this river system is all determined, if letter The return value of number RoutTest (i) is -1, indicates have river not can determine that simulation order in the scanning of this river system in this river system；

A5, the return value of this river system function RoutTest (i) is added in scanning termination mark iFinish；

A6, river system ordinal number i is enabled to add 1；

A7, judge the river system number whether i is greater than in concentration of river network structure, if then entering step A8, otherwise return step A4；

A8, judge that scanning terminates whether mark iFinish is equal to 0, if then indicating the simulation order in all rivers in this river system It is determined, terminates river system scanning process, otherwise enter step A9；

A9, scanning times iSumScan is enabled to add 1；

A10, judge whether scanning times iSumScan is greater than unit total number in river in concentration of river network structure, if then indicating river The physical structure of net is unreasonable, and abnormal, end river system scanning process occurs in simulation, otherwise return step A2；

River system traversal the following steps are included:

The stream structure body ordinal number that B1, setting currently traverse is j；

B2, judge whether j-th of stream structure body has determined simulation order in river system scanning process before, if then setting The return value for setting recursive function RoutTest (j) is 0, enters step B13, otherwise enters step B3；

B3, judge whether the upstream units array of pointers of j-th of stream structure body has value, if then indicating j-th of stream structure There is tributary stream structure body in body upstream, enters step B4, otherwise indicates that j-th of stream structure body upstream does not have tributary stream structure Body enters step B11；

B4, Initialize installation end value iResult=0；

B5, Initialize installation upstream branching flow stream structure body ordinal number k=1；

B6, recursive function RoutTest (k) is called to traverse k-th of upstream branching flow stream structure body；

B7, the return value for obtaining function RoutTest (k) simultaneously accumulate it in end value iResult；If function RoutTest (k) return value is 0, indicates that the simulation order of the upstream branching flow stream structure body is all determined, if function RoutTest (k) Return value be -1, indicate that the upstream branching flow stream structure body not can determine that simulation order；

B8, upstream branching flow stream structure body ordinal number k is enabled to add 1；

B9, judge whether k is greater than upstream branching flow stream structure body number, if then entering step B10, otherwise return step B6；

Whether B10, judging result value iResult are equal to 0, if then indicating to have determined that the simulation of j-th of stream structure body is suitable Sequence, the return value of setting RoutTest (j) are 0, enter step B13, otherwise indicate the mould that can not determine j-th of stream structure body Quasi- sequence, the return value of setting RoutTest (j) are -1, enter step B13；

B11, judge whether the dividing cell array of pointers of j-th of stream structure body has value, if then indicating j-th of stream structure Body is shunted from other stream structure bodies, is entered step B12, is otherwise indicated j-th of stream structure body not from other stream structure bodies It shunts, is the upstream initial stage river of this river system, and then determine the simulation order of j-th of stream structure body, setting The return value of RoutTest (j) is 0, enters step B13；

Whether B12, the stream structure body for judging that j-th of stream structure body shunts have determined in river system ergodic process before Simulation order, if then determining the simulation order of j-th of stream structure body, the return value of setting RoutTest (j) is 0, is entered Step B13；Otherwise can not determine that the simulation order of j-th of stream structure body, the return value of setting RoutTest (j) they are -1, into Enter step B13；

7. interaction prediction method between ephemeral stream according to claim 2 and underground water, which is characterized in that described Flow rate calculation formula in step S7 between river and underground water are as follows:

The wherein flow between QRIV ' expression river and underground water is stored in field STREAM, and value is that timing indicates river The leakage for flowing to underground reservoir indicates that underground reservoir flows to the excretion in river when value is negative；CRIV indicates river- Hydraulic conductivity between underground reservoir is stored in field COND；h_i,j,kUnit where indicating this section calculates water Head, RBOT indicate the elevation of certain point of riverbed bases, are stored in field STRBTM；HRIV indicates the water level of section, storage In field STRSTAGE.

8. interaction prediction method between ephemeral stream according to claim 7 and underground water, which is characterized in that described The determination method of the water level HRIV of section are as follows:

If the value of field BCALSTAGE is 0, the water level HRIV of section is specified by user；If the value of field BCALSTAGE is 1, Section water level, calculation formula are then calculated according to Manning formula are as follows:

Wherein Q be each section in river inbound traffics, n is Manning roughness coefficient, be stored in field STRNDC, c be section with Hydraulic conductivity between underground reservoir is stored in field COND；W is the riverbed width of section, is stored in field STRWDT In；S is the gradient of river bed of section, is stored in field STRSLP.