diff --git a/src/main/java/com/gunshi/project/xyt/controller/ForecastUseparamController.java b/src/main/java/com/gunshi/project/xyt/controller/ForecastUseparamController.java index 8bcdd5e..e0ebdc5 100644 --- a/src/main/java/com/gunshi/project/xyt/controller/ForecastUseparamController.java +++ b/src/main/java/com/gunshi/project/xyt/controller/ForecastUseparamController.java @@ -1,12 +1,16 @@ package com.gunshi.project.xyt.controller; import com.baomidou.mybatisplus.core.conditions.query.QueryWrapper; +import com.baomidou.mybatisplus.core.conditions.update.UpdateWrapper; +import com.baomidou.mybatisplus.core.toolkit.CollectionUtils; import com.baomidou.mybatisplus.core.toolkit.IdWorker; import com.baomidou.mybatisplus.core.toolkit.ObjectUtils; import com.baomidou.mybatisplus.core.toolkit.StringUtils; import com.baomidou.mybatisplus.extension.plugins.pagination.Page; import com.gunshi.core.result.R; +import com.gunshi.project.xyt.model.ForecastU; import com.gunshi.project.xyt.model.ForecastUseparam; +import com.gunshi.project.xyt.service.ForecastUService; import com.gunshi.project.xyt.service.ForecastUseparamService; import com.gunshi.project.xyt.validate.markers.Insert; import com.gunshi.project.xyt.validate.markers.Update; @@ -23,9 +27,13 @@ import org.springframework.web.bind.annotation.RequestMapping; import org.springframework.web.bind.annotation.RestController; import java.io.Serializable; +import java.math.BigDecimal; +import java.util.ArrayList; import java.util.Date; import java.util.List; +import java.util.Map; import java.util.Objects; +import java.util.stream.Collectors; /** * 描述: 预报_通用参数管理 @@ -40,6 +48,9 @@ public class ForecastUseparamController { @Autowired private ForecastUseparamService service; + @Autowired + private ForecastUService forecastUService; + @Operation(summary = "新增") @PostMapping("/insert") @@ -59,6 +70,37 @@ public class ForecastUseparamController { dto.setChtm(oldData.getChtm()); dto.setUpdateTime(new Date()); boolean result = service.updateById(dto); + // 修改计算单位线的参数时，需要重新计算一次 + if (result && (dto.getParamCode().equals("cymj") || dto.getParamCode().equals("hdc") || dto.getParamCode().equals("hdpj") || dto.getParamCode().equals("dt") || dto.getParamCode().equals("h")) && !dto.getParamValue().equals(oldData.getParamValue())) { + // 承雨面积areaF："cymj", 河道长lengthL："hdc", 河道坡降j："hdpj", 时段∆t："dt", h毫米净雨："h" + QueryWrapper qw = new QueryWrapper().in("param_code", "cymj", "hdc", "hdpj", "dt", "h"); + List> uMaps = service.listMaps(qw); + Map uMap = uMaps.stream().collect(Collectors.toMap(map -> (String) map.get("param_code"), map -> (String) map.get("param_value"))); + if (uMap.containsKey("cymj") && StringUtils.isNotEmpty(uMap.get("cymj")) + && uMap.containsKey("hdc") && StringUtils.isNotEmpty(uMap.get("hdc")) + && uMap.containsKey("hdpj") && StringUtils.isNotEmpty(uMap.get("hdpj")) + && uMap.containsKey("dt") && StringUtils.isNotEmpty(uMap.get("dt"))) { + List uParam = service.calcU(new BigDecimal(uMap.get("cymj")), + new BigDecimal(uMap.get("hdc")), + new BigDecimal(uMap.get("hdpj")), + Double.valueOf(uMap.get("dt")), + StringUtils.isNotEmpty(uMap.get("h")) ? new BigDecimal(uMap.get("h")) : BigDecimal.ONE, StringUtils.isNotEmpty(uMap.get("swfq")) ? uMap.get("swfq") : "I"); + if (CollectionUtils.isNotEmpty(uParam)) { + forecastUService.remove(new UpdateWrapper<>()); + List uList = new ArrayList<>(); + Date date = new Date(); + for (BigDecimal u : uParam) { + ForecastU forecastU = new ForecastU(); + forecastU.setId(IdWorker.getId()); + forecastU.setUValue(u); + forecastU.setChtm(date); + forecastU.setUpdateTime(date); + uList.add(forecastU); + } + forecastUService.saveBatch(uList); + } + } + } return R.ok(result ? dto : null); } diff --git a/src/main/java/com/gunshi/project/xyt/service/ForecastResultsService.java b/src/main/java/com/gunshi/project/xyt/service/ForecastResultsService.java index a1ac830..30f6e7f 100644 --- a/src/main/java/com/gunshi/project/xyt/service/ForecastResultsService.java +++ b/src/main/java/com/gunshi/project/xyt/service/ForecastResultsService.java @@ -761,106 +761,4 @@ public class ForecastResultsService extends ServiceImpl - * @auther: cxw - * @date: 2024-11-06, 周三, 10:16:26 - */ - private List calcU(BigDecimal areaF, BigDecimal lengthL, BigDecimal j, double dt, BigDecimal tl){ - // 计算m1、n - Map m1AndNMap = calcM1AndN(areaF, lengthL, j); - BigDecimal m1 = m1AndNMap.get("m1"); - BigDecimal n = m1AndNMap.get("n"); - // k = m1 / n。如果m1有非线性改正，则为改正后的值。算例保留1位小数 - BigDecimal k = m1.divide(n, 1, BigDecimal.ROUND_HALF_UP); - /* - * 瞬时单位线u(0,t)=(1/(k*Γ(n)))*((t/k)^(n-1))*e^(-t/k) - */ - // 计算伽马函数 - BigDecimal nGamma = BigDecimal.ONE.divide(k.multiply(BigdecimalUtil.bigDecimalGamma(n)), 4, BigDecimal.ROUND_HALF_UP); - List uList = new ArrayList<>(); - for (int i = 1; i <= 11; i++) { - BigDecimal t = BigDecimal.valueOf(dt * i); - if (t.compareTo(tl) > 0) { - break; - } - // 计算t、k、n - BigDecimal tkn = BigdecimalUtil.bigDecimalExponentiation(t.divide(k, 4, BigDecimal.ROUND_HALF_UP), n.subtract(BigDecimal.ONE).doubleValue()).setScale(2, BigDecimal.ROUND_HALF_UP); - // 计算 e 的指数 - BigDecimal eToThePower = BigdecimalUtil.bigDecimalExponentiation(BigDecimal.valueOf(Math.E), t.divide(k, 4, BigDecimal.ROUND_HALF_UP).multiply(BigDecimal.ONE.negate()).setScale(2, BigDecimal.ROUND_HALF_UP).doubleValue()); - // 计算u(0,t) - BigDecimal u0t = nGamma.multiply(tkn).multiply(eToThePower).setScale(2, BigDecimal.ROUND_HALF_UP); - uList.add(u0t); - } - System.out.println(uList); - // TODO S(t)=0∫t*u(0,t)dt为参数n、k的函数，有S(t)查算表？ - - // TODO 时段单位线公式u(∆t,t)=[S(t)-S(t-∆t)] - - // TODO h毫米净雨时段单位线q(∆t,t)=Fh/(3.6*∆t)*[u(∆t,t)] - - - return new ArrayList<>(); - } - - - /** - * @description: 计算m1和n参数 - * @param areaF - * @param lengthL - * @param j - * @return: void - * @auther: cxw - * @date: 2024-11-06, 周三, 13:16:24 - */ - private Map calcM1AndN(BigDecimal areaF, BigDecimal lengthL, BigDecimal j) { - BigDecimal m1 = BigDecimal.ZERO; - BigDecimal n = BigDecimal.ZERO; - // 所属片区（全省山丘区瞬时单位线分三个片进行参数的地区综合，第一片包括水文分区 1、2、4区即京广线两侧及鄂东黄冈、咸宁地区一带,(江汉平原湖区在外);第二片包括水文分区6、8、9、11区即鄂北，鄂西北及宜昌地区长江以北一带;第三片包括7、10水文分区即清江流域，恩施地区） - String slice = "I"; - if (slice.toUpperCase().equals("I")) {// I片（1、2、4区） - if (areaF.compareTo(new BigDecimal("30")) > 0) {// areaF > 30 - // m1 = 0.82*F^0.29*L^0.23*j^-0.20 - m1 = new BigDecimal("0.82").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.29)).multiply(BigdecimalUtil.bigDecimalExponentiation(lengthL, 0.23)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, -0.20)); - } else if (areaF.compareTo(new BigDecimal("30")) <= 0) {// areaF <= 30 - // m1 = 1.38*F^0.27*L^0.216*j^-0.185 - m1 = new BigDecimal("1.38").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.27)).multiply(BigdecimalUtil.bigDecimalExponentiation(lengthL, 0.216)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, -0.185)); - } - if (j.compareTo(new BigDecimal("5")) > 0) {// j > 5 - // n = 0.34*F^0.35*j^0.1 - n = new BigDecimal("0.34").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.35)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, 0.1)); - } else if (j.compareTo(new BigDecimal("5")) <= 0) {// j <= 5 - // n = 1.04*F^0.3*j^0.1 - n = new BigDecimal("1.04").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.3)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, 0.1)); - } - } else if (slice.toUpperCase().equals("II")) {// II片（6、8、9、11区） - // m1 = 1.64*F^0.231*L^0.131*j^-0.08 - m1 = new BigDecimal("1.64").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.231)).multiply(BigdecimalUtil.bigDecimalExponentiation(lengthL, 0.131)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, -0.08)); - // n = 0.529*F^0.25*j^0.20 - n = new BigDecimal("0.529").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.25)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, 0.20)); - } else if (slice.toUpperCase().equals("III")) {// III片（7、10区） - // m1 = 0.8*F^0.3*L^0.1*j^-0.06 - m1 = new BigDecimal("0.8").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.3)).multiply(BigdecimalUtil.bigDecimalExponentiation(lengthL, 0.1)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, -0.06)); - // n = 0.69*F^0.224*j^0.092 - n = new BigDecimal("0.69").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.224)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, 0.092)); - } - // TODO 特殊流域适用的公式 - - // TODO 当设计雨强超过10mm/h时，除熔岩地区外一般应作参数m1外延的非线性改正 - - BigDecimal finalM1 = m1.setScale(1, BigDecimal.ROUND_HALF_UP); - BigDecimal finalN = n.setScale(1, BigDecimal.ROUND_HALF_UP); - return new HashMap<>() {{ - put("m1", finalM1); - put("n", finalN); - }}; - } } diff --git a/src/main/java/com/gunshi/project/xyt/service/ForecastUseparamService.java b/src/main/java/com/gunshi/project/xyt/service/ForecastUseparamService.java index 22a3b6e..f4cdbbf 100644 --- a/src/main/java/com/gunshi/project/xyt/service/ForecastUseparamService.java +++ b/src/main/java/com/gunshi/project/xyt/service/ForecastUseparamService.java @@ -3,10 +3,18 @@ package com.gunshi.project.xyt.service; import com.baomidou.mybatisplus.extension.service.impl.ServiceImpl; import com.gunshi.project.xyt.mapper.ForecastUseparamMapper; import com.gunshi.project.xyt.model.ForecastUseparam; +import com.gunshi.project.xyt.util.BigdecimalUtil; import lombok.extern.slf4j.Slf4j; import org.springframework.stereotype.Service; import org.springframework.transaction.annotation.Transactional; +import java.math.BigDecimal; +import java.util.ArrayList; +import java.util.HashMap; +import java.util.List; +import java.util.Map; +import java.util.stream.Collectors; + /** * 描述: 预报_通用参数管理 * author: cxw @@ -18,6 +26,122 @@ import org.springframework.transaction.annotation.Transactional; public class ForecastUseparamService extends ServiceImpl { + /** + * @description: 计算时间单位线u参数 + * @param areaF 流域面积 k㎡（承雨面积） + * @param lengthL 主河道长度 km + * @param j 主河道加权平均比降 千分率‰（河道坡降） + * @param dt 时段∆t，时间间隔，例如0.5 + * @param h h毫米净雨，例如1、10等，默认为1 + * @param slice 水文分区 + * @return: java.util.List + * @auther: cxw + * @date: 2024-11-06, 周三, 10:16:26 + */ + public List calcU(BigDecimal areaF, BigDecimal lengthL, BigDecimal j, double dt, BigDecimal h, String slice) { + // 计算m1、n + Map m1AndNMap = calcM1AndN(areaF, lengthL, j, slice); + BigDecimal m1 = m1AndNMap.get("m1"); + BigDecimal n = m1AndNMap.get("n"); + // k = m1 / n。如果m1有非线性改正，则为改正后的值。算例保留1位小数 + BigDecimal k = m1.divide(n, 3, BigDecimal.ROUND_HALF_UP); + /* + * 瞬时单位线u(0,t)=(1/(k*Γ(n)))*((t/k)^(n-1))*e^(-t/k) + */ + List> mapList = new ArrayList<>(); + Map initMap = new HashMap<>(); + initMap.put("t", BigDecimal.ZERO); + initMap.put("tDevideK", BigDecimal.ZERO); + initMap.put("st", BigDecimal.ZERO); + initMap.put("utt", BigDecimal.ZERO); + initMap.put("qtt", BigDecimal.ZERO); + mapList.add(initMap); + for (int i = 1; i <= 11; i++) { + Map result = new HashMap<>(); + // 表1-5地表径流过程计算表：t + BigDecimal t = BigDecimal.valueOf(dt * i); + result.put("t", t); + + // 表1-5地表径流过程计算表：t/k + BigDecimal x = t.divide(k, 10, BigDecimal.ROUND_HALF_UP);// x + result.put("tDevideK", x); + + BigDecimal alpha = n;// α + BigDecimal beta = BigDecimal.ONE;// β + // 表1-5地表径流过程计算表：s(t) + BigDecimal gammaDist = BigDecimal.valueOf(BigdecimalUtil.gammaDist(x.doubleValue(), alpha.doubleValue(), beta.doubleValue(), true)).setScale(10, BigDecimal.ROUND_HALF_UP);// 伽马分布 + result.put("st", gammaDist); + // 表1-5地表径流过程计算表：s(t-∆t)。stt为上一行的st + BigDecimal stt = BigDecimal.ZERO; + stt = mapList.get(i - 1).get("st").setScale(10, BigDecimal.ROUND_HALF_UP); + result.put("stt", stt); + + // 表1-5地表径流过程计算表：u(∆t,t)。utt=st-stt + BigDecimal utt = gammaDist.subtract(stt).setScale(10, BigDecimal.ROUND_HALF_UP); + result.put("utt", utt); + + // 表1-5地表径流过程计算表：q(∆t,t)。qtt=Fh*utt/(3.6*∆t) + BigDecimal qtt = areaF.multiply(h).multiply(utt).divide(new BigDecimal(3.6).multiply(t), 3, BigDecimal.ROUND_HALF_UP); + result.put("qtt", qtt); + mapList.add(result); + } + mapList.remove(0);// 去除initMap + List uList = mapList.stream().map(map -> map.get("qtt")).collect(Collectors.toList()); + return uList; + } + + + /** + * @description: 计算m1和n参数 + * @param areaF + * @param lengthL + * @param j + * @param slice + * @return: void + * @auther: cxw + * @date: 2024-11-06, 周三, 13:16:24 + */ + private Map calcM1AndN(BigDecimal areaF, BigDecimal lengthL, BigDecimal j, String slice) { + BigDecimal m1 = BigDecimal.ZERO; + BigDecimal n = BigDecimal.ZERO; + // 所属片区（全省山丘区瞬时单位线分三个片进行参数的地区综合，第一片包括水文分区 1、2、4区即京广线两侧及鄂东黄冈、咸宁地区一带,(江汉平原湖区在外);第二片包括水文分区6、8、9、11区即鄂北，鄂西北及宜昌地区长江以北一带;第三片包括7、10水文分区即清江流域，恩施地区） + if (slice.toUpperCase().equals("I")) {// I片（1、2、4区） + if (areaF.compareTo(new BigDecimal("30")) > 0) {// areaF > 30 + // m1 = 0.82*F^0.29*L^0.23*j^-0.20 + m1 = new BigDecimal("0.82").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.29)).multiply(BigdecimalUtil.bigDecimalExponentiation(lengthL, 0.23)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, -0.20)); + } else if (areaF.compareTo(new BigDecimal("30")) <= 0) {// areaF <= 30 + // m1 = 1.38*F^0.27*L^0.216*j^-0.185 + m1 = new BigDecimal("1.38").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.27)).multiply(BigdecimalUtil.bigDecimalExponentiation(lengthL, 0.216)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, -0.185)); + } + if (j.compareTo(new BigDecimal("5")) > 0) {// j > 5 + // n = 0.34*F^0.35*j^0.1 + n = new BigDecimal("0.34").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.35)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, 0.1)); + } else if (j.compareTo(new BigDecimal("5")) <= 0) {// j <= 5 + // n = 1.04*F^0.3*j^0.1 + n = new BigDecimal("1.04").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.3)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, 0.1)); + } + } else if (slice.toUpperCase().equals("II")) {// II片（6、8、9、11区） + // m1 = 1.64*F^0.231*L^0.131*j^-0.08 + m1 = new BigDecimal("1.64").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.231)).multiply(BigdecimalUtil.bigDecimalExponentiation(lengthL, 0.131)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, -0.08)); + // n = 0.529*F^0.25*j^0.20 + n = new BigDecimal("0.529").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.25)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, 0.20)); + } else if (slice.toUpperCase().equals("III")) {// III片（7、10区） + // m1 = 0.8*F^0.3*L^0.1*j^-0.06 + m1 = new BigDecimal("0.8").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.3)).multiply(BigdecimalUtil.bigDecimalExponentiation(lengthL, 0.1)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, -0.06)); + // n = 0.69*F^0.224*j^0.092 + n = new BigDecimal("0.69").multiply(BigdecimalUtil.bigDecimalExponentiation(areaF, 0.224)).multiply(BigdecimalUtil.bigDecimalExponentiation(j, 0.092)); + } + // TODO 特殊流域适用的公式 + + // TODO 当设计雨强超过10mm/h时，除熔岩地区外一般应作参数m1外延的非线性改正 + + BigDecimal finalM1 = m1.setScale(3, BigDecimal.ROUND_HALF_UP); + BigDecimal finalN = n.setScale(3, BigDecimal.ROUND_HALF_UP); + return new HashMap<>() {{ + put("m1", finalM1); + put("n", finalN); + }}; + } } diff --git a/src/main/java/com/gunshi/project/xyt/util/BigdecimalUtil.java b/src/main/java/com/gunshi/project/xyt/util/BigdecimalUtil.java index cf1284f..4d047da 100644 --- a/src/main/java/com/gunshi/project/xyt/util/BigdecimalUtil.java +++ b/src/main/java/com/gunshi/project/xyt/util/BigdecimalUtil.java @@ -1,5 +1,7 @@ package com.gunshi.project.xyt.util; +import org.apache.commons.math3.distribution.GammaDistribution; + import java.math.BigDecimal; import java.math.MathContext; import java.math.RoundingMode; @@ -34,18 +36,12 @@ public class BigdecimalUtil { * @auther: cxw * @date: 2024-11-06, 周三, 11:31:24 */ - public static BigDecimal bigDecimalGamma(BigDecimal x) { - BigDecimal result = BigDecimal.valueOf(gamma(x.doubleValue())); - return result.setScale(x.scale(), RoundingMode.HALF_UP); - } - public static double gamma(double x) { - return Math.exp(logGamma(x)); - } - public static double logGamma(double x) { - double tmp = (x - 0.5) * Math.log(x + 4.5) - (x + 4.5); - double ser = 1.0 + 76.18009173 / (x + 0) - 86.50532033 / (x + 1) - + 24.01409822 / (x + 2) - 1.231739516 / (x + 3) - + 0.00120858003 / (x + 4) - 0.00000536382 / (x + 5); - return tmp + Math.log(ser * Math.sqrt(2 * Math.PI)); + public static double gammaDist(double x, double alpha, double beta, boolean cumulative) { + GammaDistribution gammaDist = new GammaDistribution(alpha, 1/beta); + if (cumulative) { + return gammaDist.cumulativeProbability(x); + } else { + return gammaDist.density(x); + } } } diff --git a/src/main/resources/doc/地表径流过程计算表.jpg b/src/main/resources/doc/地表径流过程计算表.jpg new file mode 100644 index 0000000..bb6814c Binary files /dev/null and b/src/main/resources/doc/地表径流过程计算表.jpg differ diff --git a/src/main/resources/doc/时段单位线计算方法.xlsx b/src/main/resources/doc/时段单位线计算方法.xlsx new file mode 100644 index 0000000..1ab28df Binary files /dev/null and b/src/main/resources/doc/时段单位线计算方法.xlsx differ